# Formulation in compressibleInterFoam

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

 December 10, 2009, 16:51 Formulation in compressibleInterFoam #1 New Member   Scott Miller Join Date: Nov 2009 Posts: 15 Rep Power: 15 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, 13:58 "dgdt" #2 Member   Hamed Aghajani Join Date: Mar 2009 Location: London, UK Posts: 77 Rep Power: 16 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, 18:48 #3 New Member   Scott Miller Join Date: Nov 2009 Posts: 15 Rep Power: 15 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, 06:15
pEqn
#4
Member

Hamed Aghajani
Join Date: Mar 2009
Location: London, UK
Posts: 77
Rep Power: 16
Hi there,
I could work out, the 2nd part, but I am thinking to "dgdt"

Hamed
Attached Images
 compressibleInterFoam.JPG (37.6 KB, 1029 views)

 March 19, 2010, 06:59 #5 New Member   Scott Miller Join Date: Nov 2009 Posts: 15 Rep Power: 15 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, 16:19
#6
New Member

Join Date: Mar 2009
Location: Providence, RI
Posts: 18
Rep Power: 16
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
Attached Images
 pressure.jpg (18.6 KB, 839 views)

 March 25, 2010, 21:30 #7 New Member   Scott Miller Join Date: Nov 2009 Posts: 15 Rep Power: 15 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, 09:05 more light pls! #8 Member   Hamed Aghajani Join Date: Mar 2009 Location: London, UK Posts: 77 Rep Power: 16 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: 16 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.

April 1, 2010, 11:18
#10
Member

Piotr Prusinski
Join Date: Oct 2009
Location: Warsaw, Poland
Posts: 67
Rep Power: 15
Hi all,

@Scott:
Quote:
 Originally Posted by scttmllr 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.
I have the same feeling, when I look at the definition in createFields.H
Code:
volScalarField dgdt =
pos(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001));
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.

Quote:
 Originally Posted by natrask If you take the partial of that second line with respect to pressure, you get the third line.
I think you meant here 'with respect to density' or do you treat density as a function of pressure?

Quote:
 Originally Posted by natrask the ideal gas law -> \psi_g = 1/(R*T)
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

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:
Code:
psi [ 0 -2 2 0 0 ]
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.

 April 1, 2010, 12:55 #11 New Member   Nat Trask Join Date: Mar 2009 Location: Providence, RI Posts: 18 Rep Power: 16 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

April 1, 2010, 16:43
#12
Member

Piotr Prusinski
Join Date: Oct 2009
Location: Warsaw, Poland
Posts: 67
Rep Power: 15
Quote:
 Originally Posted by natrask 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
Ok, then it looks to me quite consistent. Really nice

Quote:
 Originally Posted by natrask With regards to a reference, most of this stuff is coming straight from my thesis,
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: 67 Rep Power: 15 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: 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

April 6, 2010, 15:39
#14
Member

Piotr Prusinski
Join Date: Oct 2009
Location: Warsaw, Poland
Posts: 67
Rep Power: 15
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):

instead of the last one here:
Quote:
 Originally Posted by haghajani

 April 7, 2010, 09:50 #15 New Member   Scott Miller Join Date: Nov 2009 Posts: 15 Rep Power: 15 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: 67 Rep Power: 15 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: 15 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: 67 Rep Power: 15 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: 16 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: 67 Rep Power: 15 Have you looked at this http://openfoamwiki.net/index.php/Ho...ure_to_icoFoam ??

 Tags compressible, compressibleinterfoam, theory