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

About the gammaEqn in interFoam

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

Like Tree22Likes
  • 19 Post By sega
  • 1 Post By suraj
  • 1 Post By isabel
  • 1 Post By sandy

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   January 7, 2009, 04:22
Default Hi guys, I am newbie to Ope
  #1
Senior Member
 
Sandy Lee
Join Date: Mar 2009
Posts: 213
Rep Power: 18
sandy is on a distinguished road
Hi guys,

I am newbie to OpenFOAM. Could somebody explain the code about gammaEqn in interFoam?

surfaceScalarField phic = mag(phi/mesh.magSf());
---------------
phi=uf·Sf
phic=|uf·Sf/Sf|=|uf| ???
---------------

phic = min(interface.cGamma()*phic, max(phic));
---------------
What is "cGamma"? Where can I find its definition?
---------------

surfaceScalarField phir = phic*interface.nHatf();
---------------
nHat=grad(Gamma)/mag(grad(Gamma))
nHatf=fvc::interface(nHat) ????
---------------

for (int gCorr=0; gCorr<nGammaCorr; gCorr++)
{
surfaceScalarField phiGamma =
fvc::flux
(
phi,
gamma,
gammaScheme
)
+ fvc::flux
(
-fvc::flux(-phir, scalar(1) - gamma, gammarScheme),
gamma,
gammarScheme
);
-------------------
What are "fvc::flux" and "gammarScheme"? What is the definition of "phiGamma"?
--------------------

MULES::explicitSolve(gamma, phi, phiGamma, 1, 0);

rhoPhi = phiGamma*(rho1 - rho2) + phi*rho2;
}
-----------------
rhoPhi is mass flux (rho·uf·Sf), right?

In Rusche's thesis, the indicator equation is solved by Eq.(3.58)or (3.59). Why I can not find this equation in gammaEqn.h ??

Please give me some hint.

Thank you in advance.
Sandy
sandy is offline   Reply With Quote

Old   January 7, 2009, 05:54
Default Hi guys, I am newbie to O
  #2
Senior Member
 
sega's Avatar
 
Sebastian Gatzka
Join Date: Mar 2009
Location: Frankfurt, Germany
Posts: 729
Rep Power: 20
sega is on a distinguished road
Quote:
Hi guys,

I am newbie to OpenFOAM. Could somebody explain the code about gammaEqn in interFoam?
I'm currently dealing with the same problem.
Let me try to help you.
All other readers are welcome to correct me and help both of us.

Quote:
surfaceScalarField phic = mag(phi/mesh.magSf());
---------------
phi=uf·Sf
phic=|uf·Sf/Sf|=|uf| ???
---------------
So far so good.

phi=uf & Sf

