A discussion on turbulence models


I wanted to open a topic on turbulence models since I often see questions about it but it’s not yet covered in the SimScale documentation (it’s in progress :wink: ). Here is a list of the models currently supported by SimScale

  • Laminar
  • LES Smagorinsky
  • LES Spalart-Allmaras
  • RANS k-epsilon
  • RANS k-omega
  • RANS k-omega SST

A couple of questions to get this discussion started:

  • What are the best resources for understanding these models and how they should be used?

  • How can we be sure that the turbulence model that has been selected, is in fact within an acceptable range of accuracy? Does anyone have examples?

Y+ range recommendations in documentation
Laminar-turbulent transitions models?
CFD - Wing Flow Simulation
'Performance of a split air-conditioner in a domestic setting' simulation project by pankajkumar979
Community Digest July 2016

Dear all,

I just want to add some information related to the topic. For all who are interested in turbulence modeling in more details should understand what we are doing here and hence it is necessary to understand the general equations first. The main literature that I used in that topic of turbulence modeling are [Ferziger] [Wilcox] and [Bird]. Of course there are plenty of literatures. A summary of the literature that I am using can be found on http://holzmann-cfd.de/index.php/en/literature

In the link you also find the ISBN numbers and titles of the books (maybe you have them already in your library). If you do not want to spend to much time on all the books, you also can use my publication about » Mathematics, Numerics, Derivations and OpenFOAM « that is published on my webspace or on researchgate.

Now I want to give you a brief overview, without going to deep into detail (turbulence modeling is a taff topic)

If we use laminar, the equations are valid for each kind of flow (yes - even turbulent flows). The fact is simple, laminar will solve the general Navier-Stokes equations without any need of modeling. Hence we get natural vortexes and if we make the mesh fine enough, we can resolve the turbulence in all the beauty (DNS). But using laminar in turbulent flow fields will lead to extreme fine meshes to resolve the so called Kolmogorov vortexes. The fact that this is never possible for complex geometries, we only use the laminar mode for LAMINAR flow fields. If we have turbulent flow, we need to simplify things.

RANS Models
As already mentioned, resolving the turbulence means to resolve the smallest vortexes that lead to extreme fine meshes and therefore the time step will also decrease dramatically. Thus, we introduce the so called Reynolds time-averaging method that will lead to the Reynolds-Averaged-Navier-Stokes equations. The complete derivation of these equations can be found in my publications. Deriving the RANS equations lead to unknown terms that have to be expressed with already known quantities. Hence, we introduce some hypothesis to express the unknown with other quantities that finally lead to the already known turbulence models. In the turbulence models we have a lot of approximations and assumptions. If you check out the derivations of the equations you will get a feeling how many terms are approximated and then you get a feeling that turbulence modeling is always a big assumption.

I do not want to go too deep into detail, but there are differences between compressible and incompressible RANS equations.

  • For the incompressible equations we always get the mean Reynolds averaged quantities U, p, k, epsilon, T etc. (further more, the pressure p is not the real pressure, here we get an addition kinematic pressure)
  • For the compressible equations the derivation is even more complex, hence we introduce the Favre-weighted Reynolds time-averaging method. The results are the mean Favre averaged quantities U, p, k, epsilon, T etc. Note that the Favre-Averaging concept is correct only in a mathematical point of view, but not in the physical point of view.

This models are mostly based on the assumption of isotropic turbulence, that is not always true, especially if we are talking about turbulence near walls or after some rigid body etc. .

LES Models
LES models are similar to RANS models but here, instead of averaging the whole flow, we resolve the big vortex scales of the turbulence and model the smallest turbulence scales. This is a much better approach because if we go to smaller and smaller scales, the isotropic assumption is valid in much cases. Hence, we get a better physical approach for the turbulence.

What happen to the equations if we use turbulence models
Luckily nothing special due to the fact that we make the Reynolds and Favre averaging in a way that we end up with the (more or less) same equations as we have it for the general flow field. The only important thing is, that we assume a higher viscosity (based on the turbulence), that lead to higher diffusion of momentum, concentration, temperature, enthalpy etc. . The molecular viscosity is increased by this special amount of turbulent viscosity that is also known as eddy viscosity. The turbulence models are simply calculating the eddy viscosity.

