
[Sponsors] 
December 10, 2009, 17:51 
Formulation in compressibleInterFoam

#1 
New Member
Scott Miller
Join Date: Nov 2009
Posts: 15
Rep Power: 9 
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! 

March 18, 2010, 14:58 
"dgdt"

#2 
Member
Hamed Aghajani
Join Date: Mar 2009
Location: London, UK
Posts: 77
Rep Power: 9 
Hi,
Do you found any answer for your first question. I would be thankfull if you clarify the term "dgdt" in compressibleLesInterFoam. Thanks 

March 18, 2010, 19:48 

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

March 19, 2010, 07:15 
pEqn

#4 
Member
Hamed Aghajani
Join Date: Mar 2009
Location: London, UK
Posts: 77
Rep Power: 9 
Hi there,
I could work out, the 2nd part, but I am thinking to "dgdt" Any idea, please share here, Hamed 

March 19, 2010, 07:59 

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

March 25, 2010, 17:19 

#6 
New Member
Nat Trask
Join Date: Mar 2009
Location: Providence, RI
Posts: 18
Rep Power: 9 
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 

March 25, 2010, 22:30 

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

March 26, 2010, 10:05 
more light pls!

#8 
Member
Hamed Aghajani
Join Date: Mar 2009
Location: London, UK
Posts: 77
Rep Power: 9 
Hi Nat,
would you please explain more on your derivation, I didn' get on 2nd line; Thanks, Hamed 

March 28, 2010, 20:06 

#9 
New Member
Nat Trask
Join Date: Mar 2009
Location: Providence, RI
Posts: 18
Rep Power: 9 
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*(pp_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. 

April 1, 2010, 11:18 

#10  
Member
Piotr Prusinski
Join Date: Oct 2009
Location: Warsaw, Poland
Posts: 53
Rep Power: 9 
Hi all,
@Scott: Quote:
Code:
volScalarField dgdt = pos(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001)); @Nat: I don't want to argue with your solution, I just want to make myself sure I understood everything properly. Quote:
How did you came up with this, what is your reference for it?? and how did you use it in the 4th eq.? It is not so obvious to me due to the fact that normally isothermal compressibility is expressed as following which means the unit is m^2/N. Furthermore using ideal gas law you can obtain so called compressibility factor Z (which is unitless). Moreover according to suggestion for psi definition used in transportProperties for both phases the unit for it is s^2/m^2: Code:
psi [ 0 2 2 0 0 ] 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. 

April 1, 2010, 12:55 

#11 
New Member
Nat Trask
Join Date: Mar 2009
Location: Providence, RI
Posts: 18
Rep Power: 9 
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/kgk) 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 

April 1, 2010, 16:43 

#12  
Member
Piotr Prusinski
Join Date: Oct 2009
Location: Warsaw, Poland
Posts: 53
Rep Power: 9 
Quote:
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. 

April 3, 2010, 14:29 

#13 
Member
Piotr Prusinski
Join Date: Oct 2009
Location: Warsaw, Poland
Posts: 53
Rep Power: 9 
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 Interfacecapturing numerical solution procedure, which is as follows:
As far as I see compressibleInterFoam solver follows this procedure stepbystep. 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 

April 6, 2010, 15:39 

#14 
Member
Piotr Prusinski
Join Date: Oct 2009
Location: Warsaw, Poland
Posts: 53
Rep Power: 9 

April 7, 2010, 09:50 

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

April 7, 2010, 10:17 

#16 
Member
Piotr Prusinski
Join Date: Oct 2009
Location: Warsaw, Poland
Posts: 53
Rep Power: 9 
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? Last edited by piprus; April 7, 2010 at 10:53. 

April 7, 2010, 10:29 

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

April 7, 2010, 10:56 

#18 
Member
Piotr Prusinski
Join Date: Oct 2009
Location: Warsaw, Poland
Posts: 53
Rep Power: 9 
Thanks!!! Stupid me of course :P


April 8, 2010, 07:03 
Energy eq. + Multiphase solver

#19 
Member
Hamed Aghajani
Join Date: Mar 2009
Location: London, UK
Posts: 77
Rep Power: 9 
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 

April 8, 2010, 17:02 

#20 
Member
Piotr Prusinski
Join Date: Oct 2009
Location: Warsaw, Poland
Posts: 53
Rep Power: 9 
Have you looked at this http://openfoamwiki.net/index.php/Ho...ure_to_icoFoam ??


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 