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

0003344  OpenFOAM  Contribution  public  20190909 12:23  20191001 11:24 
Reporter  Wenyuan Fan  Assigned To  
Priority  none  Severity  feature  Reproducibility  always 
Status  new  Resolution  open  
Product Version  
Fixed in Version  
Summary  0003344: Issues with turbulence modeling in isothermal VOF solvers  
Description  Under the VOF framework, the flow of the isothermal mixture belongs to the variabledensity incompressible flow category. For such flows, VOFbased solvers of OpenFOAM fail to construct the correct governing equations for turbulence modeling. varRhoTurbVOF contains a set of newly designed VOFbased solvers which could use the desired governing equations for turbulence quantities.  
Steps To Reproduce  Details could be found at https://doi.org/10.1016/j.cpc.2019.106876 and https://arxiv.org/abs/1811.12580v2 Implementations for different OpenFOAM versions are provided at https://github.com/wenyuanfan/varRhoTurbVOF  
Tags  No tags attached.  

The flow is treated as incompressible in interFoam and hence the divergence of the mixture velocity is zero so it is valid to use the incompressible form of turbulence models. There are issues with the meaning of shear and turbulence generation in the interface region due to the unphysical nature of interface capturing but using a compressible formulation for the turbulence models does not resolve this. A more general formulation is provided in compressibleInterFoam in which mixture and separate phase turbulence models are supported, however no special handling for the interface capturing region is provided. 

Thank you for the quick reply! I totally agree that the mixture flow is divergencefree, but I do not think it is appropriate to use the "incompressible form of turbulence models" in OpenFOAM. The reason is that such incompressible turbulence models are designed for strict incompressible flows, where the density field is constant (unity in OpenFOAM). However, the density of the mixture is changing with the phase distribution, meaning that the flow is incompressible but with varying density. Therefore, those strict incompressible turbulence models are no longer valid for isothermal VOF solvers. I also agree with you that there are issues with the meaning of the shear and turbulence generation in the interface region, which are caused by the interface capturing method itself. We do not think our solution could solve all these issues. But, I do believe our solution is advantageous in at least three aspects. First, it is mathematically consistent. Since the flow is variabledensity incompressible, we need to consider the varyingdensity effect in turbulence modeling, which cannot be achieved by using the "incompressible form of turbulence models" in OpenFOAM. Second, it is much less sensitive to mesh refinement, as we described in the paper. Third, it produces better results for the case we discussed in the paper. 

The incompressible turbulence models are valid for each phase away from the interface because each phase is incompressible. In reality the density is a constant in each phase and it is not mathematically consistent to treat the system as a variable density mixture; this is a poor approximation. None of the turbulence models are valid in the interface capturing region and formally special treatment is required there. The current approach in interFoam is particularly robust and has been used successfully by thousands of users for thousands of cases since I created it in the late 90's. However, if you prefer the approach in which a variable density mixture approach is used for turbulence you can run compressibleInterFoam which support both this form and separate turbulence models for the two phases. It can be run compressible or incompressible in either phase by choosing the appropriate equation of state. 

To be more accurate and precise, I will use self consistency, instead of mathematical consistency. For a VOF solver, we need to accept all advantages and flaws of VOF, e.g., the favorable mass conservation feature and the questionable interfacial region. Therefore, the twophase system should be treated as a varyingdensity mixture as long as the VOF method is used. Concerning the self consistency that I have just mentioned, it simply means we should use the same treatment/assumption in the same solver. For instance, the variabledensity effect is considered in the momentum equation, then it should be considered in the turbulence modeling as well. For the same reason, if strict incompressible turbulence models were used, they should be used together with a strict incompressible momentum equation, where the varyingdensity effect is not considered. However, neither treatment is used in interFoam. I agree with you that we need special treatments for turbulence modeling when it comes to the interfacial region. And this is still a challenging topic. You are right that one could use compressibleInterFoam as an alternative. We tried this last year and it worked. However, for other VOF solvers that do not have the corresponding compressible versions, this approach will not work. Anyway, thank you for handling this issue! I really appreciate your enormous contribution to the community. 

