[richpaj]

Just to recap before I make an observation about the code details we have the
following equations for the 2 phase mixture:

rho (mixture density) = alpha1 * rho1 + alpha2 * rho2 ......(1)

Mass continuity for each phase of the form:

ddt( rho1 * alpha1) + div( rho1 * alpha1 * U ) = 0 ...........(2)

which can be rearranged as

ddt( alpha1) + div( alpha1 * U ) = - ( alpha1 / rho1 ) * DDt( rho1 ) ..........(3)

Then, introducing a compressibility so that rho1 = rho0 + psi1 * p yields

ddt( alpha1) + div( alpha1 * U ) = - ( alpha1 * psi1 / rho1) DDt( p ) .......(4)

A similar equation exists for alpha2 so that by addition we obtain

div( U ) = - { ( alpha1 * psi1 / rho1) + ( alpha2 * psi2 / rho2) } DDt( p ) ......(5)

where alpha1 + alpha2 = 1.


So, (4) can be rewritten with a number of differing possibilities for the RHS e.g

RHS = - ( alpha1 * psi1 / rho1) DDt( p ) (a)
RHS = alpha1 * div( U) + alpha1 * dgdt (b)

where dgdt = { (alpha2 * psi2 / rho2) - (alpha1 * psi1 / rho1) } DpDt cf. pEqn.H i.e.

if (nonOrth == nNonOrthCorr)
{
dgdt =
(pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
*(pEqnComp & p);
..
}


(b) appears to be preferred owing to the use of the explicit form for "alpha1 * div(U)"
to balance the convection term in the transport equation (cf. alphaEqn.H).

It remains to treat the term "alpha1 * dgdt" appropriately and this
is handled either implicitly or explicitly
depending on the sign of dgdt. As a result I would anticipate the following:

dgdt < 0 ---> Implicit, SpTerm += dgdt[celli] ( this will appear as fvm::Sp(SpTerm, alpha1)

dgdt > 0 ----> Explicit, SuTerm += dgdt[celli] * alpha1[celli]

It is given that Sp is initialized with 0 and Su with "alpha1 * div(U)" in alphaEqns.H
in which case I'd expect the following:

forAll(dgdt, celli)
{
if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
{
Su[celli] += dgdt[celli]*alpha1[celli];
}
else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
{
Sp[celli] += dgdt[celli];
}
}

Instead we find the following in alphaEqns.H i.e.

forAll(dgdt, celli)
{
if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
{
Sp[celli] -= dgdt[celli]*alpha1[celli];
Su[celli] += dgdt[celli]*alpha1[celli];
}
else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
{
Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]);
}
}

which appears to suggest Sp should have been initialized with "alpha1 * dgdt ", or have I missed something here I wonder?
Regards,

Richard K.

[scttmllr]

Hi Richard,

In your eqn (4), you would want the LHS term to be of the form "U * div(alpha1)" in order for alpha1 to remain bounded. If you subtract an "alpha1 * divU" from each side of (4), you should get what is in the code. Note that this change will make your compressible pressure equation come out the same, but the RHS term of the alpha equation will change. You can replace the divU that we just put on the RHS of (4) via use of the pressure equation.

I hope that makes sense, and helps.

-Scott

[richpaj]

Hello there,

My understanding though is that discretization is applied to
the conserved form of the phase fraction equation for alpha1
in order to permit a conserved representation of interface compression
(not indicated in eqn. (4) but contained in code). Boundedness is then sought by sorting out the RHS
of (4) into explicit and implict terms.

Futhermore, the structure of the eqn being solved is given by MULESTemplates.C

fvScalarMatrix psiConvectionDiffusion
(
fvm::ddt(rho, psi)
+ fv::gaussConvectionScheme<scalar>(mesh, phi, UDs).fvmDiv(phi, psi)
//- fv::gaussLaplacianScheme<scalar, scalar>(mesh, CDs, snGrads)
//.fvmLaplacian(Dpsif, psi)
- fvm::Sp(Sp, psi)
- Su
);

