View Issue Details

IDProjectCategoryView StatusLast Update
0003310OpenFOAMBugpublic2019-09-18 12:48
Reporteruser136Assigned To 
Status newResolutionopen 
PlatformHPC-ClusterOSCentOSOS Versionunknown
Product Version 
Fixed in Version 
Summary0003310: Variable density not correctly treated in turbulence models?

We're doing large scale simulations with interFoam (water-air, flow around hydraulic structures like weirs). In these I observed an inplausible exchange of the variables of the turbulence models (k, omega, ...) from air to water and vice versa. For real life problems, this leads to huge errors in the turbulent viscosity.

I prepared a small testcase (interFoam, k-OmegaSST) to show the problem:

- Closed box
- Lower half filled with water, low value of k
- Upper half filled with air, high value of k
- No flow

Expected behaviour:

The high values of k in the air should slowly vanish. A certain amount of k should be transferred into the water (diffusion of k), but this is expected to produce only a very small rise of k in the water, due to the large density difference (NOT talking about interface effects, stratification etc.!)

Observed behaviour:

The high values of k in the air diffuse into the water with no visible impact of the different densities. From a physical point of view this is impossible, because it means that the total energy in the system grows significantly. Calculating the volume integral of "(rho*k)" for the whole domain gives the "real" turbulent kinetic energy in the system: It should drop over time, but instead it rises by more than two orders of magnitude. The attached pictures illustrate this.


k (turbulent kinetic energy) is not really an energy, but energy per mass unit. This it is not a quantity which follows a global conservation law. Instead, conservation must be guaranted for "(rho*k)" for advection and diffusion. Looking in the source code, we see that in the k-OmegaSST-model this seems to be done correctly by multiplying rho and k in the PDE:

     fvm::ddt(alpha, rho, k_)
      + ...

But I have the impression, that rho is _not_ the variable rho field as it should be, but a constant (geometricOneField) and thus conservation of energy is not fulfilled. The same problem is observed for other variables of the turbulence models (epsilon, omega) and for several turbulence models when using interFoam.

Best regards,

Carsten Thorenz (user136)
Steps To ReproduceExtract and run the attached test case in OpenFoam-7:

cd _bugreport
blockMesh; setFields; interFoam

TagsNo tags attached.



2019-07-19 11:43

reporter (456,395 bytes)
results.0000.jpg (56,131 bytes)
results.0000.jpg (56,131 bytes)
results.0002.jpg (76,441 bytes)
results.0002.jpg (76,441 bytes)
results.0006.jpg (65,170 bytes)
results.0006.jpg (65,170 bytes)


2019-07-19 11:54

manager   ~0010616

Last edited: 2019-07-19 12:01

View 2 revisions

Currently in interFoam the turbulence model is formulated based on the divergence free volumetric flux for incompressible flow. There is no mass-flux through the interface so k and epsilon are not transferred between the phases by a convective flux (but diffusion across the interface is permitted). Given that the density ratio is generally VERY large solving for rho*k and rho*epsilon is quite stiff and the values in the low density phase are likely to accumulate numerical error.

Have you tried setting-up interFoam with the compressible form turbulence models?


2019-07-19 12:16

manager   ~0010617

compressibleInterFoam is setup for either a mixture or individual phase turbulence models, all in fully conservative compressible form. You could run your case with this solver and incompressible equation of state for each phase and compare with the interFoam results. If the fully conservative approach turbulence proves reliable interFoam can be changed to be consistent with compressibleInterFoam.


2019-07-19 12:22

reporter   ~0010618

About stiffness: Yes, is is numerically hard to take the density jump into account. But from a physical point of view it is necessary to solve for rho*k etc., Furthermore, "d(rho'k)/dt" is what is written in many publications, so the code should reproduce it.

About compressible turbulence models: I'm not sure how to do that. I just now added "libs ("");" kinto controlDict. Is that enough? If yes: It doesn't change the behaviour.


2019-07-19 12:23

reporter   ~0010619

... sorry, didn't read your response in time. I'll try.


2019-07-19 12:29

manager   ~0010620

> About compressible turbulence models: I'm not sure how to do that. I just now added "libs ("");" kinto controlDict. Is that enough?

No, you would need to make some changes to interFoam (see how the turbulence models are instantiated in compressibleInterFoam) but it could be done.


2019-07-19 13:03

reporter   ~0010621

I changed the testcase to be run with compressibleInterFoam instead. Sorry, I'm no expert for "compressible", so maybe the setup is not perfect (setup is attached).

The results are MUCH better: k no longer diffuses horribly into the water (see attached .avi-file) . Still the total energy in the system rises, though much less than with interFoam. This might be a numerical issue, as increasing the grid resolution reduced this issue (though not completely)?