Adding another VoF solver which duplicates some of the functionality of compressibleInterFoam introduces substantial maintenance overhead. It would be better if functionality is merged rather than duplicated. Any change to interFoam must preserve the current effective, robust and well tested functionality. It would be possible to add the other turbulence options currently available in compressibleInterFoam to interFoam, i.e. variable density mixture formulation and separate phase turbulence models and this would be easier to maintain than your proposal of having a completely separate solver structure. However given that this functionality is already available in compressibleInterFoam it is not clear how much interest there would be in this work. An alternative would be to make parts of compressibleInterFoam selectable, e.g. the energy equations, and this could then replace interFoam entirely. 

> For instance, the variabledensity effect is considered in the momentum equation, then it should be considered in the turbulence modeling as well. I don't see why this is necessary given that the interface region in erroneous anyway. A better approach would be to add modeling for the interface region which might be based on the variable density or current incompressible formulation or even separate turbulence models for the two phases. Once we have this model we can decide which turbulence formulation should be used. 

It is true that there will be repetitions if we want to keep both the strict incompressible and the variabledensity incompressible solvers. Because the only difference between these two types of solvers is the turbulence modeling part, which only has several lines of codes. Therefore, one possible solution is making the turbulence modeling part selectable. Of course, we need to find a smart way to do it. Talking about compressibleInterFoam, I think it is another reason why we need to make such changes. Some solvers have two versions, say, incompressible and compressible. To my understanding, this classification enhances the solver performance, but it should not violate the consistency between solvers. For instance, we can use both simpleFoam and rhoSimpleFoam (by choosing the appropriate equation of state) to simulate a flow with constant density. There might be slight difference in the results due to roundoff errors. However, these two results should always be numerically equivalent, and the only difference should be that simpleFoam runs faster than rhoSimpleFoam. When it comes to interFoam and compressibleInterFoam, this is not the case. Also, interFoam is not the only solver affected by this issue. For instance, interMixingFoam also has this issue, but we cannot use "compressibleInterMixingFoam" to avoid this problem, because it does not exist. Regarding the necessity of using consistent formulations for both the momentum equation and the turbulence modeling part, I think it is a requirement for a model that all submodels should not contradict with each other. Otherwise, the model does not have a solid theoretical basis. I totally agree with you that we need to develop better turbulence models for twophase flows. However, the current problem is that we only have two solutions, which both look erroneous. One solution is erroneous in a more selfconsistent manner than the other in the sense that all the treatments are in line with the basic assumption of VOF. If I have to make a choice between these two, I would select the more selfconsistent one. Anyway, I think it would be beneficial for the community to understand this issue. No matter which solver they choose, it would always be nice if they know the rationale behind it. 

> but we cannot use "compressibleInterMixingFoam" to avoid this problem, because it does not exist. It would be trivial to create a compressibleInterMixingFoam solver if required, but so far no one has requested or contributed it. However, I am no longer convinced that interMixingFoam or compressibleInterMixingFoam are good approaches and it would be better to solve for mixing in terms of species instead of additional phases which mix. As it happens this would be easier in the compressibleInterFoam framework rather than interFoam because the standard thermo packages used in compressibleInterFoam have composition options. Given that interFoam has been used successfully for 20years with the current turbulence modeling framework any change to this would need to be validated on a VERY wide range of cases, particularly if a replacement of the current approach is to be considered. How many cases and for what range of problems have you tested the alternative formulations on? 

I understand the popularity of interFoam, so we need to be very careful with any changes. I also admit that we have only tested our solvers with few cases. Therefore, much more tests are needed. However, this is far beyond our ability, and can be only achieved by efforts from the community. By the way, I have just noticed Issue 0003310, which means that there is certain demand from the community. I hope our solvers could fulfill such demands. 

OK, we can wait until the community tests your changes or compares interFoam with compressibleInterFoam on a range of cases to confirm or otherwise the need for them. Can you provide a version of interFoam in which the three options for turbulence modeling are supported with a runtime selection mechanism, basically an extension of the mechanism I implemented in compressibleInterFoam? 

Yes, I will try to create a version of interFoam which supports runtime selection for turbulence modeling. I hope I can finish it in the coming days. 

You may be also interested in testing an approach consisting on turbulence damping at the interface. Here a link to a kOmega model derivation, tested in OF5 https://bugs.openfoam.org/view.php?id=3310#c10733 

