Theory: Turbulence Models¶
This section details the turbulence models implemented in Caelus 9.04. Caelus 9.04 supports both steadystate and transient turbulence models. From an averaging point of view, Caelus 9.04 supports Reynolds averaged (RANS) and large eddy simulations (LES). The current focus of turbulence modelling in Caelus 9.04 is on incompressible flows. Future versions of Caelus will support compressible flows and appropriate compressible turbulence models.
Reynolds Averaged NavierStokes Equations¶
The general form of the incompressible NavierStokes equations in absence of body forces is given by the following expression:
The momentum equation is complemented by continuity equation, which for incompressible flow is the divergence free condition for velocity:
The Reynolds averaging procedure decomposes instantaneous variables into time averaged and fluctuating components according the expression
Instantaneous quantities are averaged over a sufficiently long time such that the time scale of the fluctuating part of the field is much smaller than the averaging interval \(T\)
This definition allows for time variability in the averaged momentum and continuity equations although the structures that have characteristic time scales comparable to or larger than the filter width \(T\) could be resolved.
An important property of the fluctuating portion of the flow fields in Reynolds averaged simulations is that they have zero mean value by the definition:
Keeping in mind this property, the Reynolds averaging procedure is applied to the NavierStokes equations to obtain Reynolds Averaged NavierStokes (RANS) equations:
The new term in momentum equations, \(\partial_j\overline{u_i' u_j'}\), appears due to the nonlinear nature of the momentum equation. This term is called the Reynolds stress tensor and it represents the influence of turbulence on the mean flow. The Reynolds stress tensor contains 6 unique, unknown components. The unknowns must be modelled and are determined by relating them to time averaged quantities. The most common approach used to close the Reynolds stresses is the Boussinesq hypothesis:
where the term \(\overline{S}_{ij}\) is the symmetric part of the time averaged stress tensor
and \(k\) is turbulent kinetic energy
The Boussinesq hypothesis is also called eddy viscosity approximation and most RANS turbulence models use it in this linear form. The main difference between various turbulence models employing the Boussinesq hypothesis is how the turbulent viscosity \(\nu_t\) is determined.
The following sections describes the transport equations of the RANS turbulence models included in Caelus 9.04. Note, the overbar used to denote time averaged quantities, will be omitted for the clarity throughout this section.
SpalartAllmaras Turbulence Model¶
Model Description
The SpalartAllmaras model uses one equation to compute the turbulent viscosity. The SpalartAllmaras model used in Caelus 9.04 corresponds to the version of Shur et al. [9] with curvature correction. The transport equation in this version is the same as the original SpalartAllmaras model [10] but with rotation and curvature effects accounted for. This modification is accomplished by introducing the \(f_{r1}\) function into the transport equation for the viscositylike variable
Terms appearing in the SpalartAllmaras turbulence model have the following definitions:
Model definitions
Symbol 
Definition 
Formula 

\(\tilde{\nu}\) 
Viscositylike variable 

\(\nu_t\) 
Turbulent eddy viscosity 
\(\tilde{\nu} f_{v1}\) 
\(f_{v1}\) 
Viscous Damping Function 
\(\frac{\chi^3}{\chi^3 + c_{v1}}\) 
\(\chi\) 
Viscosity ratio 
\(\frac{\tilde{\nu}}{\nu}\) 
\(\nu\) 
Kinematic Viscosity 

\(\tilde{S}\) 
Modified vorticity 
\(\Omega + \frac{\tilde{\nu}}{\kappa^2 d^2}f_{v2}\) 
\(\Omega\) 
Vorticity magnitude 
\(\Omega = \sqrt{2W_{ij}W_{ij}}\) 
\(W_{ij}\) 
Vorticity 
\(\frac{1}{2}(\partial_j u_i  \partial_i u_j)\) 
\(d\) 
Wall distance 

