Theory: Turbulence Models

This section details the turbulence models implemented in Caelus 9.04. Caelus 9.04 supports both steady-state 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 Navier-Stokes Equations

The general form of the incompressible Navier-Stokes equations in absence of body forces is given by the following expression:

\[\partial_t u_i + \partial_j(u_i u_j) = -\frac{1}{\rho} \partial_i p + \partial_j \nu \partial_j u_i.\]

The momentum equation is complemented by continuity equation, which for incompressible flow is the divergence free condition for velocity:

\[\partial_i u_i = 0\]

The Reynolds averaging procedure decomposes instantaneous variables into time averaged and fluctuating components according the expression

\[u_i(x_i,t) = \overline{u}_i(x_i) + u_i'(x_i,t).\]

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\)

\[\overline{u}_i = \frac{1}{T} \int u_i(x_i,t) dt.\]

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:

\[\frac{1}{T} \int u_i'(x_i,t) dt = 0.\]

Keeping in mind this property, the Reynolds averaging procedure is applied to the Navier-Stokes equations to obtain Reynolds Averaged Navier-Stokes (RANS) equations:

\[\partial_i \overline{u}_i = 0\]

\[\partial_t \overline{u}_i + \partial_j (\overline{u}_i \overline{u}_j) = \frac{1}{\rho} \partial_i \overline{p} + \partial_j(\nu \partial_j \overline{u}_i ) - \partial_j \overline{u_i' u_j'}\]

The new term in momentum equations, \(-\partial_j\overline{u_i' u_j'}\), appears due to the non-linear 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:

\[-\overline{u_i' u_j'} = 2 \nu_{t} \left( \overline{S}_{ij} - \frac{2}{3} k \delta_{ij}\right),\]

where the term \(\overline{S}_{ij}\) is the symmetric part of the time averaged stress tensor

\[\overline{S}_{ij} = \frac{1}{2}(\partial_j \overline{u}_i + \partial_i \overline{u}_j),\]

and \(k\) is turbulent kinetic energy

\[k = \frac{1}{2} \sum_i \overline{(u_i')^2}.\]

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 over-bar used to denote time averaged quantities, will be omitted for the clarity throughout this section.

Spalart-Allmaras Turbulence Model

Model Description

The Spalart-Allmaras model uses one equation to compute the turbulent viscosity. The Spalart-Allmaras 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 Spalart-Allmaras 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 viscosity-like variable

\[\partial_t \tilde{\nu} + u_j \partial_j \tilde{\nu} = c_{b1}f_{r1}\tilde{S}\tilde{\nu} - c_{w1}f_{w1}\left(\frac{\tilde{\nu}}{d}\right)^2 + \frac{1}{\sigma}\left[\partial_j(\nu+\tilde{\nu})\partial_j \tilde{\nu} + c_{b2}\partial_i \tilde{\nu} \partial_i \tilde{\nu} \right]\]

Terms appearing in the Spalart-Allmaras turbulence model have the following definitions:

Model definitions

Symbol	Definition	Formula
\(\tilde{\nu}\)	Viscosity-like 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^6-r)\)
\(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 Spalart-Allmaras 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 Spalart-Allmaras 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 Spalart-Allmaras 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 Spalart-Allmaras 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 Spalart-Allmaras on meshes that are not wall-resolved, 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. non-realizable. The realizable \(k\)–\(\varepsilon\) turbulence model resolves this problem by introducing a non-constant definition of \(C_{\mu}\) and by modifying the \(\varepsilon\) equation. The transport equations for the realizable \(k\)–\(\varepsilon\) turbulence model are as follows:

\[\partial_t k + \partial_j (u_j k) = \partial_j\left[\mu + \frac{\mu_t}{\sigma_k}\partial_j k\right] + G_k -\varepsilon\]

\[\partial_t \varepsilon + \partial_j (u_j \varepsilon) = \partial_j\left[\mu + \frac{\mu_t}{\sigma_{\varepsilon}}\partial_j \varepsilon \right] + C_1 S \varepsilon - C_2 \frac{\varepsilon^2}{k+ \sqrt{\nu\varepsilon}}\]

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 high-Reynolds 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 free-stream 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 near-wall 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\) free-stream 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:

\[\partial_t k + \partial_j (u_j k) = \frac{P}{\rho} - \beta^* \omega k + \partial_j[(\nu + \sigma_k \nu_t) \partial_j k]\]

\[\partial_t \omega + \partial_j(u_j \omega) = \frac{\gamma P}{\rho \nu_t} - \beta \omega^2 + \partial_j[(\nu + \nu_t) \partial_j \omega] + 2(1-F_1) \frac{\sigma_{\omega2}}{\omega} \partial_j k \partial_j \omega\]

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,

\[P = min(P, 10 \beta^* \omega k).\]

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 near-wall 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 wall-resolved, 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.