Thank you for the information! Actually, we have been working with Egorov's turbulence damping model for some time, and have made some improvements to the model, which could be found at https://doi.org/10.3390/fluids4030136 In the model you provided, "dn" is given as a user specified scalar, instead of a volScalarField which should be calculated based on the local mesh and flow information. Also, we proposed an asymmetric damping treatment in the paper. The issue I would like to address in this thread is about self consistency. Of course, a selfconsistent model does not always guarantee correctness due to certain limitations of the model. As described in the above mentioned paper, varRhoInterFoam cannot give reasonable predictions without using Egorov's turbulence damping model. Then one may argue that if a selfconsistent model does not work well, why should we use it? The answer is simple. We can at least make sure that the problem is not caused by the selfconsistency of the model. 

Hello Wenyuan Fan. Let me understand your point of view: interFoam considers two immiscible incompressible fluids with different densities (air, water). Each phase is incompressible (and the incompressible equations apply, both for momentum and turbulence). At the interface we have a physical discontinuity. Physically we know that there exist a relaminarization of the flow at the interface. There could be also other phenomena (like interphase mass/momentum/turbulence transfer). In the real world, at the interface, the density gradient goes to infinity. A VOF approach "smears" this discontinuity trying to keep the interface compressed on, let me say 2 to 5 cells. If you consider a compressible version of the turbulence model, the only terms that differ from the incompressible version are those related to density gradients. So you are trying to compute these gradients on a region of 2 to 5 cells where large density gradients occour. Apart from the physical inchoerency of calling "density gradient" a region where, in reality, there is a phase change (the DENSITY GRADIENT is only a result of the VOF ability to compress the interface, so it HAS NO PHYSICAL MEANING). Also, due to the fact that large gradients occour on a small number of cells, have you quantified the discretization errors of such an approach? In my opinion there could be large errors. I think it may be preferrable a different approach, i.e. simply model the interface interaction effects. So I use the Egorov approach for the relaminarization terms. The characteristic lenght "dn", in strongly anisotropic cells, is simply the cell height at the interface (and not the cube root of the cell volume, as implemented in certain commercial codes). In my procedure, it should not be "automatically" computed from the mesh and flow characteristics (provided the computational mesh adheres to strict requirements in quality and geometry at the interface). At the end of the day, all is multiplied by a B factor that is a "large" number. In real world problems with complex geometries this approach proves to be accurate (and consistent with experimental validation). In my experience this approach, combined with a good mesh, adhering to quality requirements (especially at the interface), is the main ingredient for obtaining useful results. 

Hi Michele, Regarding turbulence modeling in VOF, my key point is self consistency, as I have already stressed in previous posts. Let's consider the issue in another way. I totally agree with your argument about the DENSITY GRADIENT that it's unphysical and there might be large discretization errors. Then why should we include this term in the momentum equation? If a strict incompressible momentum equation were used, I would agree the idea of using strict incompressible turbulence models. In my opinion, the statement "the DENSITY GRADIENT HAS NO PHYSICAL MEANING" is true in the real world, but is dangerous in the realm of VOF. Concerning the Egorov approach, it is useful. Therefore, we are using it and trying to make it better. I agree with you that the cubic root of cell volume is a bad treatment. "dn" could be a constant for given meshes and flow conditions. However, for general cases, where the interface is moving, "dn" should not be treated as a constant since the interface orientation is changing. Discussions about the calculation of "dn" could be found in our paper. Regarding the value of "B", we need to find an appropriate number, not a random "large" number, or a "large enough" number. This is also discussed in the paper. Anyway, the Egorov approach is useful, but it needs fine tuning of the input parameter. The method might be problematic in the absence of experimental values because we don't know how large should "B" be. 

