[scttmllr]

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]

Hi,

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

Thanks

[scttmllr]

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]

Hi there,
I could work out, the 2nd part, but I am thinking to "dgdt"


Any idea, please share here,

Hamed

[scttmllr]

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.

[natrask]

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]

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]

Hi Nat,
would you please explain more on your derivation, I didn' get on 2nd line;
Thanks, Hamed

[natrask]

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.

[piprus]

Hi all,

@Scott:

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

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.

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?

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:

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.

[natrask]

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]

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

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.

[piprus]

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

[piprus]

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:

Originally Posted by haghajani

[scttmllr]

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]

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]

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]

Thanks!!! Stupid me of course :P

[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

[piprus]

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