CFD Online Logo CFD Online URL
Home > Wiki > Two equation turbulence models

Two equation turbulence models

From CFD-Wiki

Jump to: navigation, search
Turbulence modeling
RANS-based turbulence models
  1. Linear eddy viscosity models
    1. Algebraic models
      1. Cebeci-Smith model
      2. Baldwin-Lomax model
      3. Johnson-King model
      4. A roughness-dependent model
    2. One equation models
      1. Prandtl's one-equation model
      2. Baldwin-Barth model
      3. Spalart-Allmaras model
    3. Two equation models
      1. k-epsilon models
        1. Standard k-epsilon model
        2. Realisable k-epsilon model
        3. RNG k-epsilon model
        4. Near-wall treatment
      2. k-omega models
        1. Wilcox's k-omega model
        2. Wilcox's modified k-omega model
        3. SST k-omega model
        4. Near-wall treatment
      3. Realisability issues
        1. Kato-Launder modification
        2. Durbin's realizability constraint
        3. Yap correction
        4. Realisability and Schwarz' inequality
  2. Nonlinear eddy viscosity models
    1. Explicit nonlinear constitutive relation
      1. Cubic k-epsilon
      2. EARSM
    2. v2-f models
      1. \overline{\upsilon^2}-f model
      2. \zeta-f model
  3. Reynolds stress model (RSM)
Large eddy simulation (LES)
  1. Smagorinsky-Lilly model
  2. Dynamic subgrid-scale model
  3. RNG-LES model
  4. Wall-adapting local eddy-viscosity (WALE) model
  5. Kinetic energy subgrid-scale model
  6. Near-wall treatment for LES models
Detached eddy simulation (DES)
Direct numerical simulation (DNS)
Turbulence near-wall modeling
Turbulence free-stream boundary conditions
  1. Turbulence intensity
  2. Turbulence length scale

Two equation turbulence models are one of the most common types of turbulence models. Models like the k-epsilon model and the k-omega model have become industry standard models and are commonly used for most types of engineering problems. Two equation turbulence models are also very much still an active area of research and new refined two-equation models are still being developed.

By definition, two equation models include two extra transport equations to represent the turbulent properties of the flow. This allows a two equation model to account for history effects like convection and diffusion of turbulent energy.

Most often one of the transported variables is the turbulent kinetic energy, k. The second transported variable varies depending on what type of two-equation model it is. Common choices are the turbulent dissipation, \epsilon, or the Specific turbulence dissipation rate, \omega. The second variable can be thought of as the variable that determines the scale of the turbulence (length-scale or time-scale), whereas the first variable, k, determines the energy in the turbulence.


Boussinesq eddy viscosity assumption

The basis for all two equation models is the Boussinesq eddy viscosity assumption, which postulates that the Reynolds stress tensor, \tau_{ij}, is proportional to the mean strain rate tensor, S_{ij}, and can be written in the following way:

\tau_{ij} = 2 \, \mu_t \, S_{ij} - \frac{2}{3}\rho k \delta_{ij}

Where \mu_t is a scalar property called the eddy viscosity which is normally computed from the two transported variables. The last term is included for modelling incompressible flow to ensure that the definition of turbulence kinetic energy is obeyed:

k=\frac{\overline{u'_i u'_i}}{2}

The same equation can be written more explicitly as:

 -\rho\overline{u'_i u'_j} = \mu_t \, \left( \frac{\partial U_i}{\partial x_j} + \frac{\partial U_j}{\partial x_i} \right) - \frac{2}{3}\rho k \delta_{ij}

The Boussinesq assumption is both the strength and the weakness of two equation models. This assumption is a huge simplification which allows one to think of the effect of turbulence on the mean flow in the same way as molecular viscosity affects a laminar flow. The assumption also makes it possible to introduce intuitive scalar turbulence variables like the turbulent energy and dissipation and to relate these variables to even more intuitive variables like turbulence intensity and turbulence length scale.

The weakness of the Boussinesq assumption is that it is not in general valid. There is nothing which says that the Reynolds stress tensor must be proportional to the strain rate tensor. It is true in simple flows like straight boundary layers and wakes, but in complex flows, like flows with strong curvature, or strongly accelerated or decelerated flows the Boussinesq assumption is simply not valid. This gives two equation models inherent problems to predict strongly rotating flows and other flows where curvature effects are significant. Two equation models also often have problems to predict strongly decelerated flows like stagnation flows.

Near-wall treatments

HRN (left) vs LRN(right). HRN uses log law in order to estimate gradient in the cell.

The structure of turbulent boundary layer exhibits large, compared with the flow in the core region, gradients of velocity and quantities characterising turbulence. See introduction to wall bounded turbulent flows for more detail. In a collocated grid these gradients will be approximated using discretisation procedures which are not suitable for such high variation since they usually assume linear interpolation of values between cell centres.

Moreover, the additional quantities appearing in two-equation models require specification of their own boundary conditions that on purely physical grounds cannot be specified a priori.

This situation gave rise to a plethora of near-wall treatments. Generally speaking, two approaches can be distinguished:

  • Low Reynolds number treatment (LRN)
  • High Reynolds number treatment (HRN)

