# Formulation in compressibleInterFoam

 Register Blogs Members List Search Today's Posts Mark Forums Read

 September 24, 2010, 15:36 dgdt source term in alpha equation #21 Member   Richard Kenny Join Date: Mar 2009 Posts: 59 Rep Power: 9 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. pfhan, sharonyue and RjwV like this.

 September 24, 2010, 15:55 #22 New Member   Scott Miller Join Date: Nov 2009 Posts: 15 Rep Power: 7 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

 September 25, 2010, 08:14 #23 Member   Richard Kenny Join Date: Mar 2009 Posts: 59 Rep Power: 9 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(mesh, phi, UDs).fvmDiv(phi, psi) //- fv::gaussLaplacianScheme(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.

 September 29, 2010, 02:56 #24 Member   Richard Kenny Join Date: Mar 2009 Posts: 59 Rep Power: 9 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 pfhan, RjwV, YHBen and 1 others like this.

January 24, 2011, 12:57
#25
Senior Member

Illya Shevchuk
Join Date: Aug 2009
Posts: 174
Rep Power: 8
Sorry for stupid questions, but:

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

why is
Quote:
 - ( 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

 January 28, 2011, 05:03 #26 Member   Richard Kenny Join Date: Mar 2009 Posts: 59 Rep Power: 9 >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

January 28, 2011, 09:26
#27
Senior Member

Illya Shevchuk
Join Date: Aug 2009
Posts: 174
Rep Power: 8
Quote:
 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

Last edited by linch; February 11, 2011 at 05:20.

 September 1, 2011, 04:20 #28 Senior Member   Nima Sam Join Date: Sep 2009 Location: Tehran, Iran Posts: 1,123 Blog Entries: 1 Rep Power: 14 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

September 1, 2011, 06:11
#29
Member

Richard Kenny
Join Date: Mar 2009
Posts: 59
Rep Power: 9
>>>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"

Quote:
 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)

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

hopefully that clarifies.

Rgds.

September 5, 2011, 14:18
#30
New Member

Join Date: Apr 2011
Location: Iran , tehran
Posts: 6
Rep Power: 6
Quote:
 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???

with many regards.

 September 15, 2011, 16:24 #31 Senior Member   Nima Sam Join Date: Sep 2009 Location: Tehran, Iran Posts: 1,123 Blog Entries: 1 Rep Power: 14 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();

May 21, 2012, 08:13
#32
Senior Member

Illya Shevchuk
Join Date: Aug 2009
Posts: 174
Rep Power: 8
Quote:
 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.

May 21, 2012, 11:01
#33
Member

Richard Kenny
Join Date: Mar 2009
Posts: 59
Rep Power: 9
Quote:
 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 (!)

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

May 21, 2012, 11:46
#34
Senior Member

Illya Shevchuk
Join Date: Aug 2009
Posts: 174
Rep Power: 8
Thanks Richard, but I think we're talking at cross-purposes.

Quote:
 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
Code:
```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.

Quote:
 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

May 21, 2012, 21:27
#35
Member

Richard Kenny
Join Date: Mar 2009
Posts: 59
Rep Power: 9
Quote:
 Originally Posted by linch Thanks Richard, but I think we're talking at cross-purposes. In the source codes (compressibleInterFoam) it is indeed Code: ```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)
//.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.

May 22, 2012, 04:35
#36
Senior Member

Illya Shevchuk
Join Date: Aug 2009
Posts: 174
Rep Power: 8
Hi Richard,

Quote:
 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

Quote:
 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

June 4, 2013, 09:03
#37
New Member

NaiXian Leslie Lu
Join Date: Jun 2009
Location: France
Posts: 26
Rep Power: 8
Quote:
 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

Cheers,
Leslie

June 4, 2013, 21:16
#38
Member

Richard Kenny
Join Date: Mar 2009
Posts: 59
Rep Power: 9
Quote:
 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.

June 5, 2013, 10:09
#39
New Member

NaiXian Leslie Lu
Join Date: Jun 2009
Location: France
Posts: 26
Rep Power: 8
Quote:
 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 coeﬃcients that are greater than zero, and explicit for the coeﬃcients 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

June 5, 2013, 21:15
#40
Member

Richard Kenny
Join Date: Mar 2009
Posts: 59
Rep Power: 9
Quote:
 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 coeﬃcients that are greater than zero, and explicit for the coeﬃcients 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

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

 Tags compressible, compressibleinterfoam, theory

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post Chris CFX 4 December 8, 2009 00:51 Rave Main CFD Forum 0 August 11, 2008 14:55 kulwinder FLUENT 0 May 22, 2004 18:44 Pedro Phoenics 1 July 5, 2001 12:17 Fernando Velasco Hurtado Main CFD Forum 3 January 7, 2000 17:51

All times are GMT -4. The time now is 23:35.