(It's the dot product of the velocity at the face and the face vector, I will use & as the dot product)

So phic = |phi/|Sf||

Quote:
phic = min(interface.cGamma()*phic, max(phic));
---------------
What is "cGamma"? Where can I find its definition?
---------------
cGamma is a scalar expression for limiting the artificial compression velocity.
You will find it in fvSolution in the PISO sub-dictionary.

This artificial compression velocity is Ur in (3.58) in Rusche's thesis.
So to calculate phir is to calculate phir = Ur & Sf! This is exactly what is done in the following steps!
With cGamma=0 you will simply skip the additional compression, with cGamma=1 the gamma-field will be solved like (3.58) in Rusche!

Quote:
surfaceScalarField phir = phic*interface.nHatf();
---------------
nHat=grad(Gamma)/mag(grad(Gamma))
nHatf=fvc::interface(nHat) ????
---------------
Now this is the phir mentioned above phir = Ur & Sf.
You will find the definitions for nHat and nHatf in the interfaceProperties.C and .H files.

nHatf = nHat & Sf !

nHat = grad(gamma)/(|grad(gamma)| + deltaN)

deltaN is a very small stabilization factor in case |grad(gamma)| will become 0 (the denumerator would be zero without deltaN, which will lead to termination of the program).
This is the case outside the transition region of gamma!

Quote:
for (int gCorr=0; gCorr<nGammaCorr; gCorr++)
{
surfaceScalarField phiGamma =
fvc::flux
(
phi,
gamma,
gammaScheme
)
+ fvc::flux
(
-fvc::flux(-phir, scalar(1) - gamma, gammarScheme),
gamma,
gammarScheme
);
-------------------
What are "fvc::flux" and "gammarScheme"? What is the definition of "phiGamma"?
--------------------
Step by step, I'm not sure what happens in detail here. fvc::flux will calculate explicit values of the phi and phir values defined above.

I understand it like this.
phiGamma = gamma * (uf & Sf) is the transformed (Gauss theorem) div(u*gamma) in (3.58).
When calling it with fvc::flux(phi,gamma,gammaScheme)
it will be discreizised using the gammaScheme, which is the scheme used for div(phi,gamma) in the fvSchemes dictionary. Have a look at the very first line in gammaEqn.H. The entry for div(phi,gamma) is assigned to the variable gammaScheme.

The same holds for calling -fvc::flux(-phir,scalar(1)-gamma,gammarScheme)
which will lead to a discretization of
(1-gamma)*(Ur & Sf) in (3.58) with the fvScheme for div(phirb,gamma). Calling THIS expression like
fvc::flux(THIS,gamma,gammarScheme)
leads to a discretization of
gamma*(ur & Sf)*(1-gamma)

Quote:
MULES::explicitSolve(gamma, phi, phiGamma, 1, 0);
This will give the two above calculated fluxes (the two divergences in (3.58) at the faces to MULES. MULES is a special numerical scheme for solving convective transport equations.
So after calculating the fluxes MULES will solve them explicitly in time.

Quote:
In Rusche's thesis, the indicator equation is solved by Eq.(3.58)or (3.59). Why I can not find this equation in gammaEqn.h ??
Basically the gammaEqn.H simply caluclates the discretisized values of the divergences in (3.58) and gives them to MULES. Maybe a direct representation of (3.58) can be found in MULES.

Quote:
rhoPhi = phiGamma*(rho1 - rho2) + phi*rho2;

-----------------
rhoPhi is mass flux (rho·uf·Sf), right?
I think rhoPhi can be interpreted as an averaged mass flux, weighted between the two phases for gamma? Here I'm guessing more than knowing ...

Dear Sandy.
I hope I could help you and would appreciate and further help from the message board, as this is very interesting for me too.
Greetings. Sebastian.
Bu
__________________
Schrödingers wife: "What did you do to the cat? It's half dead!"
sega is offline   Reply With Quote

Old   June 9, 2009, 15:43
Default Thanks Sebastin!
  #3
New Member
 
Suraj Deshpande
Join Date: Mar 2009
Location: Madison, WI, USA
Posts: 18
Rep Power: 17
suraj is on a distinguished road
Hello Sebastian,
Thanks for such an elaborate explanation. I was looking exactly for this!

Regards,
Suraj
D.R. likes this.
suraj is offline   Reply With Quote

Old   June 30, 2009, 06:22
Default
  #4
Senior Member
 
isabel
Join Date: Apr 2009
Location: Spain
Posts: 171
Rep Power: 16
isabel is on a distinguished road
MULES::explicitSolve(gamma, phi, phiGamma, 1, 0);


This line means this:
gamma: is the actual value to be solved
phi: is the normal convective flux
phiGamma: U*gamma + gamma*(1-gamma)*U
1, 0 : max and min gamma values

Please correct me if I am wrong.
My doubt is: How can I add a source to this equation? I want to solve this:

d(gamma)/dt + div(phigamma) = Source


How can I add the source term?

sharonyue likes this.

Last edited by isabel; June 30, 2009 at 07:49.
isabel is offline   Reply With Quote

Old   July 29, 2009, 08:48
Default
  #5
Member
 
lord_kossity's Avatar
 
Andreas Dietz
Join Date: Mar 2009
Location: Munich
Posts: 79
Rep Power: 17
lord_kossity is on a distinguished road
Hi Isabel,

did you succeed in implementing a source term into the gamma equation? And would you share you're knowledge?

Best,
Andreas
lord_kossity is offline   Reply With Quote

Old   July 31, 2009, 05:01
Default
  #6
Senior Member
 
isabel
Join Date: Apr 2009
Location: Spain
Posts: 171
Rep Power: 16
isabel is on a distinguished road
Hi kossity,

There is a tutorial in Internet called "Solve Cavitating flow around a 2D hydrofoil using a user modified version of interPhaseChangeFoam"
I have followed that tutorial and tried this:

volScalarField Su = source;
volScalarField Sp = 0;
MULES::implicitSolve(oneField(), gamma, phi, phiGamma, Sp, Su, 1, 0);

The problem runs Ok but I don't have good convergence. I don't know if it is because I have set Sp to zero.
isabel is offline   Reply With Quote

Old   July 31, 2009, 09:01
Default
  #7
Senior Member
 
Sandy Lee
Join Date: Mar 2009
Posts: 213
Rep Power: 18
sandy is on a distinguished road
Hi kossity, try to read gammaEqn.H and pdEqn.H of the interPhaseChangeFoam solver, maybe you will find how to add source term to gamma equation.

Hi isable, why you set Sp = 0? Usually it is not useful because it will lead to be disconvergent. People always try to find a way to discrete the source term linearizing (namely Sp*psi + Su) in order to get the "diagonal predominance" of the matrix.
In fact, I also never try to derive this rule , however, it was wrote by every CFD book.
vonboett likes this.
sandy is offline   Reply With Quote

Old   August 3, 2009, 03:55
Default
  #8
Senior Member
 
isabel
Join Date: Apr 2009
Location: Spain
Posts: 171
Rep Power: 16
isabel is on a distinguished road
Hi sandy,

My source is cte*(grad(T) & grad(phi)), so I set:

volVectorField gradT = fvc::grad(T);
volVectorField gradpsi = fvc::grad(psi);
Sp = 0;
Su = cte*(gradT & gradpsi);


I don't know how to discretize Sp in order to became different from zero.
isabel is offline   Reply With Quote

Old   August 3, 2009, 07:57
Default
  #9
Member
 
lord_kossity's Avatar
 
Andreas Dietz
Join Date: Mar 2009
Location: Munich
Posts: 79
Rep Power: 17
lord_kossity is on a distinguished road
Quote:
Originally Posted by sandy View Post
Hi kossity, try to read gammaEqn.H and pdEqn.H of the interPhaseChangeFoam solver, maybe you will find how to add source term to gamma equation.
Thank you - I will do so.


I suppose Isabel is right. Source got to be added in the MULES routine.

Finally added some (senseless) source - works. Now it's time to deal with the real source.

Last edited by lord_kossity; August 3, 2009 at 12:06.
lord_kossity is offline   Reply With Quote

Old   August 3, 2009, 16:40
Default
  #10
Senior Member
 
isabel
Join Date: Apr 2009
Location: Spain
Posts: 171
Rep Power: 16
isabel is on a distinguished road
Hi Kossity,

If you want more information about the gammaEqn, you must real carefully the file MULESTemplates.C
isabel is offline   Reply With Quote

Old   August 3, 2009, 21:52
Default
  #11
Senior Member
 
Sandy Lee
Join Date: Mar 2009
Posts: 213
Rep Power: 18
sandy is on a distinguished road
Quote:
Originally Posted by isabel View Post
Hi sandy,

My source is cte*(grad(T) & grad(phi)), so I set:

volVectorField gradT = fvc::grad(T);
volVectorField gradpsi = fvc::grad(psi);
Sp = 0;
Su = cte*(gradT & gradpsi);

I don't know how to discretize Sp in order to became different from zero.
Hi isabel, my supervisor just told me:

a: if source term = constant, you need to do nothing;

b: if there is a relationship between source term with the variable psi, you can do two ways:

1) if you use the value psi of the last interative step instead of the current varialbe, the source term is still equal to a constant; however, it will lead to a very very slowly iteration, so this method is not available;

2) you can linearize the source term into (Sp*psi + Su). You need to keep Sp <or= 0 in order to get the diagonal predominance.

Do you understand clearly? I seldom meet Sp = 0 because source terms are always a fuction relationship with the variable psi.

PS: I mean psi = the solved variable in a equation.
sandy is offline   Reply With Quote

Old   November 7, 2011, 10:34
Default
  #12
Member
 
bojiezhang
Join Date: Jan 2010
Posts: 64
Rep Power: 16
bojiezhang is on a distinguished road
Quote:
Originally Posted by sandy View Post
Hi isabel, my supervisor just told me:

a: if source term = constant, you need to do nothing;

b: if there is a relationship between source term with the variable psi, you can do two ways:

1) if you use the value psi of the last interative step instead of the current varialbe, the source term is still equal to a constant; however, it will lead to a very very slowly iteration, so this method is not available;

