
[Sponsors] 
September 24, 2010, 15:36 
dgdt source term in alpha equation

#21 
Member
Richard Kenny
Join Date: Mar 2009
Posts: 64
Rep Power: 10 
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. 

September 24, 2010, 15:55 

#22 
New Member
Scott Miller
Join Date: Nov 2009
Posts: 15
Rep Power: 9 
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: 64
Rep Power: 10 
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. 

September 29, 2010, 02:56 

#24 
Member
Richard Kenny
Join Date: Mar 2009
Posts: 64
Rep Power: 10 
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 

January 24, 2011, 12:57 

#25  
Senior Member
Illya Shevchuk
Join Date: Aug 2009
Location: Darmstadt, Germany
Posts: 176
Rep Power: 9 
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:
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: 64
Rep Power: 10 
>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
Location: Darmstadt, Germany
Posts: 176
Rep Power: 9 
Thank you, I'll follow your advices!
Quote:
best, Ilya Last edited by linch; February 11, 2011 at 05:20. 

September 1, 2011, 04:20 

#28 
Senior Member
Nima Samkhaniani
Join Date: Sep 2009
Location: Tehran, Iran
Posts: 1,211
Blog Entries: 1
Rep Power: 17 
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: 64
Rep Power: 10 
>>>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:
 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. 

September 5, 2011, 14:18 

#30  
New Member
Reza khodadadi
Join Date: Apr 2011
Location: London, UK
Posts: 18
Rep Power: 7 
Quote:
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. 

September 15, 2011, 16:24 

#31 
Senior Member
Nima Samkhaniani
Join Date: Sep 2009
Location: Tehran, Iran
Posts: 1,211
Blog Entries: 1
Rep Power: 17 
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
Location: Darmstadt, Germany
Posts: 176
Rep Power: 9 
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: 64
Rep Power: 10 
Quote:
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. 

May 21, 2012, 11:46 

#34  
Senior Member
Illya Shevchuk
Join Date: Aug 2009
Location: Darmstadt, Germany
Posts: 176
Rep Power: 9 
Thanks Richard, but I think we're talking at crosspurposes.
Quote:
Code:
if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0) { Sp[celli] = dgdt[celli]*alpha1[celli]; Su[celli] += dgdt[celli]*alpha1[celli]; } Quote:
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: 64
Rep Power: 10 
Quote:
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. 

May 22, 2012, 04:35 

#36  
Senior Member
Illya Shevchuk
Join Date: Aug 2009
Location: Darmstadt, Germany
Posts: 176
Rep Power: 9 
Hi Richard,
Quote:
Quote:
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: 9 
Quote:
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 

June 4, 2013, 21:16 

#38  
Member
Richard Kenny
Join Date: Mar 2009
Posts: 64
Rep Power: 10 
Quote:
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: 9 
Quote:
Hi Richard, This is my confusion, on programmersGuide P37 "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: 64
Rep Power: 10 
Quote:
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(Up1, .. ) 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 

Tags 
compressible, compressibleinterfoam, theory 
Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Low Reynolds Number kepsilon formulation CFX 10.0  Chris  CFX  4  December 8, 2009 00:51 
Immersed Boundary Formulation  Rave  Main CFD Forum  0  August 11, 2008 14:55 
DPM Steady formulation with collisions  kulwinder  FLUENT  0  May 22, 2004 18:44 
energy equation formulation  Pedro  Phoenics  1  July 5, 2001 12:17 
Compressible vs. Incompressible formulations  Fernando Velasco Hurtado  Main CFD Forum  3  January 7, 2000 17:51 