Ok thanks for the clarification. My conclusion is the following. Regarding the declared "self consistency" of using a compressible turbulence model in interFoam. Are you aware that you are implicitly assuming that a general compressible turbulence model, that has been developed for substantially different types of flow, with coefficients calibrated on typical compressible flow benchmarks, has no specific modelling/calibration for the VOF interface? A turbulence model is never an exact closure: it has a series of coefficients that are tuned on specific scenarios, and is capable to capture only the phenomena that are included in its terms. At least a recalibration of the coefficients on known benchmarks should be performed before using blindly a turbulence model on a completely different scenario. Moreover, are you sure that the nonuniform density diffusion term is a correct model for the relaminarization mechanism? This term (in the k equation) is just a consequence of density variation in compressible flow, but does not "inform" the turbulence model that an interface do exist. In my opinion there must exist a specific correction in order to inform the turbulence model that there is an interface, the same way that in more advanced compressible turbulence models there are terms that account for the presence of shock waves interfaces (typically in the dissipation equation). Regarding the shock waves interface terms, consider that we have typical density ratio of about 5 in (hypersonic) strong shocks, whereas in a interFoam simulation with airwater the density ratio is about 800, so these specific terms are really needed. And there we return to the Egorov damping term: it "informs" the turbulence model of the presence of an interface and provides a way to deal with it. Regarding the coefficients, B should be "large enough" to suppress the turbulence: provided is large enough, there is a large range of variation (exactly like for the coefficient of the solidwall function). Before stating that a model should be replaced in complex industrial applications, a lot of more work should be done. I would be more than happy that, after 15 years of successful usage of the Egorov model, some advancement is made. Personally, if I had to construct a new turbulence model for VOF, I would go back to the relaminarization mechanism at the interface: i.e. the suppression of the interfacenormal fluctuations (v2). This mechanism is strongly anisotropic and consequently the Boussinesq hypothesis is debatable. So I would try to derive a model like the elliptic relaxation model (Durbinâ€™s v2f), with specific interfaces sources/sinks on the equations, linked to a nonlinear (or algebraic) Reynolds stress tensor. If such a model proves to be robust, numerically efficient, and accurate (on a sufficient range of benchmarks, with different meshing strategies, and with a demonstration of its mesh convergence capabilities), then I would consider to replace the simple Egorov dissipation source. 

I am quite aware that we are using singlephase turbulence models to model multiphase flows, which is definitely problematic. However, this is the approach that has been widely adopted in different software and studies. Your comments on calibration and recalibration are a very good point, because we even have such problems for singlephase flows. However, I don't think such comments are fair to the variabledensity incompressible turbulence models, because the strict incompressible turbulence models have the same issues as well. Therefore, I don't see the link between the calibration/recalibration issue and the selfconsistent formulation. In addition, the nonuniform density diffusion term purely comes from the mathematics derivation, we have never claimed that it is a model for the relaminarization mechanism. Also, I would like to emphasize that a flow with variable density is unnecessarily compressible. I agree with you that the desired VOF turbulence models should be able to detect the location of the interface. However, I still believe that there is no large enough value for B. As we discussed in the paper, the nonzero behavior has been observed for k near the interface both in experiments and DNS, meaning that there is certain level of turbulence around the interface. A large enough B simply kills all the turbulence. Thank you for sharing your thoughts about the turbulence modeling for VOF. I agree with you that the anisotropic Reynolds stress tensor should be an important ingredient. However, adding such terms is not an easy work, so many people may use LES directly. Of course, there should be additional treatments in LES as well. 

Hello Henry, I have created a new version of varRhoInterFoam which supports runtime selection of turbulence modeling. Please find it in the attachment. The separatephase turbulence modeling is not included since we are not working with it. For the moment, the new turbulence class is written for immiscibleIncompressibleTwoPhaseMixture transport model. However, it can be extended to other transport models (and VOF solvers) by using templates. OpenFOAM7selectable.tar.gz (17,700 bytes) 

> The separatephase turbulence modeling is not included since we are not working with it. It would be much better if all three options were handled in the same system based on the approach in compressibleInterFoam. The merge into OpenFOAMdev can wait until this is complete. Following feedback from the community it appears that interface turbulence damping is more important than changing the turbulence model formulation. Currently we do not provide interface damping for the turbulence but it would seem quite straightforward and could be done in an fvOption. Can you provide an fvOption for this? 