where here psi<=>alpha1 with the flux terms generated according to alphaEqns.H

Eventually, I intend to compare (A) (original code) and (B) below to see
if they yield similar results. The related experiments we've done
might also possibly guide the more suitable choice of representation,



(A)

Initialization
Su = alpha1 * div( U)
Sp = 0

forAll(dgdt, celli)
{
if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
{
Sp[celli] -= dgdt[celli]*alpha1[celli];
Su[celli] += dgdt[celli]*alpha1[celli];
}
else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
{
Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]);
}
}


(B)

Initialization
Su = alpha1 * div( U) + alpha1 * dgdt
Sp = 0

forAll(dgdt, celli)
{

if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
{
Sp[celli] += dgdt[celli];
Su[celli] -= dgdt[celli]*alpha1[celli];
}
}

Regards,

Richard K.

[richpaj]

I see now I made a mistake in the identification of dgdt it should be:

dgdt = { ( psi2 / rho2) - (psi1 / rho1) } DpDt

The RHS of (4) can now be written as:

RHS = - ( alpha1 * psi1 / rho1) DDt( p )
= - alpha1 * ( div( U) + psi1 / rho1) DDt( p ) ) + alpha1 * div(U)
= alpha1 * alpha2 * dgdt + alpha1 * div( U)


and where the following are noted:

1) treat "alpha1 * div (U)" explicitly to balance the explicit form of the transport equation.
(used to initialize Su)


2) for "dgdt" term

dgdt < 0 ---> Implicit, Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli])
( this will appear as fvm::Sp(SpTerm, alpha1)

dgdt > 0 ----> Explicit, SuTerm += dgdt[celli] * alpha1[celli],
Implicit, SpTerm -= dgdt[celli] * alpha1[celli]

Regards,

RGK

[linch]

Sorry for stupid questions, but:

what's the difference between DDt( p ) and DpDt? Both seem to be Lagrangian derivatives of pressure...

why is

- ( alpha1 * psi1 / rho1) DDt( p ) = - alpha1 * ( div( U) + psi1 / rho1) DDt( p ) ) + alpha1 * div(U)

alpha1 * div(U) was added and alpha1 * div(U)DDt(p) substracted

and from the next step it comes out that:
-(div( U) + psi1 / rho1) DDt( p ) = alpha2 * { ( psi2 / rho2) - (psi1 / rho1) } DpDt
could somebody please explain it in detail?

and why there must be a differentiation between dgdt>0 and dgdt<0?

I'll appreciate any help to clarify the derivation.

Regards,
Ilya

[richpaj]

>what's the difference between DDt( p ) and DpDt? Both seem to be Lagrangian >derivatives of pressure...

DDt is the convective derivative and Ddt is the lagrangian temporal derivative.

>alpha1 * div(U) was added and alpha1 * div(U)DDt(p) substracted

one part will be treated explicitly and the other implicitly

>and why there must be a differentiation between dgdt>0 and dgdt<0?

take a look at P37 of the programmer's manual which explains
how explicit and implicit source terms work.

Good luck,

RGK

[linch]

Thank you, I'll follow your advices!

Originally Posted by richpaj

>what's the difference between DDt( p ) and DpDt? Both seem to be Lagrangian >derivatives of pressure...

DDt is the convective derivative and Ddt is the lagrangian temporal derivative.

Isn't it the same?

best,
Ilya

[nimasam]

hi friends thanks for ur all great discussion!!!
but i still cant understand from where dgdt comes from?
LHS of eq 4 is : - ( alpha1 * psi1 / rho1) DDt( p ) but how we can add alpha*div(U) and subtract alpha*div(U)*ddt(p) ? it would not remain the same!!!!
im totally confused please explain more

[richpaj]

>>>DDt is the convective derivative and Ddt is the lagrangian temporal derivative.

>Isn't it the same?

>>best,
>>Ilya