2) you can linearize the source term into (Sp*psi + Su). You need to keep Sp <or= 0 in order to get the diagonal predominance.

Do you understand clearly? I seldom meet Sp = 0 because source terms are always a fuction relationship with the variable psi.

PS: I mean psi = the solved variable in a equation.
Hello sandy:
I have a problem. I want to add a source like source*psi, and the source is the funcion of runtime. When I use the form like
"MULES::explicitSolve(geometricOneField(), alpha1, phi, phiAlpha, source , geometricZeroField() , 1, 0);"

the result comes not true, the source fraction become complex at the souce place and the surface become inconsistent. I do not know what is wrong? can you give me some advice! Thank you!

bojiezhang
bojiezhang is offline   Reply With Quote

Old   February 24, 2012, 03:14
Default different between mDotp() and mDotAlphal()??
  #13
Member
 
vahid
Join Date: Feb 2012
Location: Mashhad-Iran
Posts: 80
Rep Power: 13
vahid.najafi is an unknown quantity at this point
Hello OpenFoam users ,
I'm studying on the Kunz model.we know in this model we have two term for mass dest and prod.(m+ and m-)
What mean mDotAlphal() and mDotP() In interphasechangeFoam?
What is the difference between these two?

vahid.najafi is offline   Reply With Quote