Which model should be used
That depends on the problem you are focusing on. You will realize that if you use two different turbulence models with the same settings, you will end up with two different solution. As I already mentioned, most models are based on the isotropic turbulence (energy cascade). If you are interested in regions near the surface you should choose a model that has some corrections for the near wall field or even uses an anisotropic hypothesis. If you know that your turbulence is really directed only in two directions (if we assume a 3D flow), you should also use a turbulence model that is not based on the isotropic assumption (there was a nice demonstration video but I could not find it now).

Of course, the more complex the models get, the more effort has to be done to calculate the stuff.
In a very general way we can say for the three simplest models:

  • k-epsilon for far fields
  • k-omega for near fields
  • k-Omega-SST is the combination of both, it switches between them using the dimensionless y+ distance to the wall

I did not check the source code of FOAM for the last arguments but it should be somehow like that. Why is it like that? Its simple, the dissipation of the turbulent kinetic energy field epsilon can be calculated very well in the far field whereas close to the walls it will calculate to high turbulence. Here it is better to use the frequency omega.

Hopefully this can help you to get the first overview of turbulence modeling.


Very interesting and useful homepage @TobiasHolzmann!

Since I want to learn more about CFD and turbulence in general your post comes along just at the right time. Hope to benifit a lot from your knowledge, experience as well as your videos :slight_smile:




@TobiasHolzmann great intro!


since a lot of people are going to be doing a lot of simulations involving walls (eg flow over a car body, flow through pipes etc), it is essential to understand the relevant near wall meshing criteria depending on the turbulence model used.

First and foremost a prism (boundary) layer mesh is essential to get realistic results in all wall-bounded cases (ordinary cartesian meshes are useful only in the simplest of cases). Then, if y is the wall-normal coordinate, it is paramount to determine what value of y1+ ie the distance (non-dimensionalised by the viscous length scale) of the first wall-normal grid cell centre will be. In most cases, the user determines this choice and usually, there are 2 options; either the first cell centre lies in the viscous sublayer of the boundary layer (which in general means y1+ < 5) or it lies in the log-layer of the boundary layer (which means in general y1+ > 30). The choice of y1+ and turbulence model must be made in unison. In general, a y1+ > 30 would be employed with k-epsilon models (with wall-functions calculating the flow for y1+ < 30) whereas y1+ < 5 would be employed with k-omega models. These ideas add to those mentioned @TobiasHolzmann in ‘Which model should be used’.

In case of advanced techniques such LES or DNS, the criteria on near-wall meshing become increasingly strict. In most of these cases, one will see y1+ < 5 and in many cases even < 1. The near-wall meshing criteria are extremely rigorous here as considering the physics, it is important to resolve all the high gradients that are present in the flow close to the wall.



This is a great topic, @AnnaFless, and one that will probably be in discussion forever. Commercial code CFD help documentation is a great starting point, as it is more application based than theory based.

Also, NAFEMS is a wonderful unbiased resource. They have a list of links here: https://www.nafems.org/about/technical-working-groups/cfd/links/ and a list of publications here: https://www.nafems.org/about/technical-working-groups/cfd/publications/

How to ensure the turbulence model is within an acceptable range of accuracy? Start by converging Turbulent Energies/Minimizing Residuals in the computation, and finish by matching your flow results to Test Data!

Is there a public area designated for “Benchmark Projects”? This would be the most direct way to share this knowledge with Simscale users, IMHO.


@mananthakkar87 - thanks for adding onto the topic, great information!

Thanks @fastwayjim, turbulence models seem to be a challenging topic so it’s good to have a place to talk about it! We currently have validation cases in our documentation. However, this is a good idea to create a space in the forum to share this information as well. Maybe you could create a topic to get this started in #using-simscale?


@TobiasHolzmann - I think we need a “Tweet this” button here in the forum. Fantastic post!


@dheiny Thanks for the replay and the feedback. Maybe I will extend the post by the y+ as already mentioned above but then we had one post. Even I want to add some more literature for people who are interested int the topic. »Tweet it«


I think that the most hard to grasp aspect of turbulence modeling in SimScale is the use of automatic wall functions. I would like to focus on two turbulence models: the RANS k-omega SST and the DES of Spalart-Allmaras. In both it seems that automatic wall functions have been extensively used to promote an almost monotonic error decrease as the meshes are refined near the walls, bringing the adjacent to wall cell centers from y+ > 30 (inertial sublayer where the standard wall function would work) to the y+ < 1 range (laminar sublayer, where these turbulence models could be used without any wall function). And here came my questions (not very good ones, I realized later, see my next post below).