yes, you're quite right, I should've written "DDt is the convective derivative and Ddt the Eulerian time derivative"

Originally Posted by nimasam

hi friends thanks for ur all great discussion!!!
but i still cant understand from where dgdt comes from?
LHS of eq 4 is : - ( alpha1 * psi1 / rho1) DDt( p ) but how we can add alpha*div(U) and subtract alpha*div(U)*ddt(p) ? it would not remain the same!!!!
im totally confused please explain more

Sorry, there was a "hanging" right ")", the following in (4) namely,

- alpha1 * ( div( U) + psi1 / rho1) DDt( p ) ) + alpha1 * div(U)

should read:

- alpha1 * ( div( U) + psi1 DDt( p ) / rho1 ) + alpha1 * div(U)


hopefully that clarifies.

Rgds.

[reza_65]

Originally Posted by haghajani

Hi there,

I am trying to add energy eq. to compressibleInterFoam.
To do this, enthalpy should be evaluated using T field. Getting inspiration from compressible solvers, thermoType and mixture properties should be given in thermophysicalProperites dictionary; I have no idea how could it be done for multiphase solvers?

I tried to define the enthalpy in createField, i.e. h=Cp.T, but then it throws error; I think it would need default boundary condition, because hEqn is being solved.

Any idea or anyother approach is appreciated?

best regards,
Hamed

...................


i have same problem, too. i have no idea what should i do, for adding energy equation in CompressibleInterFoam solver.
do you solve this problem???
please help me..

with many regards.

Reza khodadadi.

[nimasam]

hi friends
could you tell me why in the pEqn.H in compressibleInterFoam for calculation of phi , only incompressible part is considered?

phi += pEqnIncomp.flux();

[linch]

Originally Posted by richpaj

The RHS of (4) can now be written as:

RHS =
....
= alpha1 * alpha2 * dgdt + alpha1 * div( U)

and where the following are noted:

2) for "dgdt" term
....
dgdt > 0 ----> Explicit, SuTerm += dgdt[celli] * alpha1[celli],
Implicit, SpTerm -= dgdt[celli] * alpha1[celli]

Explicit would be ----> SuTerm += dgdt[celli] * alpha1[celli] * (1.0 - alpha1[celli]) , or am I wrong? Besides of it, there would be no implicit part if the source would be treated explicitly. I still don't understand, where the if dgdt>0 part comes from.

[richpaj]

Originally Posted by linch

Explicit would be ----> SuTerm += dgdt[celli] * alpha1[celli] * (1.0 - alpha1[celli]) , or am I wrong?

yes, what you write looks correct, must've composed the original in
haste or late at night (!)

Originally Posted by linch

I still don't understand, where the if dgdt>0 part comes from.

dgdt is evaluated in pEqn cf.