Dear colleagues, here is an update to this and the issue reported here: https://bugs.openfoam.org/view.php?id=3310 We checked the impact of Wenyuan's implementation of incorporating rho consistently into the turbulence models for interFoam. Here I focus _only_ on incorporating the density in the equation, not on damping etc. as implemented by Michele or Brecht Devolder. Damping clearly enhances the physics if done correctly, but I think first the numerical part should be "best". Our investigations show that incorporating rho in the equations clearly enhances the modeling of turbulence in the "heavy" phase of interFoamcomputations. The attached pictures show a testcase. It is the flow through a weir complex (flow is left to right). Picture 1 shows the geometry, 2 shows velocity right below the water surface, 3 shows velocity on a vertical slice through the domain in flow direction. This slice is regarded in detail in the following pictures. Picture 4 shows the development of k computed with "normal" interFoam. The white line is the water surface. The upper image part shows the slice directly behind the weir. K should increase in the water due to local dissipation. But instead a lot of k is generated in the air and this effortlessly diffuses into the water. This is due to the constant density, which gives water and air the same "weight" in the PDE. Further downstream (lower image) we would expect k to be highest at the bottom (friction). Again this is spoiled with k coming from the air. Picture 5 shows the same with the modified turbulence model by Brecht Devolder (https://bugs.openfoam.org/view.php?id=3310, rho in turbulence model), but without the damping he implemented . The results are much better in the water: Near the weir k is dominated by local effects in the water. In the farfield bottom friction is visible. The results in the air are not very nice (numerically), but they do not influence the water too much. Picture 6 shows the results with varRhoInterFoam by Wenyuan. Again the results in the water are better than in the original version, but here the air looks numerically better, too. Comment: The simulations were still running when I made the pictures, so they do not all show the same time. In my oppinion incorporating the density in the equation enhances the quality for people that are interested in simulating the water phase. But for people interested in the air it might (!) be worse. Actually I think the usergroup of "water" guys is much bigger, so I vote for using the correct density. Maybe as a flag? Based on that, additionally a damping near the water surface (fvOptions) would help. Best, Carsten 2_velocity.jpg (392,201 bytes) 3_slice_velo.jpg (319,880 bytes) 

Hello Henry and Carsten. Thank you for the comments and feedback! Regarding the separatephase turbulence modeling, I think it might violate the selfconsistency that we have been talking about. Since such an approach is designed for the twofluid method where each phase has its independent momentum equation, it will be confusing to use this approach in VOF where we only have one momentum equation for the mixture. This is my main concern about adding this type of turbulence modeling. Regarding the fvOptions for turbulence damping, I will provide one in the coming days. Actually, I have been seeking a better way to add the damping term such that the users don't need to keep a copy of the fvOptions file or use a separate turbulence model. Any suggestions? By the way, I really like the following sentence by Carsten. Damping clearly enhances the physics if done correctly, but I think first the numerical part should be "best". 

Please find the fvOption for turbulence damping in the attachment. I have tested it with OpenFOAM7 and OpenFOAM1706. It should work with other recent versions as well. The usage could be found in "turbulenceDamping.H". 

turbulenceDamping.tar.gz (3,446 bytes) 

Thanks for the turbulence damping implementation, I had a quick look and there are some issues. Is it applicable to both komega and kepsilon models? Will it work for interFoam and compressibleInterFoam? Also the layout of the code is rather different to the conventions in OpenFOAM: https://openfoam.org/dev/codingstyleguide/ Could you update it before it is included in OpenFOAMdev? It would save a lot of time. Do you have results with this damping which shows the effect on the two approaches to mixture turbulence? 

Please find the updated source file in the attachment. The original model needs to work with omega equations, and my implementation only handles omegabased turbulence models. I think such a treatment is also adopted in ANSYS FLUENT. This implementation deals with mixture turbulence models, both strict incompressible and full forms. Therefore, it will always work for interFoam. However, it will only work for compressibleInterFoam provided that the separatephase turbulence modeling is not used. We tested this model with varRhoInterFoam, the results are given in the paper. Of course, one could test the combination of the damping and interFoam, since, as I mentioned above, the implementation could handle strict incompressible turbulence models as well. turbulenceDamping.tar2.gz (3,510 bytes) 

Is there a fundamental reason why this form of damping cannot be applied via an epsilon equation? Many many users use kepsilon forms with interFoam and would also benefit from turbulence damping. I am glad that the implementation is general enough to work with compressibleInterFoam with mixture turbulence and I understand that the model is only for mixture turbulence modeling. A completely different and probably more physical approach would be needed with separate turbulence models in the two phases. 