For the Spalart-Allmaras (under the LES label) SimScale uses the Spalding wall function ( http://cpp.openfoam.org/v3/a01670.html#details ) for wall shear stress calculation (when Wall function BC option is selected)? This wall shear stress is used to apply no-slip boundary condition (BC) on the momentum equation while nuTilda = 0 is taken as BC (at walls) for the nuTilda transport equation, right?

For the k-omega SST we need BC’s for momentum, k and omega. SimScale uses the Spalding wall function for obtaining wall shear stress (used in momentum BC) and k at the center of the cell adjacent to the wall (through the third equation in https://www.cfd-online.com/Wiki/Near-wall_treatment_for_k-omega_models ) while the omega value is taken from square root of the sum of squares of the values calculated (as also shown in this same link). Right or wrong?

Y+ range recommendations in documentation

Hi Luciano

No-slip is enforced by applying velocity of the wall on wall patches, because flow velocity immediately next to the wall has the velocity of the wall. The Spalding wall function is continuous and therefore allows nut of the wall patches to be consistent with non-dimensional wall units. I don’t know what nuTilda value is applied on wall patches on SimScale - only the staff members would know. But you are right about nuTilda = 0 as wall BC.

Again, I don’t have access to the SimScale code, so not sure what is implemented.


I have just learned that (at least when using one of the OpenFOAM solvers available inside SimScale) we do not need access to the SimScale source code in order to find out what are exactly the boundary conditions used when the “Wall function” option is selected. This information is available within the results compressed file that we can download after each case run (it is at the “0” subdirectory in that compressed file, we can look at the “nuTilda” file to know the BC for Spalart-Allmaras and at the “k” and “omega” files for k-omega SST BC; “alphat” file is also of interest if the problem involves heat transfer). In these files there are keywords and data entries read in OpenFOAM and it is in OpenFOAM source code and documentation that we finally find very detailed and precise information (although somewhat dispersed) about those boundary conditions.

For Spalart-Allmaras it can be seen that the automatic wall function in use in SimScale is based on an if statement that activates the Launder-Spalding* wall function when y+ is greater than a limit value (that can change with the wall roughness parameter E, but for smooth walls is 11.5301). At cells (close to walls) where y+ is below this limit value, no wall function is used (nu_t = 0). * Meaning that y+ = C_mu^0.25k^0.5y/nu_w definition is used. I have made a confusion with the continuous Spalding wall function (that is not used by SimScale) in my previous post.

In k-omega SST it seems that zero gradient BC is used for k. At the cells close to walls, omega is set equal to the square root of the squares of omega_Log and omega_vis, pretty much as explained in https://www.cfd-online.com/Wiki/Near-wall_treatment_for_k-omega_models . The nu_t calculation, however, is performed by the same method used for the Spallart-Allmaras model (in OpenFOAM I did not found any use of the fourth root of fourth powers of friction velocities that appears in the link above). Finally, it should be noticed that the production source term in the k transport equation is modified, when the wall function is used, to deal with the increased shear rate in the boundary layer (recognizing that the shear taking place in the laminar sublayer do not contribute to turbulence production).


Thanks for clarifying those BC’s.

Edit: I don’t think it is possible to tell whether ‘the fourth root of fourth powers of shear velocities’ is used in the computation of wall shear, as wall shear is not used in the turbulence model, ie., calculating wall shear is more likely a utility. The BC’s you have seen apply on the nut boundary field, whereas the blending of friction velocity relies on known yStar values. They are two separate things as far as I understand.


In order to complete my previous post about what I have learned reading the OpenFOAM source code (and did not found clearly stated anywhere else), it is worthy mentioning that OpenFOAM implementation of the Jayatilleke wall function for turbulent thermal boundary layers also uses an if statement, much in the same fashion used in Fluent according to equation (10.8.5) in http://jullio.pe.kr/fluent6.1/help/html/ug/node451.htm (with the difference in notation of yStar, identified as yPlus in OpenFOAM).


Hi Luciano

I don’t think you have understood what I said. You can’t tell if OpenFOAM uses ‘the fourth root of fourth powers of friction velocities’ to compute wall shear by merely looking at BC’s, because the BC and the wall shear calculation are two separate things.

In terms of yPlus and yStar. I believe OpenFOAM has got both in its wall functions. Depending on whether non-dimensional wall distance is computed by wall-normal velocity gradient or turbulent kinetic energy, the wall function uses the appropriate formulas.


Yes, Dylan, I have been experiencing some difficulty in learning (again or by the first time) what I have not looked at as much as I would like to in these 20 years past since I took a formal course in Turbulence Modeling - where we could not simulate very practical situations using an open source Fortran77 2D code (similar to a famous one used by Patankar in teaching, restricted to cartesian grids) and an under-stimulated access to Fluent. Not to speak about my poor communication in English skills…Thanks for worrying about my difficulties and helping me!

I think that, In some contexts, the phrase “shear stress calculation” may be not just referring to that one used for obtaining a friction drag coefficient, some kind of head loss, or the hydrodinamical load over a surface (at the post-processing stage of a CFD simulation). It may refers also to that of shear stresses implicitly involved in the diffusive terms of the RANS equations and in the production “source” terms appearing in the k and epsilon/omega transport equations (in an internal product by the deformation rate tensor that amounts to power/volume). Since we are talking about subjects scarcely found in a few textbooks I want to praise a recent and comprehensive one, by Moukalled et alii, that I just knew http://www.springer.com/gp/book/9783319168739

About that fourth root of fourth powers of friction velocities, I would like to say that I have nothing against it (is a recommendation of the greats Menter and Esch), but I feel that OpenFOAM developers did not deemed its implementation advantageous.



I have thought what the title of the page https://www.cfd-online.com/Wiki/Near-wall_treatment_for_k-omega_models said exactly - wall treatments of omega-based models. Under no condition, would I know that you also referred to the source terms, and I have explicitly used the term ‘wall shear’. But…thanks for clearing that out.

By the way, I wouldn’t call it implicit if it is a source. Do you mean implicit source as in the source term is solved iteratively?

I have stopped thinking ( and see no point ) why something is not implemented in OF. I simply implement what I believe is right. It is good, though, that you have clarified the BC’s on SimScale.

I have to say I am clearly less experienced ( and have spent less time ) in CFD than you. The only kind of stuff I do is write solvers in OF.

Coming back to the important stuff:

  1. I wouldn’t solve SA the way it is done on SimScale from what you have described - I wouldn’t use k to predict y+ in this case.
  2. I would use either the Spalding law or the improved wall function for nut by Allmaras et al (2012)
  3. I am hardly a SST user, maybe you can point out the advantage of ‘fourth root of fourth powers of friction velocities’.


I have to say that presently my intention is basically to use SimScale to give a taste of CFD to the undergraduate students to whom I give Transport Phenomena classes. Thanks to SimScale I was able to simulate experiments that these students perform in laboratory and I am satisfied by the agreement of these simulations with the measurements (of wall pressure, mean velocity using Pitot tube) and calculations (of Nusselt number) made by them. The question about the automatic wall function was important to give the students correct information about what happens as the grid is refined, especially near walls (and why people worries less about the y+ now than in the past). Precise information will be even more important for reporting this didactic experience to my peers. So I am not doing any ambitious research in turbulence modeling.

I believe that the reason for using the “fourth root of fourth powers of friction velocities” blending formula is essentially the same for using the Spalding continuous wall function or the analytic solution of the SA equation: to obtain smoothly varying (and agreeing with DNS and experimental results) quantities as we refine the grid and have cell centers in the range 1< y+ < 30. My guess is that numerical tests revealed that the simple switching (at the intersection found near y+ = 11.6) between the wall function and the viscous formula for shear stress, is already sufficiently smooth and accurate (agrees satisfactorily with DNS and experiments) when used with k-omega SST or SA turbulence models. Only the omega calculation (in k-omega SST) really demanded the use of a special blending formula. What should be deemed as sufficiently accurate is a rather subjective matter, but it can be argued that the wall function approach is not for those requiring a very high accuracy - they have to go for the full resolution of the boundary layer.


It is interesting that your computations have shown y+ insensitive results on wall pressure, even when the wall unit is in the buffer layer.