dgdt =
(pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
*(p_rghEqnComp & p_rgh);

and where the sign variation derives from the coefficients multiplying
'p_rgh' in p_rghEqnComp. If you examine the definition of fvm::Sp(..)
(programmer's manual) you'll see how these terms reinforce the diagonal
coefficients of the fvMatrix.

Rgds,

Richard K.

[linch]

Thanks Richard, but I think we're talking at cross-purposes.

Originally Posted by richpaj

yes, what you write looks correct, must've composed the original in
haste or late at night (!)

In the source codes (compressibleInterFoam) it is indeed

if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
{
    Sp[celli] -= dgdt[celli]*alpha1[celli];
    Su[celli] += dgdt[celli]*alpha1[celli];
}

So it has nothing to do with your haste.

Originally Posted by richpaj

dgdt is evaluated in pEqn cf.

dgdt =
(pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
*(p_rghEqnComp & p_rgh);

and where the sign variation derives from the coefficients multiplying
'p_rgh' in p_rghEqnComp. If you examine the definition of fvm::Sp(..)
(programmer's manual) you'll see how these terms reinforce the diagonal
coefficients of the fvMatrix.

For this one, I wasn't wondering about where dgdt is evaluated. Rather I was wondering about the origin and the meaning of the peace of code quoted above. Sp += foo; Su -= foo; looks "wrong" for me, because Su must be something like Su~Sp*X, where X is the field we're solving for. If alpha1 would have dimensions() different from (0 0 0 0 0 0 0), this wouldn't work. But since alpha1 is dimless, it works.

I hope now you understand my confusion. The quoted code seems to be wrong, but it also seems to perform right. Until this point, I can follow your derivation, but I'm still missing the last step to derive Sp & Su for (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)

Best regards,
Illya

[richpaj]

Originally Posted by linch

Thanks Richard, but I think we're talking at cross-purposes.


In the source codes (compressibleInterFoam) it is indeed

if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
{
    Sp[celli] -= dgdt[celli]*alpha1[celli];
    Su[celli] += dgdt[celli]*alpha1[celli];
}

So it has nothing to do with your haste.


For this one, I wasn't wondering about where dgdt is evaluated. Rather I was wondering about the origin and the meaning of the peace of code quoted above. Sp += foo; Su -= foo; looks "wrong" for me, because Su must be something like Su~Sp*X, where X is the field we're solving for. If alpha1 would have dimensions() different from (0 0 0 0 0 0 0), this wouldn't work. But since alpha1 is dimless, it works.

I hope now you understand my confusion. The quoted code seems to be wrong, but it also seems to perform right. Until this point, I can follow your derivation, but I'm still missing the last step to derive Sp & Su for (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)

Best regards,
Illya

Ok now I understand, I think you appear to have overlooked the factor "alpha1" that will multiply 'Sp'
in fvm::Sp (Sp, alpha1) cf. MULESTemplates.C where the matrix equation to be solved is:

is the following:

fvScalarMatrix psiConvectionDiffusion
(
fvm::ddt(rho, psi)
+ fv::gaussConvectionScheme<scalar>(mesh, phi, UDs).fvmDiv(phi, psi)
//- fv::gaussLaplacianScheme<scalar, scalar>(mesh, CDs, snGrads)
//.fvmLaplacian(Dpsif, psi)
- fvm::Sp(Sp, psi)
- Su
);

Also,

alpha1 * alpha2 * dgdt = ( alpha1 - alpha1^2 ) *dgdt

so that if dgdt> 0 we get

Sp[celli] -= dgdt[celli]*alpha1[celli]; // "-alpha1^2 *dgdt " contribution
Su[celli] += dgdt[celli]*alpha1[celli]; // alpha1 * dgdt contribution

and where we remember to extract the factor alpha1 from Sp, this is inserted later
in the implicit term "fvm::Sp(Sp, alpha1)". Check out the programmer's manual
for reference to fvm::Sp(..).

Good luck,

Regards,

Richard K.

[linch]

Hi Richard,

Originally Posted by richpaj

Ok now I understand, I think you appear to have overlooked the factor "alpha1" that will multiply 'Sp'
...
and where we remember to extract the factor alpha1 from Sp, this is inserted later
in the implicit term "fvm::Sp(Sp, alpha1)". Check out the programmer's manual
for reference to fvm::Sp(..).

No, I haven't overlooked it. Moreover, this fact has initiated my doubts. Because(!) Sp is being multiplied with the primary variable we're solving for, I was expecting it to have different form than Su. So don't worry, I know what Su's & Sp's are

Originally Posted by richpaj

Also,

alpha1 * alpha2 * dgdt = ( alpha1 - alpha1^2 ) *dgdt

so that if dgdt> 0 we get

Sp[celli] -= dgdt[celli]*alpha1[celli]; // "-alpha1^2 *dgdt " contribution
Su[celli] += dgdt[celli]*alpha1[celli]; // alpha1 * dgdt contribution

Yes, that's the answer I was looking for! Thanks a lot. alpha1 * alpha2 is a quadratic term, while ( alpha1 - alpha1^2 ) is linear-quadratic. That was the reason for my confusion.

Thank you very much and have a nice day!

Best regards,
Illya

[LESlie]

Originally Posted by richpaj

2) for "dgdt" term

dgdt < 0 ---> Implicit, Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli])
( this will appear as fvm::Sp(SpTerm, alpha1)

dgdt > 0 ----> Explicit, SuTerm += dgdt[celli] * alpha1[celli],
Implicit, SpTerm -= dgdt[celli] * alpha1[celli]

Regards,

RGK

Thank you for the discussion! Really beneficial!

I have a confusion about the implicit/explicit treatment of the source discretisation.


As:
ddt(alpha1) + U*div(alpha1) = alpha1*alpha2*dgdt
According to the manual when coefficient is greater than zero it is implicit and vice versa.
So shouldn't it be
dgdt > 0 ---> Implicit
dgdt < 0 ---> Explicit
instead?


Cheers,
Leslie

[richpaj]

Originally Posted by LESlie

Thank you for the discussion! Really beneficial!

I have a confusion about the implicit/explicit treatment of the source discretisation.


As:
ddt(alpha1) + U*div(alpha1) = alpha1*alpha2*dgdt
According to the manual when coefficient is greater than zero it is implicit and vice versa.
So shouldn't it be
dgdt > 0 ---> Implicit
dgdt < 0 ---> Explicit
instead?


Cheers,
Leslie

By way of example, in

D[ alpha ]/Dt = rho * alpha

rho < 0 : treated implicitly i.e. reinforces diagonal coefficient of alpha,
rho > 0 : treated explicitly

which informs the handling of "ddt(alpha1) + U*div(alpha1) = alpha1*alpha2*dgdt".

I hope this helps,

Rgds,

Richard.

[LESlie]

Originally Posted by richpaj

By way of example, in

D[ alpha ]/Dt = rho * alpha

rho < 0 : treated implicitly i.e. reinforces diagonal coefficient of alpha,
rho > 0 : treated explicitly

which informs the handling of "ddt(alpha1) + U*div(alpha1) = alpha1*alpha2*dgdt".

I hope this helps,

Rgds,

Richard.

Hi Richard,

This is my confusion, on programmersGuide P-37
"Therefore OpenFOAM provides a mixed source discretisation procedure
that is implicit when the coefficients that are greater than zero, and explicit for the coefficients less than zero."

By the way do you have any idea why they take away the
adjustPhi(phiHbyA, U, p_rgh) in the pEqn.H which for example exists in interFoam? Shouldn't the flux still be adjusted to obey the continuity despite that it's compressible?

Best Regards,
Leslie

[richpaj]

Originally Posted by LESlie

Hi Richard,

This is my confusion, on programmersGuide P-37
"Therefore OpenFOAM provides a mixed source discretisation procedure
that is implicit when the coefficients that are greater than zero, and explicit for the coefficients less than zero."

By the way do you have any idea why they take away the
adjustPhi(phiHbyA, U, p_rgh) in the pEqn.H which for example exists in interFoam? Shouldn't the flux still be adjusted to obey the continuity despite that it's compressible?

Best Regards,
Leslie

Hello there,

1) but the source terms it references are located on the RHS of the given
equation in the following canonical manner

L[U] = SP + SU

(where L is some differential operator and SP&SU are implicit and explicit
source terms for example)

so that when you discretize you get something like

Ap * Up = F(Up-1, .. )

and where Ap is the diagonal coefficient.

2) as for adjustPhi(phiHbyA, U, p_rgh)

this is applied in the case of incompressible problems when the pressure
reference (or level) is unknown. If you look in the source code
(src/finiteVolume/cfdTools/general/adjustPhi/adjustPhi.C) you'll find:

if (p.needReference())
{

}

By contrast, compressible problems have an equation of state or
compressibility relation which acts to constrain the total pressure.

Hopefully this clarifies,

Regards,

Richard