CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

Formulation in compressibleInterFoam

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

Like Tree9Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   December 10, 2009, 17:51
Default Formulation in compressibleInterFoam
  #1
New Member
 
Scott Miller
Join Date: Nov 2009
Posts: 15
Rep Power: 7
scttmllr is on a distinguished road
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!
scttmllr is offline   Reply With Quote

Old   March 18, 2010, 14:58
Default "dgdt"
  #2
Member
 
Hamed Aghajani
Join Date: Mar 2009
Location: London, UK
Posts: 77
Rep Power: 8
haghajani is on a distinguished road
Hi,

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

Thanks
haghajani is offline   Reply With Quote

Old   March 18, 2010, 19:48
Default
  #3
New Member
 
Scott Miller
Join Date: Nov 2009
Posts: 15
Rep Power: 7
scttmllr is on a distinguished road
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
scttmllr is offline   Reply With Quote

Old   March 19, 2010, 07:15
Exclamation pEqn
  #4
Member
 
Hamed Aghajani
Join Date: Mar 2009
Location: London, UK
Posts: 77
Rep Power: 8
haghajani is on a distinguished road
Hi there,
I could work out, the 2nd part, but I am thinking to "dgdt"


Any idea, please share here,

Hamed
Attached Images
File Type: jpg compressibleInterFoam.JPG (37.6 KB, 212 views)
haghajani is offline   Reply With Quote

Old   March 19, 2010, 07:59
Default
  #5
New Member
 
Scott Miller
Join Date: Nov 2009
Posts: 15
Rep Power: 7
scttmllr is on a distinguished road
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.
scttmllr is offline   Reply With Quote

Old   March 25, 2010, 17:19
Default
  #6
New Member
 
Nat Trask
Join Date: Mar 2009
Location: Providence, RI
Posts: 18
Rep Power: 8
natrask is on a distinguished road
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
File Type: jpg pressure.jpg (18.6 KB, 323 views)
natrask is offline   Reply With Quote

Old   March 25, 2010, 22:30
Default
  #7
New Member
 
Scott Miller
Join Date: Nov 2009
Posts: 15
Rep Power: 7
scttmllr is on a distinguished road
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
scttmllr is offline   Reply With Quote

Old   March 26, 2010, 10:05
Default more light pls!
  #8
Member
 
Hamed Aghajani
Join Date: Mar 2009
Location: London, UK
Posts: 77
Rep Power: 8
haghajani is on a distinguished road
Hi Nat,
would you please explain more on your derivation, I didn' get on 2nd line;
Thanks, Hamed
haghajani is offline   Reply With Quote

Old   March 28, 2010, 20:06
Default
  #9
New Member
 
Nat Trask
Join Date: Mar 2009
Location: Providence, RI
Posts: 18
Rep Power: 8
natrask is on a distinguished road
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.
natrask is offline   Reply With Quote

Old   April 1, 2010, 11:18
Default
  #10
Member
 
Piotr Prusinski
Join Date: Oct 2009
Location: Warsaw, Poland
Posts: 53
Rep Power: 7
piprus is on a distinguished road
Hi all,

@Scott:
Quote:
Originally Posted by scttmllr View Post
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 View Post
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 View Post
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.
piprus is offline   Reply With Quote

Old   April 1, 2010, 12:55
Default
  #11
New Member
 
Nat Trask
Join Date: Mar 2009
Location: Providence, RI
Posts: 18
Rep Power: 8
natrask is on a distinguished road
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
natrask is offline   Reply With Quote

Old   April 1, 2010, 16:43
Default
  #12
Member
 
Piotr Prusinski
Join Date: Oct 2009
Location: Warsaw, Poland
Posts: 53
Rep Power: 7
piprus is on a distinguished road
Quote:
Originally Posted by natrask View Post
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 View Post
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 is offline   Reply With Quote

Old   April 3, 2010, 14:29
Default
  #13
Member
 
Piotr Prusinski
Join Date: Oct 2009
Location: Warsaw, Poland
Posts: 53
Rep Power: 7
piprus is on a distinguished road
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 is offline   Reply With Quote

Old   April 6, 2010, 15:39
Default
  #14
Member
 
Piotr Prusinski
Join Date: Oct 2009
Location: Warsaw, Poland
Posts: 53
Rep Power: 7
piprus is on a distinguished road
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 View Post
piprus is offline   Reply With Quote

Old   April 7, 2010, 09:50
Default
  #15
New Member
 
Scott Miller
Join Date: Nov 2009
Posts: 15
Rep Power: 7
scttmllr is on a distinguished road
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
scttmllr is offline   Reply With Quote

Old   April 7, 2010, 10:17
Default
  #16
Member
 
Piotr Prusinski
Join Date: Oct 2009
Location: Warsaw, Poland
Posts: 53
Rep Power: 7
piprus is on a distinguished road
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.
piprus is offline   Reply With Quote

Old   April 7, 2010, 10:29
Default
  #17
New Member
 
Scott Miller
Join Date: Nov 2009
Posts: 15
Rep Power: 7
scttmllr is on a distinguished road
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.
scttmllr is offline   Reply With Quote

Old   April 7, 2010, 10:56
Default
  #18
Member
 
Piotr Prusinski
Join Date: Oct 2009
Location: Warsaw, Poland
Posts: 53
Rep Power: 7
piprus is on a distinguished road
Thanks!!! Stupid me of course :P
piprus is offline   Reply With Quote

Old   April 8, 2010, 07:03
Default Energy eq. + Multiphase solver
  #19
Member
 
Hamed Aghajani
Join Date: Mar 2009
Location: London, UK
Posts: 77
Rep Power: 8
haghajani is on a distinguished road
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
haghajani is offline   Reply With Quote

Old   April 8, 2010, 17:02
Default
  #20
Member
 
Piotr Prusinski
Join Date: Oct 2009
Location: Warsaw, Poland
Posts: 53
Rep Power: 7
piprus is on a distinguished road
Have you looked at this http://openfoamwiki.net/index.php/Ho...ure_to_icoFoam ??
piprus is offline   Reply With Quote

Reply

Tags
compressible, compressibleinterfoam, theory

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Low Reynolds Number k-epsilon 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


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