CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM Running, Solving & CFD (
-   -   Nonconservative gammaEqn (

elisabet January 27, 2009 09:53

Hi all OpenFOAMers, I want
Hi all OpenFOAMers,

I want to deal with density currents (water/water with salt) simulations and, to do so, I depart from the twoLiquidMixingFoam and, as a first step, include the RANS libraries (OF-1.5.x). If turbulence has to be taken into account, the gamma equation has to be modified including the turbulent diffusion coefficient, obtained from the turbulent viscosity. The new form of the gammaEqn is then:

fvScalarMatrix gammaEqn (
+ fvm::div(phi, gamma)
- fvm::laplacian(Deff, gamma)
- (fvc::grad(gamma) & fvc::grad(Deff)) //turbulence

where Deff=Dab + turbulence->nut()/Prc; where Prc (or particle Schmidt number) has to be defined by the user).

The surprising thing is that, during the simulations, the total salt concentration (gamma.weightedAverage(mesh.V()).value()) increases, let's say, 5%. However, when I delete from gammaEqn the new term corresponding to variable diffusivity, the total salt mass is conserved.

I tried to modify the numerical schemes for gammaEqn, the default ones are:
grad(gamma) Gauss linear;
grad(Deff) Gauss linear;
div(phi,gamma) Gauss vanLeer;
but OF only allows me to modify the div scheme. I tested also limitedLinear01 but the imbalance was worst, of about 11% of mass increase.

The mesh is structured and orthogonal.

Has anyone any idea of what's going on?


ngj January 27, 2009 16:43

Hi Elisabet I did a similar
Hi Elisabet

I did a similar exercise some time ago, thus I went into my implementation. Apparently the only difference is that I didn't have

- (fvc::grad(gamma) & fvc::grad(Deff))

in my gamma equation, as I took the programmersguide quite literally (correct, I do not know), but I had mass conservation nevertheless.
Adding the above made it becomes non-conserving.

Good luck,


elisabet January 27, 2009 19:01

Hi Niels, Thanks for your a
Hi Niels,

Thanks for your answer.
I cannot omit this term as, when you develop the equations taking into account that Deff is not constant, this new term (or other similar expressions) appears.

Any suggestion?


ngj January 28, 2009 02:31

Just for a clarification, so w
Just for a clarification, so when the programmers guide states, that

laplacian(gamma,U) = Nab dot (gamma Nab U)

then the implementation does only return

gamma laplace U

which is why one have to manually add

Nab U dot Nab gamma?

If that is the case, I have to re-visit the implementation my-self.

Thanks and I will keep my head spinning


elisabet January 28, 2009 09:04

Hi Niels, I guess I'm wrong
Hi Niels,

I guess I'm wrong.
Ok, in momentum equation and for incompressible fluids you have:
div(tau) = div(nu_eff grad(u)) + grad(nu_eff):grad(u)

hence, I just implemented the same procedure for gammaEqn, what I guess it's not correct.
I've found some references for these kinds of flows and no one says anything of this extra term. However, the diffusion coefficient is not constant and takes the form Deff = D_laminar + nu/Sc, where Sc is around 0.7-0.9
If I keep this notation, then, the algorithm is conservative.

Sorry for the confusion and thanks for making me think more about it. Tell me, please, if I'm missing anything.


ngj January 28, 2009 09:26

Hi Elisabet Well, you have
Hi Elisabet

Well, you have made me unsure, thus I went into the source, and at least for fvc:laplacian, see l. 114, it is calculated by considering the diffusion as constant, however I have a hard time understanding whether or not the fvm:laplacian does that as well.

For something quite different:
You say that Sc = 0.7-0.9, but since you are dividing nu_t with this value, you state that the diffusion of the agent is more effective than the diffusion of the water. For turbulent flows it seems counter-intuitive to me?

I considering a small test, as my company needs to clearify the above, thus I will come back on the question of the laplacian.


elisabet January 28, 2009 12:22

Hi Niels, Check in src/fini
Hi Niels,

Check in src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplaci anScheme.C

About the particle Sc number, it is defined as: Sc=nu_t/D_t, where subscript t stands for turbulent. And the same occurs with the thermal Prandtl number in the thermal diffusion coefficient of energy equation. However, in momentum equation your effective viscosity is directly: nu_eff=nu_lam + nu_t. All this comes from the Eddy Viscosity Model (or Boussinesq assumption).


bjj February 2, 2009 08:56

Here are a few thoughts for th
Here are a few thoughts for the above discussion.

After using the version of twoLiquidMixingFoam modified to include the RANS libraries by Niels I have been wondering why the gammeEqn should be modified to include the turbulent diffusivity?

As I see it, when mixing fresh and salt water the diffusion between the two phases is far most dominated by the turbulent spreading. If you imagine a tank with still water and a dense bottom layer and a light top layer I would believe that the mixing of the two phases would be very small.

Therefore when the turbulence model is included this should take care of the spreading/mixing of the two phases through the UEqn. At this point I see the gamma field as a passive scalar used to specify the combined density of the two phases thus any additional diffusion (beside the turbulence diffusion in the UEqn) is not included. If more diffusion is needed i.e. due to molecular diffusion between the two phases a value can be assigned to the Dab factor in the gammaEqn.

Niels, you mention that the only difference is that your version didn't have
If you remember after testing your version we realized that the diffusion between the two phases was significantly too high. A simple lock-exchange test ended up with total mixing of the two phases. Based on the thoughts about where to include the turbulence diffusivity I removed the turbulent diffusivity from the gammaEqn and thereby having the original formulation (setting the Dab factor = 0) and including the turbulence as turbulence->nut() only in the UEqn. This gave a remarkable improvement of the mixing results.

Elisabeth, have you tried to make a test (lock-exchange or similar) with the version where you have deleted the line:
It would be interesting to see if you will experience the same over predicted mixing as we have, as that was what made me think about the issue of including turbulence diffusivity in both gammaEqn and UEqn.


ngj February 11, 2009 16:07

Hi I have revisited the ab

I have revisited the above problem. I have run a lock-flow which measures 2 times 10 m i 2D. The densities of the two liquids are 1000 and 990 respectively.

The initial setup is two identically sized columns which are released. The problem is solved as a test to look at the mixing rate with a k-epsilon model on a structured and uniformly distributed mesh, thus the effect of boundaries is neglected at present.

I have solved it without the additional term discussed above as it should not be there (one can show that the term in the momentum equation shows up because it is a vector field).

Further the eddy viscosity is present in the gamma equations, as it needs to be there to mimic the turbulent mixing.

The results look rather encouraging with a nice lockflow, see figure, where the above is discretized using 112*512 cells and the one below is discretized using 56*256 cells.

So, apparently the mixing seems to be qualitatively correct and I think the chosen way of implementing the mixing of two liquids is the correct way to go.

Best regards,


ngj February 11, 2009 16:12

Hmmm, it didn't work. Though t
Hmmm, it didn't work. Though this does:

Best regards,


elisabet February 17, 2009 05:55

Hi Niels, I did a similar r
Hi Niels,

I did a similar run some weeks ago and I got also nice results. However, we are planning to compare them with experimental data and see exactly how accurate it is. Moreover, I would need to have a correct interface tracking in the top, and I still have to learn more about it. Oh! and don't forget about a new turbulence model that includes the Richardson number! So.... still lots of pendent work...

Best regards,


ngj February 24, 2009 12:03

Hi Elisabet So what you bas
Hi Elisabet

So what you basically need to do is couple this advection-diffusion problem with a VOF-method or the interTrackFoam and let the solution from the advection-diffusion effect the density? However dangerous a word like basic might be

I am not working on it myself, but I now that those who do have had a problem finding experimental data regarding stratified flows, thus do you have any references you wish to share, or are you going to do the experiments by yourself?

Good luck,


P.S. From Milan, I recall you talked about fusion reactors ... isn't a large step to begin looking at stratified flow due to salt and I assume significantly smaller temperatures?

elisabet February 25, 2009 07:02

Hi Niels, As you've said de
Hi Niels,

As you've said density currents is not my field, but I know some people working on it in Barcelona and Rome and I'm trying to collaborate with them, at least in a long term. They are planning to do some experiments and to analyze some experimental results, but still anything is published, sorry. By the way, you remember well: my research field is MHD applied to fusion.

About interTrackFoam I experienced some stability problems for unstructured meshes when the Froude number is very low (in version 1.4.1-dev) so I will repeat the same calculations (not right now) with the new version. Also, the implementation of the Richardson approximation in k-epsilon turbulence model is not straight forward and I should spend more time on it. If you have any experience in one of these two fields, please, do not hesitate sharing it!

Best regards,


ngj February 28, 2009 20:42

Hi Elisabet Just a quick qu
Hi Elisabet

Just a quick question / comment regarding the turbulence modeling:

1. With respect to the momentum equation, you need the uncorrected eddy viscosity (un-corrected meaning no effect of Richardson number).

2. With respect to the gamma-equation, you need the corrected eddy viscosity.

So why is it you want to change the turbulence model and not simply modify the eddy viscosity before applying it to the gamma-equation (i.e. outside the turbulence library), as you need the non-corrected for the momentum?!?
I am not strong in mixing, thus I might be on the wrong track.

By the way, is the modification of the eddy viscosity based on a local Richardson number or a global one? And how do you define the turbulent diffusion coefficient D_t, if Sc is local and thereby unknown? Does it imply the use of an additional equation?

I was at a seminar this Friday, where Håkan Nilsson noted that the standard kOmegaSST does not give reasonable results for swirling flows. The comparison I have made with the same model for mixing does produce to much mixing (to large nu_t I guess in the interior of the domain), thus Håkan introduced a filtered version and achieved significantly better results. It might be interesting with respect to the present problem!?!

Have a nice weekend,


ngj March 1, 2009 13:39

P.S. Send me an email as I hav
P.S. Send me an email as I have some information on the filtered version, I can forward to you.


elisabet March 11, 2009 13:03

Hi Niels, Sorry for the del
Hi Niels,

Sorry for the delay, I was quite busy with other things.

First of all, there are some nowadays public experimental results with density currents. Take a look at:


Michele La Rocca, Claudia Adduce, Giampiero Sciortino, and Allen Bateman Pinzon.Experimental and numerical simulation of three-dimensional gravity currents on smooth and rough bottom. Phys. Fluids 20 106603 (2008)
By the way, related with the fluid equations, I don't agree on the assumption that the turbulence model for N-S equations does not have to be modified (corrected). If you take the full momentum equation, with the Boussinesq hypothesis, and develop the corresponding epsilon and k equations, new terms will appear: a body force production term in k equation (normally called G) and a modified generation-destruction term in the epsilon equation. A good review of possible models could be:


J. Worthy, V. Sanderson, and P. Rubini. A Comparison of Modi ed k-epsilon Turbulence Models
for Buoyant Plumes. School of Mecanical Engineering, Cran eld University, UK.
However, as you know, all these RANS models rely on new experimental constrains that have to be adjusted for each case, hence, it is not an exact science... but DNS is too expensive!



ngj March 18, 2009 03:56

Additional Term in Turbulence Formulation
Hi and good-morning

I have read the second reference by now, and I agree that one need to modify the turbulence formulation. Thanks for shedding some light on it.

I do not fully understand the thermodynamic approximations / relations in the appendix and "T" is left undefined (I am guessing temperature?). However it is merely details. What I am thinking is that one need to decide whether to model turbulence generation due to the density fluctuations or the temperature (?) fluctuations.
This is because I can easily imagine flows where T is constant and the density is varying, and one wants to make the implementation as generic as possible!?

Enjoy your day,


elisabet March 18, 2009 14:17

density fluctuations
Hi Niels,

The 'original' buoyancy term is just (\rho\ \vec{g}), where the density \rho is a function of any other scalar, usually temperature, but it could also be pressure or concentration.
From this point, the Boussinesq hypothesis stands that for moderate temperature (or p or C) variations, the density can be considered as constant in all the terms except for the buoyancy term. What considerably simplifies the momentum equation.
One step further is to choose how the density depends on other variables. The most simple, and common, way is to apply a Taylor expansion, for example:
\rho\ =\ \rho_0\  +\ \left( \frac{\partial \rho}{\partial T} \right)_{p_0}   \left( T-T_0 \right)\ +\ \left( \frac{\partial \rho}{\partial p} \right)_{T_0} \left( p-p_0 \right)\ +\ \frac{1}{2} \left( \frac{\partial^2 \rho}{\partial T^2} \right)_{p_0} \left( T-T_0 \right)^2\ +\  ...

Usually, only the first term of the Taylor expansion is considered, which reduces even more the applicability of the approximation. Moreover, if the compressibility coefficient is considered negligible, the approximation remains:
\rho\ \sim\ \rho_0\ +\ \left( \frac{\partial \rho}{\partial T} \right)_{p_0} \left( T-T_0 \right)
what can be written as:
\rho\ \sim\ \rho_0\ \left( 1\ -\ \beta_0 \left( T-T_0 \right) \right)

and this is what we normally have. However, if the density varies depending on both, the temperature and the concentration, you will obtain an equivalent expression but including a new \beta_C corresponding to the concentration expansion coefficient:
\beta_C=-\frac{1}{\rho_0} \left( \frac{\partial \rho}{\partial C} \right)_0

Sorry if I extended too long but I've seen that very often the Boussinesq hypothesis is not well understood. Hope this will help to someone!

If in your case you have both terms (\beta_T and \beta_C) in your momentum equation, you would have to develop your own k equation and deduce the new terms, it is straight forward as both terms have the same structure. The main drawback, from my point of view, is the definition of each Prandtl number for the Boussinesq approximation of the averaged terms (I'm talking about turbulence now, not about the buoyancy term).



ngj March 20, 2009 05:24

Hi Elisabet

Thank you for the thorough explanation.

Have a nice weekend,


Khelian973 April 21, 2009 05:46

Hi Elisabeth, Niels

The explanation about boussinesq assumption was really interesting... Im a new foamers and im a little confused about the equations...
Indeed i work with buoyantSimpleFoam solver and i though that the boussinesq assumption was considered but i saw on the forum that there is a "special" boussinesq-approximation solver. So I guess this approximation is not used in buoyantSimpleFoam.
If im right, i would like to know which density approximation or density law is used in buoyantSimpleFoam to solve the equations
It's maybe a stupid question...

All times are GMT -4. The time now is 11:51.