\(f_{v2}\) 
Function 
\(1\frac{\chi}{1+\chi f_{v1}}\) 
\(f_{w}\) 
Function 
\(g\left[\frac{1+c_{w3}^6}{g^5+c_{w3}^6}\right]^{1/6}\) 
\(g\) 
Function 
\(r+c_{w2}(r^6r)\) 
\(r\) 
Variable 
\(r=min\left[\frac{\tilde{\nu}}{\tilde{S} \kappa^2 d^2},10\right]\) 
It should be noted that the Caelus 9.04 implementation of the SpalartAllmaras turbulence model omits the rotational and turbulence tripping terms as the flow is assumed fully turbulent.
Model constants
\(c_{b1} = 0.1355\)
\(\sigma = \frac{2}{3}\)
\(c_{b2} = 0.622\)
\(c_{w2} = 0.3\)
\(c_{w3} = 2.0\)
\(c_{v1} = 7.1\)
\(c_{w1} = \frac{c_{b1}}{\kappa^2} + \frac{1+c_{b2}}{\sigma}\)
\(\kappa = 0.41\)
The SpalartAllmaras model neglects the contribution of the turbulent kinetic energy to the eddy viscosity. This is deemed acceptable for many incompressible and low Mach number compressible flows.
Boundary conditions
The SpalartAllmaras turbulence model was originally developed for external aerodynamics of streamlined bodies. Typically, mesh resolutions used for external aerodynamics have \(y^+\) values close to 1. Therefore, the SpalartAllmaras model is often referred to as a low Reynolds number model when it comes to near wall treatment. In the case of meshes resolved near the wall, i.e. \(y^+ \sim 1\), the following set of boundary conditions are recommended:
\(\tilde{\nu}_{wall} = 0\)
\(\nu_t_{wall} = 0\)
However, in Caelus 9.04 it is possible to use the SpalartAllmaras on meshes that are not wallresolved, i.e., meshes that have the first cell height next to the wall \(y^+ \gg 1\). In such cases, a wall function for \(\nu_t\) based on the velocity is recommended.
\(\nu_t_{wall} =\) nutUWallFunction (\(y^+ > 30\))
Inlet values of \(\tilde{\nu}\) at the far field boundary for external flows is suggested as follows:
\(\tilde{\nu}_{inlet} = 3 \nu\) to \(5 \nu\).
\(\nu_t\) should be set to “calculated”.
It is also possible to specify \(\tilde{\nu}\) from the turbulence intensity and length scale as follows:
\(\tilde{\nu}_{inlet} = \sqrt{\frac{3}{2}} u_i I l\)
where \(I = \frac{u_i'}{u_i}\) is turbulence intensity and
\(l\) is a turbulence length scale (measure of largest eddies in the flow field).
In the case of internal flows, depending on the level of turbulence intensity, \(\tilde{\nu}\) can be determined from the turbulence viscosity ratio
\(\tilde{\nu} = 10\nu\) or higher and
\(\nu_t\) should be set to “calculated”.
At the outlet, zero gradient an acceptable choice.
Realizable \(k\)–\(\varepsilon\) Turbulence Model¶
Model description
The realizable \(k\)–\(\varepsilon\) turbulence model is an improvement on the standard \(k\)–\(\varepsilon\) model as it is modified to guarantee realizability of normal Reynolds stresses [13]. In the standard \(k\)–\(\varepsilon\) turbulence model it is possible, under certain conditions, to have negative normal Reynolds stresses, i.e. nonrealizable. The realizable \(k\)–\(\varepsilon\) turbulence model resolves this problem by introducing a nonconstant definition of \(C_{\mu}\) and by modifying the \(\varepsilon\) equation. The transport equations for the realizable \(k\)–\(\varepsilon\) turbulence model are as follows:
The terms appearing in the realizable \(k\)–\(\varepsilon\) turbulence model have the following definitions:
Model definitions
Symbol 
Definition 
Formula 

\(\nu_t\) 
Turbulent viscosity 
\(C_{\mu}\frac{k^2}{\varepsilon}\) 
\(C_1\) 
Model parameter 
\(max \left(0.43,\frac{\eta}{\eta+5}\right)\) 
\(\eta\) 
Strain magnitude over time scale 
\(S\frac{k}{\varepsilon}\) 
\(S\) 
Strain magnitude 
\(S = \sqrt{2S_{ij}S_{ij}}\) 
\(S_{ij}\) 
Strain tensor 
\(\frac{1}{2}(\partial_j u_i + \partial_i u_j)\) 
\(C_{\mu}\) 
Model parameter 
\(\frac{1}{A_0 + A_S \frac{k U^*}{\varepsilon}}\) 
\(U^*\) 
Velocity 
\(\sqrt{S_{ij}S_{ij}+\tilde{\Omega_{ij}}\tilde{\Omega_{ij}}}\) 
\(\overline{\Omega}_{ij}\) 
Rate of rotation tensor 
\(\frac{1}{2}(\partial_j u_i  \partial_i u_j)\) 
\(\omega_k\) 
Angular velocity 

\(\Omega_{ij}\) 
Rotation tensor 
\(\overline{\Omega}_{ij}  \varepsilon_{ijk} \omega_k\) 
\(\tilde{\Omega}_{ij}\) 
Modified rotation tensor 
\(\Omega_{ij}  2 \varepsilon_{ijk} \omega_k\) 
\(A_0\) 
Model constant 
\(4.04\) 
\(A_S\) 
Function 
\(\sqrt{6} cos(\phi)\) 
\(\phi\) 
Function 
\(\frac{1}{3} cos^{1} (\sqrt{6} W)\) 
\(W\) 
Variable 
\(\frac{S_ijS_{jk}S{ki}}{S^3}\) 
Model constants
\(C_{1\varepsilon} = 1.44\)
\(C_2 = 1.9\)
\(\sigma_k = 1.0\)
\(\sigma_{\varepsilon} = 1.2\)
Boundary conditions
The realizable \(k\)–\(\varepsilon\) was developed as a highReynolds number model. In other words, the wall treatment is such that wall functions for \(\nu_t\), \(k\), and \(\varepsilon\) fields must be always used. It is recommended that velocity based wall functions for the \(\nu_t\) field is used. On the inlet, the following relations can be used to determine the values of \(k\) and \(\varepsilon\) fields:
\(k_{inlet} = \frac{3}{2} (U_{inlet} I)^2\),
\(\varepsilon = C_{\mu}^{3/4} \frac{k^{3/2}}{l}\),
\(\nu_t\) should be set to “calculated” or computed from the viscosity ratio.
\(I = \frac{u_i'}{u_i}\) is the turbulence intensity
and \(l\) is a turbulence length scale.
On outlets, zero gradient is an acceptable choice.
\(k\)–\(\omega\) SST Model¶
Model description
The \(k\)–\(\omega\) shear stress transport (SST) turbulence model was introduce by Menter [16] as an improvement on the original two equation \(k\)–\(\omega\) model of Wilcox [14]. The SST formulation uses the \(k\)–\(\omega\) formulation in the inner parts of the boundary layer and it switches to a \(k\)–\(\varepsilon\) behaviour in the freestream portions of the flow domain. Thus, the \(k\)–\(\omega\) SST turbulence model does not require a damping function close to the wall enabling its use with nearwall resolved meshes without modification. At the same time, the \(k\)–\(\omega\) SST models reverts to a \(k\)–\(\varepsilon\) model behaviour away from the walls thus removing the sensitivity of the original \(k\)–\(\omega\) model to the \(\omega\) freestream boundary value.
The implementation of the \(k\)–\(\omega\) SST model in Caelus 9.04 is that of Kuntz et al [3]. The transport equations for turbulence kinetic energy and specific dissipation are:
It should be noted that this set of transport equations only applicable for incompressible flows as density was assumed constant. Terms appearing in the \(k\)–\(\omega\) SST model have the following definitions:
Model definitions
Symbol 
Definition 
Formula 

\(\nu_t\) 
Turbulent viscosity 
\(\frac{a_1 k}{max (a_1 \omega, S F_2)}\) 
\(S_{ij}\) 
Symmetric part of stress tensor 
\(\frac{1}{2}(\partial_j u_i + \partial_i u_j)\) 
\(S\) 
Strain invariant 
\(S=\sqrt{S_{ij}S_{ij}}\) 
\(F_1\) 
Blending function 
\(tanh(arg_1^4)\) 
\(arg_1\) 
Argument for \(F_1\) function 
\(min\left[max \left(\frac{\sqrt{k}}{\beta^* \omega d},\frac{500 \nu}{d^2 \omega}\right), \frac{4 \sigma_{\omega 2} k}{CD_{k \omega} d^2}\right]\) 
\(CD_{k \omega}\) 
Function 
\(max\left(2 \sigma_{\omega 2}\frac{1}{\omega} \frac{\partial k}{\partial x_j} \frac{\partial \omega}{\partial x_j}, 10^{10}\right)\) 
\(F_2\) 
Blending function 
\(tanh(arg_2^2)\) 
\(arg_2\) 
Argument for \(F_2\) function 
\(max \left(2 \frac{\sqrt{k}}{\beta^* \omega d},\frac{500 \nu}{d^2 \omega}\right)\) 
This variant of the \(k\)–\(\omega\) SST turbulence model uses a production limiter in both \(k\) and \(\omega\) equations,
Model constants
\(\gamma_1 = \frac{5}{9}\)
\(\gamma2 = 0.44\)
\(\sigma_{k 1} = 0.85\)
\(\sigma_{\omega 1} = 0.5\)
\(\beta_1 = 0.075\)
\(\sigma_{k 2} = 1.0\)
\(\sigma_{\omega 2} = 0.856\)
\(\beta_2 = 0.0828\)
\(\beta^* = 0.09\)
\(\kappa = 0.41\)
\(a_1 = 0.31\)
Boundary conditions
As described above, the \(k\)–\(\omega\) SST turbulence model does not require damping functions near the wall. Thus it can be used for both resolved and unresolved nearwall meshes. In the case of meshes resolved near the wall, i.e. \(y^+ \sim 1\), the following set of boundary conditions are recommended:
\(k_{wall} = 0\)
\(\omega_{wall} =\) omegaWallFunction
\(\nu_t_{wall} = 0\)
For meshes that are not wallresolved, i.e., meshes that have the first cell height next to the wall \(y^+ \gg 1\), wall functions will be needed. In such cases, a wall function for \(\nu_t\) based on the velocity is recommended.
\(k_{wall} =\) kqRWallFunction
\(\omega_{wall} =\) omegaWallFunction
\(\nu_t_{wall} =\) nutUSpaldingWallFunction or alternatively nutUWallFunction if \(y^+ \gt 30\)
Inlet values of can be determined as follows:
\(k_{inlet} = \frac{3}{2} (U_{inlet} I)^2\),
\(\omega = C_{\mu}^{1/4} \frac{\sqrt{k}}{l}\),
\(\nu_t\) should be set to “calculated” or computed from the turbulent viscosity relation.
\(I = \frac{u_i'}{u_i}\) is the turbulence intensity
and \(l\) is a turbulence length scale.
On outlets, zero gradient is an acceptable choice.