I'll try over the weekend, how one of our "real world" cases reacts on running with compressibleInterFoam.

_bugreport_compressible.avi (138,558 bytes) (1,083,703 bytes)


2019-07-19 16:09

reporter   ~0010624

Hello, in my opinion this is not a bug, but just an expected behaviour of turbulence models at the interface.
I would suggest implementing the interface turbulence damping on k-w (or, if you prefer, k-wSST) model, as described here:

I made many simulations of hydraulic structures (and often we make also experimental validations of the numerical results): in my experience the workhorse for these simulations is the standard k-w with the turbulence damping.
Just my two cents.


2019-07-21 13:41

reporter   ~0010628

Hi, the current implementation of k-omega SST model will generate over-predicted turbulent viscosity (eddy viscosity) near the interface of two-phase flow.

One workaround is introducing buoyancy effect into k-omega SST model, which requires solving for 'd(rho*k)/dt' etc and add a buoyancy generation term for k-equation, rather than the current implementation of solving for 'dk/dt' etc.

Further information can be referred to the following two papers:

[1] Devolder, B., Troch, P., & Rauwoens, P. (2018). Performance of a buoyancy-modified k-ω and k-ω SST turbulence model for simulating wave breaking under regular waves using OpenFOAM®. Coastal Engineering, 138, 49–65.

[2] Larsen, B. E., & Fuhrman, D. R. (2018). On the over-production of turbulence beneath surface waves in Reynolds-averaged Navier–Stokes models. Journal of Fluid Mechanics, 853, 419–460.


2019-07-22 07:59

reporter   ~0010630

I completely agree with michele and wwzhao, that the results of k-Omega SST (and others) would be much better at the interface, if a buoyancy (or stratification) term would be introduced. But in my oppinion, this would be a second step on top of the first:

The implementation should firstly reproduce what was written down by Menter ( And that includes taking the density into account in the PDE as a variable.

Experience from a users point of view: Over the weekend I computed one of our "real world" cases both with compressibleInterFoam and with interFoam with the modified turbulence model by Devolder (but without the buoyancy term): The results are better than with the standard implementation, but not perfect. As michele and wwzhao said this is due to the buoyancy/stratification effects. I can even assume systems where the results would be worse: If you are interested in the air in an air-water-system, the density might be detrimental. But again: I think the model should reproduce what was written by Menter.


2019-07-22 08:35

manager   ~0010631

It is not clear to me that using the physical buoyancy correction is appropriate in this case as the density distribution is a consequence of the numerical approximation and not a representation of a true distribution as it is for example in driftFluxFoam. However it is clear that when using a mixture turbulence model in VoF some correction term is needed to compensate for the unphysical generation of turbulence in the interface region. If this term is included is then it is not clear if formulating the mixture equations in energy conservative form is necessary or a good idea and some research is needed. The alternative would be to use a full two-phase formulation for the turbulence in which separate incompressible form equations are solved for the two phase which is now an option in compressibleInterFoam. The problem with this approach is interface coupling and far-field stabilisation terms will be required so that the turbulence fields are defined in the limit of the phase fraction -> 0. Again this is a research topic.


2019-07-22 11:52

reporter   ~0010632

I think we must take care not to mix up several aspects.

- Buoyancy corrections for k-Omega-SST are necessary (but dirty) tricks to compensate for several problems
- Buoyancy correction is not only needed for the unphysical generation of turbulence, but mainly due to the fact that turbulence is "destroyed" due to stabilization effects at the water surface
- Buoyancy correction is only a dirty helper, because at the water surface the turbulence loses its isotropic properties and becomes partially anisotropic. So, strictly spoken k-Omega-SST is just not able to reproduce that.
- Apart from those considerations I think, that the implementation should reflect the original paper by Menter


2019-07-22 12:18

manager   ~0010634

Last edited: 2019-07-22 12:19

View 2 revisions

- Buoyancy corrections for k-Omega-SST are necessary (but dirty) tricks to compensate for several problems

I agree that some kind of compensation term is need at the interface but it is not clear that this as a "buoyancy" correction; some people use a physical buoyancy model to correct the numerical error but others use a specially formulated numerical correction.

- Buoyancy correction is not only needed for the unphysical generation of turbulence, but mainly due to the fact that turbulence is "destroyed" due to stabilization effects at the water surface

It is not clear to me that a standard buoyancy model is suitable or applicable for this.

- Buoyancy correction is only a dirty helper, because at the water surface the turbulence loses its isotropic properties and becomes partially anisotropic. So, strictly spoken k-Omega-SST is just not able to reproduce that.

It is likely that a two-phase formulation in which separate models are used for each phase would be more appropriate if suitable interface coupling terms can be formulated.

- Apart from those considerations I think, that the implementation should reflect the original paper by Menter

The original paper by Menter is for a low density variation single phase compressible fluid, there is no discussion of the applicability of the model to an two-phase incompressible VoF system in which the density variation is very large and numerical in origin.

Given that in interFoam the fluids are incompressible an incompressible formulation can be used. However in this there is a non-conservative transfer of the energy error at the surface rather than a conservative transfer of the energy error. If interface correction terms are added such that the turbulence goes to 0 at the interface there would be no transfer of energy across the interface and either the compressible or incompressible formulations would produce the same result but the incompressible is preferable for numerical reasons. If interface correction terms are added which do not reduce the turbulence to 0 at the interface then we would also need a explicit model for the potential transfer of this energy across the interface rather than relying on the diffusion terms in the turbulence models doing the right thing. So while I do not have a problem in principle with the turbulence model being changed to conservative form this change cannot be undertaken on the basis of the results of 1 case and without any interface correction term. interFoam has been used by thousands of people for thousands of cases over the many years since I wrote it in the late 90's and we cannot make this kind of fundamental change with knowing for what range of cases is brings benefit. Also more research is needed into the interface correction terms in both the compressible and incompressible formulations so that the best combination for a wide range of cases can be selected. It may be that for some cases having two turbulence models is better and support for this can be added if suitable coupling models are available or can be created.


2019-07-22 12:38

reporter   ~0010635

I see the point in your considerations. And I agree that taking the density into account doesn't heal every aspect and would result in an inconsistency to the old behaviour. What about introducing a density-aware variant of k-OmegaSST as a starter for further enhancements (damping etc.)?


2019-07-22 12:46

manager   ~0010636

Changing interFoam to support both incompressible and compressible forms of turbulence models in a general way would be a substantial task and add a lot to the complexity and maintenance overhead; surely we should first find out if this brings sufficient benefit?
Have you tried the interface correction proposed by @michele? It makes sense to me and has the character of a wall-function for the interface.


2019-07-22 14:58

reporter   ~0010638

I tried the implementaion by Devolder ( That one incorporates density plus a sink term for k. Works much better than the original model, but I can not judge if the sink term produces the "right" quantities. Didn't try the formulation proposed by michele.

The implementation by Devolder is available here:

OF222 - OF6:


2019-07-22 15:31

manager   ~0010639

That is a very odd hacked implementation, the lookup of rho should not be necessary as it is an argument to the model, the mass flux is not correct and non-conservative. It would also make sense that the buoyancy term were implemented as an fvOption rather than hard-coded into the model. My feeling is that it is inappropriate to use a buoyancy model for this correction given that there is nothing physical about the density distribution, it is purely a consequence of the numerics of interface capturing.


2019-09-11 12:14

reporter   ~0010733

For those interested in testing the turbulence damping on k-w (Egorov 2004), i.e. the approach suggested in here attached a derivation of the kOmega with that implementations. The model is tested on OpenFoam 5.x.

kOmegaAlpha.tgz (4,269 bytes)


2019-09-18 12:48

reporter   ~0010762

Please check this thread for a follow-up discussion:

Issue History

Date Modified Username Field Change
2019-07-19 11:43 user136 New Issue
2019-07-19 11:43 user136 File Added:
2019-07-19 11:43 user136 File Added: results.0000.jpg
2019-07-19 11:43 user136 File Added: results.0002.jpg
2019-07-19 11:43 user136 File Added: results.0006.jpg
2019-07-19 11:54 henry Note Added: 0010616
2019-07-19 12:01 henry Note Edited: 0010616 View Revisions
2019-07-19 12:16 henry Note Added: 0010617
2019-07-19 12:22 user136 Note Added: 0010618
2019-07-19 12:23 user136 Note Added: 0010619
2019-07-19 12:29 henry Note Added: 0010620
2019-07-19 13:03 user136 File Added: _bugreport_compressible.avi
2019-07-19 13:03 user136 File Added:
2019-07-19 13:03 user136 Note Added: 0010621
2019-07-19 16:09 michele Note Added: 0010624
2019-07-21 13:41 wwzhao Note Added: 0010628
2019-07-22 07:59 user136 Note Added: 0010630
2019-07-22 08:35 henry Note Added: 0010631
2019-07-22 11:52 user136 Note Added: 0010632
2019-07-22 12:18 henry Note Added: 0010634
2019-07-22 12:19 henry Note Edited: 0010634 View Revisions
2019-07-22 12:38 user136 Note Added: 0010635
2019-07-22 12:46 henry Note Added: 0010636
2019-07-22 14:58 user136 Note Added: 0010638
2019-07-22 15:31 henry Note Added: 0010639
2019-09-11 12:14 michele File Added: kOmegaAlpha.tgz
2019-09-11 12:14 michele Note Added: 0010733
2019-09-18 12:48 user136 Note Added: 0010762