LRN integrates every equation up to the viscous sublayer and therefore the first computational cell must have its centroid in y^{+}\sim 1. This results in very fine meshes close to the wall. Additionally, for some models additional treatment (damping functions) of equations is required to guarantee asymptotic consistency with the turbulent boundary layer behaviour. This often makes the equations stiff and further increases computation time.

HRN also known as wall functions approach relies on log-law velocity profile and therefore the first computational cell must have its centroid in the log-layer. Use of HRN enhances convergence rate and often numerical stability.

Interestingly, none of the current approaches can deal with buffer layer i.e. the layer in which both viscous and Reynolds stresses are significant. The first computational cell should be either in viscous sublayer or in log-layer -- not in-between. Automatic wall treatments, available in some codes, are an ad hoc solution but the blending techniques employed there are usually arbitrary and though they can achieve the switching between HRN and LRN treatments they cannot be regarded as the correct representation of the buffer layer.

There are two[1] possible ways of implementing wall functions in a finite volume code:

  • Additional source term in the momentum equations.
  • Modification of turbulent viscosity in cells adjacent to solid walls.

The source term in the first approach is simply the difference between logarithmic and linear interpolation of velocity gradient multiplied by viscosity (the difference between shear stresses). The second approach does not attempt to reproduce the correct velocity gradient. Instead, turbulent viscosity is modified in such a way as to guarantee the correct shear stress.

Standard wall functions

Using the compact version of log-law

\frac{U}{u_\tau} = \frac{1}{\kappa} \ln E y^{+}

where E=9.8 is equivalent to additive constants in the logarithmic law of the wall. Now using \tau_w = \rho u_\tau^2 we obtain:

  \tau_w = \frac{\rho u_\tau \kappa U}{\ln Ey^{+}},

On the other hand, the linear interpolation for shear stress, rembembering that U|_{y=0}=0, is:

  \tau_w = \rho(\nu_t + \nu)\frac{U_p}{y_p}.

Comparing the above equations we obtain an expression for turbulent viscosity can be obtained:

  \nu_t = \nu\left( \frac{y^{+}\kappa}{\ln Ey^{+}} - 1\right).

Note that u_\tau has been been incorporated in y^+. The latter remains the only unknown in the equation and has to be estimated for the current velocity field. In the standard approach this cannot be done explicitly and instead an implicit way of obtaining y^{+} has to be employed.

After multiplying log law (1) by y_p/\nu and after reorganising some terms we get:

  \frac{\kappa U_p y_p}{\nu} = y^{+}\ln{Ey^{+}}.

This equation can be solved numerically with respect to y^+ for example via root searching algorithms e.g. Newton method for specified U_p, y_p and \nu. One iteration in a Newton method for (3) is

  y^{+}_{n+1} = \frac{\frac{\kappa U_p y_p}{\nu} + y^{+}_{n}}{1 + \ln E y_n^{+}}.

Thus obtained y^{+} is then substituted to (2) Eventually the estimated u_\tau serves also to define the values of turbulent quantities in the cell adjacent to the wall:

  k_p = \frac{u_\tau^2}{\sqrt{C_\mu}}

and for e.g. k-\omega model:

  \omega_p = \frac{\sqrt{k_p}}{{C_\mu^{1/4}}\kappa y_p},

which are the values for k and \omega resulting from asymptotic analysis[2] of log layer. These wall functions for k and \omega are the results of the solution of model equation for the logarithmic layer.

The above methodology is known to produce spurious results in separated flows, where, by definition, u_\tau = 0 at the separation and reattachment point. Many extension of this approach has been proposed.

Launder-Spalding wall functions

This approaches utilises the relation

  u_\tau = C_\mu^{1/4} \sqrt{k}.

A three step procedure[3] is adopted:

  • Calculate wall viscosity.
  • Modify the production term in k equation in the wall adjacent cell.
  • Specify the value of the second variable in the wall adjacent cell.

After substituing (4) to log law relation we obtain:

  U = \frac{C_\mu^{1/4}\sqrt{k}}{\kappa}\ln E y^*,


  y^* = \frac{C_\mu^{1/4}\sqrt{k} y}{ \nu},

The production term (just in the first cell!) can be approximated as follows:

  {P}_p = -\overline{u_1^'u_2^'}\frac{\partial U_1}{\partial x_2} = 
  \underbrace{(\nu_t + \nu)\frac{\partial U_1}{\partial x_2}}_{\text{Boussinesq assumption}}
  \overbrace{\frac{C_\mu^{1/4} \sqrt{k_p}}{\kappa y_p}}^{\text{Logarithmic velocity}}

Lastly, we specify the remaining turbulence quantities in the same manner as in the previous approach.

See Also


  1. Bredberg, J. (2000), On the Wall Boundary Condition for Turbulence Models, Chalmers University of Technology Goteborg, Sweden
  2. Wilcox, D.C. (1993), Turbulence Modeling for CFD, DCW Industries, Inc.
  3. Launder, B.E. & D.B. Spalding (1974), "The numerical computation of turbulent flows", Computer Methods in Applied Mechanics and Engineering 3 (2): 269-289, ISSN 00457825

My wiki