Regarding the fundamental reason, I think it is related to the different nearwall behaviors of omega and epsilon. Because this damping modification is based on the phenomenon that, to the gas, the liquid behaves like a solid wall. For omega, there is an asymptotic solution that we can use. And this is why the parameter "B" is included in the model. However, there is no such solution for epsilon. From this point of view, this model can only work with omegabased models. Of course, this is a way to introduce similar corrections to epsilonbased models, which is less physically sound. I will try to implement it and evaluate the performance. 

Given the number of users who run interFoam and compressibleInterFoam with models other than kOmega we will need an implementation which also works for epsilon. Also given that omega and epsilon are related and each can be obtained from the other given k and that it is possible to convert epsilon equations into omega equations and vice versa (and this was done in the creation of the kOmegaSST model) I see no reason why the turbulence damping omega source cannot be converted into the equivalent epsilon damping source. It is not clear why different near wall behavior between different models is relevant here; I get the same near wall solutions for k epsilon and k omega models using corresponding wall function choices. 

You are correct. This is the way I mentioned. As I said, I need to implement it and evaluate the performance. However, this task has relatively low priority on my to do list, since I think including density is the basic requirement. Introducing damping may improve results. However, it should work with a selfconsistent turbulence model. 

OK, no problem, this can wait until the turbulence modelling is complete and support all three options and the turbulence damping support epsilon. Currently there is not much information about the difference in behavior of the current interFoam and compressibleInterFoam turbulence modelling with the turbulence damping, the indication given to me in the past is that the damping it the critical component and the difference in turbulence model form less important. Either way the community can test and use either option in the current OpenFOAM releases and try the damping model implementation you provided and report their findings. 
Date Modified  Username  Field  Change 

20190909 12:23  Wenyuan Fan  New Issue  
20190909 12:38  henry  Note Added: 0010708  
20190909 12:39  henry  Priority  immediate => none 
20190909 12:39  henry  Severity  major => feature 
20190909 13:23  Wenyuan Fan  Note Added: 0010709  
20190909 13:38  henry  Note Added: 0010710  
20190909 16:17  Wenyuan Fan  Note Added: 0010711  
20190909 16:30  henry  Note Added: 0010712  
20190909 16:34  henry  Note Added: 0010713  
20190910 10:33  Wenyuan Fan  Note Added: 0010714  
20190910 11:08  henry  Note Added: 0010716  
20190910 16:38  Wenyuan Fan  Note Added: 0010721  
20190910 16:45  henry  Note Added: 0010722  
20190910 17:10  Wenyuan Fan  Note Added: 0010723  
20190911 12:18  michele  Note Added: 0010734  
20190911 12:55  Wenyuan Fan  Note Added: 0010735  
20190912 10:06  michele  Note Added: 0010736  
20190912 11:52  Wenyuan Fan  Note Added: 0010737  
20190912 19:48  michele  Note Added: 0010738  
20190912 22:23  Wenyuan Fan  Note Added: 0010739  
20190917 10:16  Wenyuan Fan  File Added: OpenFOAM7selectable.tar.gz  
20190917 10:16  Wenyuan Fan  Note Added: 0010743  
20190917 10:42  henry  Note Added: 0010745  
20190918 12:48  user136  File Added: 1_geometry.jpg  
20190918 12:48  user136  File Added: 2_velocity.jpg  
20190918 12:48  user136  File Added: 3_slice_velo.jpg  
20190918 12:48  user136  File Added: 4_slice_interFoam_k.jpg  
20190918 12:48  user136  File Added: 5_slice_interFoam_rhoTurb_k.jpg  
20190918 12:48  user136  File Added: 6_slice_varRhoInterFoam_k.jpg  
20190918 12:48  user136  Note Added: 0010761  
20190918 14:58  Wenyuan Fan  Note Added: 0010764  
20190919 13:39  Wenyuan Fan  Note Added: 0010765  
20190919 13:42  Wenyuan Fan  File Added: turbulenceDamping.tar.gz  
20190919 13:59  henry  Note Added: 0010766  
20190919 15:15  Wenyuan Fan  File Added: turbulenceDamping.tar2.gz  
20190919 15:15  Wenyuan Fan  Note Added: 0010767  
20190919 15:31  henry  Note Added: 0010768  
20190919 16:52  Wenyuan Fan  Note Added: 0010769  
20190930 14:45  henry  Note Added: 0010783  
20191001 10:03  Wenyuan Fan  Note Added: 0010784  
20191001 11:24  henry  Note Added: 0010785 