Formulation in compressibleInterFoam
 

Hi,

I am wondering if someone can help me out with the derivation of the implementation used in compressibleInterFoam.

1. Where does the `dgdt' term come from that is used as a source for the alpha equation?

2. How do you get the coefficient that multiplies pEqnComp in pEqn.H?

Thanks!

"dgdt"
 

Hi,

Do you found any answer for your first question. I would be thankfull if you clarify the term "dgdt" in compressibleLesInterFoam.

Thanks

Hi Hamed,

I'm still not positive what it is and where it comes from, but I'm up for a discussion on it.

-STM

pEqn
 

Hi there,
I could work out, the 2nd part, but I am thinking to "dgdt"
http://img718.imageshack.us/img718/1...einterfoam.jpg

Any idea, please share here,

Hamed

Hi,

My derivations look the same as yours. I also do not agree with the coefficient of pEqnComp in the solve() function (your [7]). I think it should just be (1/rho)*d(rho)/d(p) = (alpha1*psi1 + alpha2*psi2)/rho.

As for the dgdt term, I am still missing where it comes from. The units of dgdt are [s^-1], so it appears as if it might somehow be d(alpha)/dt. Recall that in earlier versions of OF alpha was called gamma.

This term comes from the attached equations. In my notation, \tilde{Y} is the liquid mass fraction and \bar{Y} is the liquid volume fraction (which is the same as gamma in interfoam).

Best, Nat

Nat,

Thanks for the reply. In your equations, do you still express density as \rho = \sum_i \bar{Y}_i \rho_i ? If so, I cannot see how the two equations on your second line are consistent. I must be missing something.

-STM

more light pls!
 

Hi Nat,
would you please explain more on your derivation, I didn' get on 2nd line;
Thanks, Hamed

Aha, those equations certainly aren't consistent. :)

Where I wrote:

\tilde{Y} = \frac{\rho}{\rho_L} \bar{Y}

it should've read

\tilde{Y} = \frac{\rho_L}{\rho} \bar{Y}

Also, sorry to post those without an explanation. The first line comes from the continuity equation, where the total derivative of density has been expanded through the chain rule. The second line is an equivalent form of the definition of density for a mixture if you're working with mass fractions instead of volume fractions. If you take the partial of that second line with respect to pressure, you get the third line. \psi is the equivalent isothermal compressibility of the mixture, and \psi_i is the isothermal compressibilities of the individual components.

If \psi_l is a constant (i.e. rho_l = rho_0 + psi*(p-p_0)) and \psi_g is obtained from the ideal gas law -> \psi_g = 1/(R*T), substituting everything back in and converting from mass fraction to volume fraction will give you the fourth line, which looks like whats in compressibleInterFoam.

That was my best guess of what was going on when I looked at the code. I'd love to get that double checked by someone who actually works with that solver a lot.

Hi all,

@Scott:
I have the same feeling, when I look at the definition in createFields.H
So indeed, it seems to be something like the alpha (phase fraction) over the time.

@Nat: I don't want to argue with your solution, I just want to make myself sure I understood everything properly.

I think you meant here 'with respect to density' or do you treat density as a function of pressure?

How did you came up with this, what is your reference for it?? and how did you use it in the 4-th eq.? It is not so obvious to me due to the fact that normally isothermal compressibility is expressed as following
http://upload.wikimedia.org/math/6/0...a4ee1a8c28.png
which means the unit is m^2/N. Furthermore using ideal gas law you can obtain so called compressibility factor Z (which is unit-less). Moreover according to suggestion for psi definition used in transportProperties for both phases the unit for it is s^2/m^2:
and therefore is consistent with the idea: rho = rho0 + psi*p.

What you suggest by your definition of psi_g is mole/J.

Probably I missed something, so I would appreciate, if you could explain that.

Hi Piotr,

Ah yes, I did mean w.r.t. density.

The definition I've got for isothermal compressibility is

\psi = \frac{\partial \rho}{\partial \p} at constant T

So the gas phase gives you

rho = p/RT
psi = 1/RT

I think the confusion is that my R is the specific gas constant (J/kg-k) and yours is the universal gas constant (so there should be a mol/kg in there somewhere).

With regards to a reference, most of this stuff is coming straight from my thesis, which should be finished by next week. If you (or anyone else) want a copy of it, I'd be happy to send it to you if you send me your email address.

-Nat

Ok, then it looks to me quite consistent. Really nice :)

That's funny, cause I do my master thesis now either basing on this solver :) but mine can't contribute (yet) much to the forum due it's in my native language ;/ I will send you my email through PM.

Hi again, sorry for bothering you during Easter (BTW Happy Easter), but I need an answer for a simple (I think), maybe even silly, question and you might know it. So I've read the PhD thesis of Henrik Rusche, I think all of you know it, and found a part about Interface-capturing numerical solution procedure, which is as follows:
http://c7.asteroid.pl/a.eu.fotka.pl....503836_640.jpg
As far as I see compressibleInterFoam solver follows this procedure step-by-step. Is that right?? Can I say this methodology is entirely valid for this solver? Or maybe even a basis for it? Or do you see any exceptions in case of those general rules?

I can't say that for sure, cause the source code still looks to me mysterious sometimes.

//Of course gamma here is our alpha

That's me again :)

I digged deeper and analayzed the pEqn.H code today once again and the main equation looks to me like this (reading from the file):
http://c4.asteroid.pl/a.eu.fotka.pl....765603_640.jpg

instead of the last one here:

Hi,

Yes, it looks like you have it almost correct now. The last term where you have div(U*grad(P)) isn't quite there. The pEqnIncomp terms come from a div(U) in the continuity equation. The U we substitute in here comes from the discretization of the momentum equation, which gives us rUAf* \Delta P. The moral of the story is that it's not `U' as a coefficient now.

-Scott

Actually U* was used as an indicator of the fact that this is not the regular U field any more. This was just asterisk and not the matematical operator :P Sorry, my fault. Anyway thank you.

There is one thing that still bothers me. I mean why there is a difference in the second bracket? Why do I get three components instead of two suggested by Hamed. It doesn't look like I could simply get rid off one of them, in order to get the same form. Or maybe it is the same form, but I can't see it now. Why is that?

The only difference comes from writing:

U \cdot \nabla P = \nabla \cdot (U P) - P \nabla \cdot U

Why is this done? I suppose it is to write it into more of a "standard" form (i.e. div(U*P) ) and give the matrix more diagonal dominance.

Thanks!!! Stupid me of course :P

Energy eq. + Multiphase solver
 

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. :confused:

Any idea or anyother approach is appreciated?

best regards,
Hamed

Have you looked at this http://openfoamwiki.net/index.php/Ho...ure_to_icoFoam ??

dgdt source term in alpha equation
 

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.

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

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.

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

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 * 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

>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

Thank you, I'll follow your advices!
Isn't it the same?

best,
Ilya

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

>>>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"


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.

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


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.

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();

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.

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

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.

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

In the source codes (compressibleInterFoam) it is indeed
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.

Hi Richard,

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 ;)

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

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.

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