CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (http://www.cfd-online.com/Forums/openfoam-solving/)
-   -   Formulation in compressibleInterFoam (http://www.cfd-online.com/Forums/openfoam-solving/70965-formulation-compressibleinterfoam.html)

 scttmllr December 10, 2009 17:51

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!

 haghajani March 18, 2010 14:58

"dgdt"

Hi,

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

Thanks

 scttmllr March 18, 2010 19:48

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

 haghajani March 19, 2010 07:15

pEqn

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

Hamed

 scttmllr March 19, 2010 07:59

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.

1 Attachment(s)
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

 scttmllr March 25, 2010 22:30

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

 haghajani March 26, 2010 10:05

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}

\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.

 piprus April 1, 2010 11:18

Hi all,

@Scott:
Quote:
 Originally Posted by scttmllr (Post 250803) 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 (Post 252064) 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 (Post 252064) 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.

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

 piprus April 1, 2010 16:43

Quote:
 Originally Posted by natrask (Post 252774) 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 (Post 252774) 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.

 piprus April 3, 2010 14:29

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

 piprus April 6, 2010 15:39

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:
Quote:
 Originally Posted by haghajani (Post 250787) http://img718.imageshack.us/img718/1...einterfoam.jpg

 scttmllr April 7, 2010 09:50

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

 piprus April 7, 2010 10:17

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?

 scttmllr April 7, 2010 10:29

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.

 piprus April 7, 2010 10:56

Thanks!!! Stupid me of course :P

 haghajani April 8, 2010 07:03

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

 piprus April 8, 2010 17:02

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

All times are GMT -4. The time now is 07:53.