Old   February 21, 2014, 05:31
Default
  #14
New Member
 
Join Date: Nov 2010
Posts: 10
Rep Power: 15
Friederike is on a distinguished road
Hello FOAMers,
this discussion is old but i hope that someone is still interested...
I have a question about the deltaN. Sega said:
Quote:
deltaN is a very small stabilization factor
.
in interfaceProperties.C deltaN is defined as:
1e-8/pow(average(alpha1.mesh().V()), 1.0/3.0)
so it depends on the Volume. With decresing Volume i.e. with a mesh refinement deltaN gets bigger!? Does that make sense? To me it looks like it could become a not so small constant that could influence the simulation?
Has anyone more details about the backround of that formula for deltaN and its influence?
Best regards
Friederike
Friederike is offline   Reply With Quote

Old   March 19, 2014, 10:30
Default
  #15
Member
 
Join Date: Aug 2011
Posts: 89
Rep Power: 14
idefix is on a distinguished road
Hello,

Quote:
in interfaceProperties.C deltaN is defined as:
1e-8/pow(average(alpha1.mesh().V()), 1.0/3.0)
so it depends on the Volume. With decresing Volume i.e. with a mesh refinement deltaN gets bigger!? Does that make sense? To me it looks like it could become a not so small constant that could influence the simulation?
Has anyone more details about the backround of that formula for deltaN and its influence?
I just have the same thoughts about deltaN.
Could anybody comment on this please?
It would be very helpful.

Anyway: I did two calculations
1 ) rough grid: the average cell volume is 10 m³ -> deltaN = 4.6e-9
2 ) fine grid: the average cell volume is 0.0001 m³ -> deltaN = 2.2e-7

So maybe the factor is still small enough not to influence the calculation - but I am not sure about it, because the gradients of alpha are also smaller in finer grid...

Thanks a lot for your help
idefix is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
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 Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Transient interFoam with MRF jaswi OpenFOAM Running, Solving & CFD 18 October 21, 2011 00:29
Nonconservative gammaEqn elisabet OpenFOAM Running, Solving & CFD 20 April 22, 2009 12:50
InterFoam floooo OpenFOAM Running, Solving & CFD 0 November 3, 2008 12:00
About interfoam solver qiu OpenFOAM Running, Solving & CFD 0 May 6, 2007 23:48
GammaEqn and UEqn in interFoam adekian OpenFOAM Running, Solving & CFD 1 April 11, 2007 03:03


All times are GMT -4. The time now is 01:21.