View Issue Details
ID  Project  Category  View Status  Date Submitted  Last Update 

0003310  OpenFOAM  Bug  public  20190719 11:43  20190723 14:22 
Reporter  user136  Assigned To  
Priority  normal  Severity  major  Reproducibility  always 
Status  new  Resolution  open  
Platform  HPCCluster  OS  CentOS  OS Version  unknown 
Product Version  
Fixed in Version  
Summary  0003310: Variable density not correctly treated in turbulence models?  
Description  Hi! We're doing large scale simulations with interFoam (waterair, 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, kOmegaSST) 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. Explanation: 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 kOmegaSSTmodel 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 Reproduce  Extract and run the attached test case in OpenFoam7: unzip bugreport.zip cd _bugreport blockMesh; setFields; interFoam  
Tags  No tags attached.  

bugreport.zip (456,395 bytes) 

Currently in interFoam the turbulence model is formulated based on the divergence free volumetric flux for incompressible flow. There is no massflux 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 settingup interFoam with the compressible form turbulence models? 

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. 

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 ("libcompressibleTurbulenceModels.so");" kinto controlDict. Is that enough? If yes: It doesn't change the behaviour. 

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

> About compressible turbulence models: I'm not sure how to do that. I just now added "libs ("libcompressibleTurbulenceModels.so");" 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. 

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 .avifile) . 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) bugreport_compressible.zip (1,083,703 bytes) 

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 kw (or, if you prefer, kwSST) model, as described here: https://www.sharcnet.ca/Software/Ansys/16.2.3/enus/help/flu_th/flu_th_komega_turb_damp.html 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 kw with the turbulence damping. Just my two cents. 

Hi, the current implementation of komega SST model will generate overpredicted turbulent viscosity (eddy viscosity) near the interface of twophase flow. One workaround is introducing buoyancy effect into komega SST model, which requires solving for 'd(rho*k)/dt' etc and add a buoyancy generation term for kequation, 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 buoyancymodified kω and kω SST turbulence model for simulating wave breaking under regular waves using OpenFOAM®. Coastal Engineering, 138, 49–65. https://doi.org/10.1016/j.coastaleng.2018.04.011 [2] Larsen, B. E., & Fuhrman, D. R. (2018). On the overproduction of turbulence beneath surface waves in Reynoldsaveraged Navier–Stokes models. Journal of Fluid Mechanics, 853, 419–460. https://doi.org/10.1017/jfm.2018.577 

I completely agree with michele and wwzhao, that the results of kOmega 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 (https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19930013620.pdf). 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 airwatersystem, the density might be detrimental. But again: I think the model should reproduce what was written by Menter. 

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 twophase 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 farfield 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. 

I think we must take care not to mix up several aspects.  Buoyancy corrections for kOmegaSST 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 kOmegaSST is just not able to reproduce that.  Apart from those considerations I think, that the implementation should reflect the original paper by Menter 

 Buoyancy corrections for kOmegaSST 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 kOmegaSST is just not able to reproduce that. It is likely that a twophase 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 twophase 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 nonconservative 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. 

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 densityaware variant of kOmegaSST as a starter for further enhancements (damping etc.)? 

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 wallfunction for the interface. 

I tried the implementaion by Devolder (https://www.researchgate.net/publication/324798236_Performance_of_a_buoyancymodified_ko_and_ko_SST_turbulence_model_for_simulating_wave_breaking_under_regular_waves_using_OpenFOAMR). 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: https://github.com/BrechtDevolder/buoyancyModifiedTurbulenceModels OF1906: https://gitlab.com/brecht.devolder/buoyancyModifiedTurbulenceModels 

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 nonconservative. It would also make sense that the buoyancy term were implemented as an fvOption rather than hardcoded 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. 
Date Modified  Username  Field  Change 

20190719 11:43  user136  New Issue  
20190719 11:43  user136  File Added: bugreport.zip  
20190719 11:43  user136  File Added: results.0000.jpg  
20190719 11:43  user136  File Added: results.0002.jpg  
20190719 11:43  user136  File Added: results.0006.jpg  
20190719 11:54  henry  Note Added: 0010616  
20190719 12:01  henry  Note Edited: 0010616  View Revisions 
20190719 12:16  henry  Note Added: 0010617  
20190719 12:22  user136  Note Added: 0010618  
20190719 12:23  user136  Note Added: 0010619  
20190719 12:29  henry  Note Added: 0010620  
20190719 13:03  user136  File Added: _bugreport_compressible.avi  
20190719 13:03  user136  File Added: bugreport_compressible.zip  
20190719 13:03  user136  Note Added: 0010621  
20190719 16:09  michele  Note Added: 0010624  
20190721 13:41  wwzhao  Note Added: 0010628  
20190722 07:59  user136  Note Added: 0010630  
20190722 08:35  henry  Note Added: 0010631  
20190722 11:52  user136  Note Added: 0010632  
20190722 12:18  henry  Note Added: 0010634  
20190722 12:19  henry  Note Edited: 0010634  View Revisions 
20190722 12:38  user136  Note Added: 0010635  
20190722 12:46  henry  Note Added: 0010636  
20190722 14:58  user136  Note Added: 0010638  
20190722 15:31  henry  Note Added: 0010639 