Turbulence Modeling

AcuSolve supports a variety of turbulence models for fluid flow simulations. For steady state simulations the Reynolds-averaged Navier-Stokes equations are solved to arrive directly at the time averaged flow field. In case of transient simulations the governing equations are integrated in time to yield an accurate description of the flow field. There are also many different turbulence closure methods available for each type of simulation.

The various turbulence models used in AcuSolve are:

One Equation and Two Equation Models: Spalart-Allmaras (SA) turbulence model; Shear Stress Transport (SST) turbulence model; $k - ω$ turbulence model; $k - ε$ model; RNG $k - ε$ model; Realizable $k - ε$ model
Detached Eddy Simulation (DES) Models: SA – DES; SST – DES; Dynamic DES (DDES)
Large Eddy Simulation (LES) Models: Classical (Smaroginsky) model; Dynamic subgrid LES model

Reynolds-Averaged Navier-Stokes (RANS) Simulations

For incompressible turbulent flow the instantaneous velocity field can be decomposed into a time averaged velocity and its corresponding fluctuation.

This is commonly called the Reynolds’ decomposition, which is shown in Figure 1.

u_{i} = \bar{u_{i}} + u_{i}^{'}

where

$u_{i}$ : instantaneous velocity,
$\bar{u_{i}} = \frac{1}{Δ t} \int_{t}^{t + Δ t} u_{i} d t$ : time averaged velocity over $Δ t$ ,
$u_{i}^{'}$ velocity fluctuation.

Figure 1. Turbulent Velocity Signals as a Function of Time

Similarly, instantaneous pressure can be decomposed into the time averaged and fluctuation terms.

p = \bar{p} + p^{'}

where $p$ is the instantaneous pressure.

Once this concept is substituted into the instantaneous Navier-Stokes equations and then time averaging is performed the RANS equations are obtained. The following equations are, respectively, the Reynolds-averaged continuity and Reynolds-averaged momentum equations:

\frac{\partial ρ}{\partial t} + \frac{\partial (ρ {\bar{u}}_{j})}{\partial x_{j}} = 0

\frac{\partial (ρ {\bar{u}}_{i})}{\partial t} + \frac{\partial (ρ {\bar{u}}_{i} {\bar{u}}_{j})}{\partial x_{j}} = - \frac{\partial (\bar{p} δ_{i j})}{\partial x_{j}} + 2 μ \frac{\partial \bar{S_{i j}}}{\partial x_{j}} + \frac{\partial τ_{i j}^{R}}{\partial x_{j}}

where

$\bar{S_{i j}} = \frac{1}{2} (\frac{\partial \bar{u_{i}}}{\partial x_{j}} + \frac{\partial \bar{u_{j}}}{\partial x_{i}})$ is the mean strain rate tensor,
$τ_{i j}^{R} = - ρ \bar{u_{i}^{'} u_{j}^{'}}$ is the Reynolds stress tensor.

It is noted that the momentum equation has an unsteady term (the first term in the left hand side of the momentum equation). Wegner et al. (2004) argued that the RANS equations with the time term (unsteady RANS) can be applied to simulate unsteady flows when a time scale (T2) for large scale motions is larger than a time averaging period for RANS (T1) (see Figure 2). It is worth mentioning that RANS time averaging period (T1) should be considerably larger than the largest turbulence (integral) time scale, as URANS can only give a time averaged mean value for the velocity field, not solving turbulence. Thus, a representative URANS solution is demonstrated by the smoothed red line in the image below. Whereas the blue line represents the velocity profile with instantaneous fluctuations.

Steady RANS can be obtained by removing the time term when there is no large scale unsteadiness observed in turbulent flow.

Figure 2. Time Scales for T2 for Large Unsteady Oscillation vs T1 for Turbulent Time Scale

Although Reynolds-averaging eliminates a need to compute the instantaneous flow field, it introduces a new unknown term in comparison to the Navier-Stokes equations. This unknown term is referred to as the Reynolds stress tensor, which represents an added stress due to the turbulent motions. This term arises from the decomposition of the convective term using the Reynolds decomposition. Since your goal is to eliminate any dependence on the instantaneous flow field in these simplified equations, this term should be modeled in terms of the mean flow. This challenge is known as the closure problem of turbulence.

One method that is available for computing the Reynolds Stresses that appear in the RANS equations is to use the Boussinesq approximation. This approach assumes a linear relationship between the turbulent Reynolds stresses and the mean rate of strain tensor:

τ_{i j}^{R} = - ρ \bar{u_{i}^{'} u_{j}^{'}} = μ_{t} [\frac{\partial \bar{u_{i}}}{\partial x_{j}} + \frac{\partial \bar{u_{j}}}{\partial x_{i}} - \frac{2}{3} \frac{\partial \bar{u_{k}}}{\partial x_{k}} δ_{i j}] - \frac{2}{3} ρ k δ_{i j}

where

$μ_{t}$ : the eddy viscosity,
$k = \frac{1}{2} \bar{u_{i}^{'} u_{j}^{'}}$ : the turbulent kinetic energy per unit mass,
$δ_{i j}$ : the Kronecker delta.

Since $\frac{\partial \bar{u_{k}}}{\partial x_{k}}$ = 0 for incompressible flow the isotropic eddy viscosity becomes

τ_{i j}^{R} = μ_{t} [\frac{\partial \bar{u_{i}}}{\partial x_{j}} + \frac{\partial \bar{u_{j}}}{\partial x_{i}}] - \frac{2}{3} ρ k δ_{i j}

The only unknown in this equation is the eddy viscosity, which must be determined by the turbulence model. There are hundreds of available turbulence models that have been developed to provide this eddy viscosity. Some common Boussinesq models include

Mixing length model
Spalart-Allmaras model
$k - ε$ model, RNG (Renormalized Group) $k - ε$ model, Realizable $k - ε$ model
$k - ω$ model, SST (Shear Stress Transport) model
$v 2 f$ model, zeta-f model

This approach is valid for many turbulent flows, for example, boundary layer, channel flow, round jets, turbulent shear flow and mixing layer, but not for all, for example, impinging flow and swirling flow. Under these circumstances, a model that is capable of predicting anisotropic Reynolds stresses is more appropriate. This can include a full Reynolds stress model (RSM), which determines the turbulent stresses by solving a transport equation for each stress component, or an eddy viscosity model that uses a non linear stress strain relationship to compute each Reynolds stress component separately. Although the improvement in accuracy achieved from this type of modeling is appealing, the introduction of anisotropic stresses is often accompanied by increased compute cost and a decrease in robustness.

The following section provides a review of some commonly used RANS models, ranging from algebraic models (zero equation models) up to the seven equation Reynolds Stress Model.

General Form of Turbulence Models

The general formation of turbulence models can be written as

\underset{Unsteady term}{\underset{}{\frac{\partial (ρ ϕ)}{\partial t}}} + \underset{Convection term}{\underset{}{\frac{\partial (ρ {\bar{u}}_{j} ϕ)}{\partial x_{j}}}} = \underset{Diffusion term}{\underset{}{\frac{\partial}{\partial x_{j}} [(μ + \frac{μ_{t}}{σ}) \frac{\partial ϕ}{\partial x_{j}}]}} + \underset{Production term}{\underset{}{P_{ϕ}}} + \underset{Dissipation term}{\underset{}{D_{ϕ}}}

where $ϕ$ is a turbulence model variable, $μ_{t}$ is the eddy viscosity, $μ$ is the material viscosity and $σ$ is a constant.

The left-hand side of the equation describes advection, which includes two terms, the acceleration term and convective term. The right-hand side of the equation represents the summation of the diffusion, production and dissipation terms. The unsteady term represents the time dependence of turbulence model variables, while the convective term accounts for the rate of change of variables due to convection by the mean flow. The diffusion term on the right side of the equation describes the transport of turbulent variables due to the summation of material viscosity and eddy viscosity. The production term indicates the production rate of turbulent variables from the mean flow gradient, while the dissipation term represents the dissipation rate of turbulent variables due to viscous stresses.

One Equation Eddy Viscosity Models

One equation Reynolds-averaged Navier-Stokes (RANS) models solve a single scalar transport equation to compute the eddy viscosity.

Common one equation models include the Spalart-Allmaras (SA) model and the Nut-92 model. The SA model is discussed in this manual due to its application in general purpose CFD codes and its popularity for the simulation of external flows and internal flows. The details of the Nut-92 model can be found in Shur et al. (1995).

Spalart-Allmaras (SA) Model

The SA model uses a transport equation to solve for a modified kinematic eddy viscosity, $\hat{υ}$ , as a function of the kinematic eddy viscosity ( $ν_{t}$ ) (Spalart and Allmaras, 1992).

In this model, a length scale (d) in a dissipation term of the modified kinematic eddy viscosity transport equation is specified to determine the dissipation rate. This model has an advantage of having economic solutions for attached flows and moderately separated flows, but it is not recommended for massively separated flows, free shear flows and decaying turbulence.

Transport Equations

\frac{\partial (ρ \hat{υ})}{\partial t} + \frac{\partial (ρ \hat{υ} \bar{u_{j}})}{\partial x_{j}} = \frac{1}{σ} (\frac{\partial}{\partial x_{j}} [(μ + ρ \hat{υ}) \frac{\partial \hat{υ}}{\partial x_{j}}] + ρ C_{b 2} \frac{\partial \hat{υ}}{\partial x_{j}} \frac{\partial \hat{υ}}{\partial x_{j}}) + P + D

where $σ$ and $C_{b 2}$ are constants and $μ$ is the fluid dynamic viscosity. P and D are the production term and destruction term of the modified turbulent viscosity, respectively.

Production of $\hat{υ}$

P = ρ C_{b 1} \hat{S} \hat{υ}

where

$\hat{S} = \sqrt{2 Ω_{i j} Ω_{i j}} + \frac{\hat{υ}}{κ^{2} d^{2}} f_{v 2}$ ,
$Ω_{i j} = \frac{1}{2} (\frac{\partial \bar{u_{i}}}{\partial x_{j}} - \frac{\partial \bar{u_{j}}}{\partial x_{i}})$ is the rotation tensor,
$f_{v 2} = 1 - \frac{χ}{1 + χ f_{v 1}}$ ,
$d$ is the distance from the nearest wall,
$κ$ is the Von Kármán constant,
$C_{b 1}$ is a constant.

Destruction of $\hat{υ}$

D = - ρ C_{w 1} f_{w} {(\frac{\hat{υ}}{d})}^{2}

where

$f_{w} = g {(\frac{1 + C_{w 3}^{6}}{g^{6} + C_{w 3}^{6}})}^{1 / 6}$
$f_{t 2} = C_{t 3} e x p (- C_{t 4} χ^{2})$
$g = r + C_{w 2} (r^{6} - r)$
$r = \frac{\hat{υ}}{\hat{S} κ^{2} d^{2}}$
$C_{w 1}$ is a constant.

Modeling of Turbulent Viscosity $μ_{t}$

The kinematic eddy viscosity for the Spalart-Allmaras model is computed using the following relationship

μ_{t} = ρ \hat{υ} f_{v 1}

where

$f_{v 1} = \frac{χ^{3}}{χ^{3} + C_{v 1}^{3}}$ is the viscous damping function.
$χ = \frac{\hat{υ}}{ν}$ .

Model Coefficients

$C_{w 1} = \frac{C_{b 1}}{κ^{2}} + \frac{1 + C_{b 1}}{σ^{2}}$ , $C_{w 2} = 0.3$ , $C_{w 3} = 2.0$ , $C_{b 1} = 0.1355$ , $C_{b 2} = 0.622$ , $C_{v 1} = 7.1$ $σ$ = $\frac{2}{3}$ , $κ = 0.41$ .

Spalart-Allmaras (SA) Model with Rotation/Curvature Correction

The effects of system rotation and streamline curvature are present in turbomachinery components. Some examples include axial turbines, radial turbines, axial fans, compressors and centrifugal impellors.

This is where most linear eddy viscosity models fail. In order to incorporate the rotational and curvature effects for the SA model, Shur et al. (2000) introduced a version with the modified production term of the transport equation by multiplying the rotation function $f_{r 1}$ .

Modified Production of $\hat{υ}$

P = ρ f_{r 1} C_{b 1} \hat{S} \hat{υ}

where

$f_{r 1} = (1 + C_{r 1}) \frac{2 r *}{1 + r *} (1 - C_{r 3} t a n^{- 1} [C_{r 2} \hat{r}]) - C_{r 1}$
$r * = \frac{S}{\hat{ω}}$
$S = \sqrt{2 S_{i j} S_{i j}}$ is the strain rate magnitude.
$S_{i j} = \frac{1}{2} (\frac{\partial \bar{u_{i}}}{\partial x_{j}} + \frac{\partial \bar{u_{j}}}{\partial x_{i}})$ is the strain rate tensor.
$\hat{ω} = \sqrt{2 {\hat{ω}}_{i j} {\hat{ω}}_{i j}}$ is the modified vorticity magnitude.
${\hat{ω}}_{i j} = \frac{1}{2} ([\frac{\partial \bar{u_{i}}}{\partial x_{j}} - \frac{\partial \bar{u_{j}}}{\partial x_{i}}] + 2 ε_{m j l} {Ω'}_{m})$ is the rotation tensor.
$\hat{r} = \frac{2 ω_{i k} S_{j k}}{D^{4}} (\frac{\partial S_{i j}}{\partial t} + \bar{u_{i j}} \frac{\partial S_{i j}}{\partial x_{j}} + [ε_{i m n} S_{j n} + ε_{j m n} S_{i n}] {Ω'}_{m})$
$D = \sqrt{\frac{1}{2} (S^{2} + {\hat{ω}}^{2})}$
$Ω'$ is the rotation rate.

Model Coefficients

$C_{r 1}$ =1.0, $C_{r 2}$ =12.0, $C_{r 3}$ =1.0

Two Equation Eddy Viscosity Models

Two equation turbulence models are perhaps the most commonly used in industrial flow applications.

Although there are many two equation models available, some of the more popular ones include:

$k - ε$ model, Renormalized Group (RNG) $k - ε$ model, Realizable $k - ε$ model,
$k - ω$ model, Shear Stress Transport (SST) model.

Standard k-ε Model

Launder and Spalding (1974) proposed the standard k-ε turbulence model utilizing the relationships described below.

The model has been shown to be relatively accurate for high Reynolds number flows in which the turbulence behavior is close to homogeneous, and the turbulence production is nearly balanced by dissipation. Nevertheless, its performance deteriorates when predicting boundary layers under adverse pressure gradients (Bradshaw, 1997). It also has difficulty in predicting the viscous sublayer. To resolve this issue, it is suggested that a correction be made to reproduce the law of the wall for incompressible flat-plate boundary layers (Wilcox, 2000). For low Reynolds number flows, the difference between the turbulent kinetic energy production and the dissipation rates may depart from their equilibrium value of zero, thus, ad-hoc adjustments of empirical parameters are inevitable. Furthermore, the standard k-ε turbulence model does not perform well in shear layers and jets, where the turbulent kinetic energy is not balanced with the dissipation rates (Versteeg and Malalasekera, 2007). Finally, this model is not recommended for high swirling/curvature flows, diverging passage flows as well as flows with a body force under the influence of a rotating reference frame.

Transport Equations

Turbulent Kinetic Energy k

\frac{\partial (ρ k)}{\partial t} + \frac{\partial (ρ \bar{u_{j}} k)}{\partial x_{j}} = \frac{\partial}{\partial x_{j}} [(μ + \frac{μ_{t}}{σ_{k}}) \frac{\partial k}{\partial x_{j}}] + P_{k} + D_{k}

Turbulent Dissipation Rate ε

\frac{\partial (ρ ε)}{\partial t} + \frac{\partial (ρ \bar{u_{j}} ε)}{\partial x_{j}} = \frac{\partial}{\partial x_{j}} [(μ + \frac{μ_{t}}{σ_{ε}}) \frac{\partial ε}{\partial x_{j}}] + P_{ε} + D_{ε}

Production Modeling

Turbulent Kinetic Energy k

P_{k} = μ_{t} S^{2}

where

$S = \sqrt{2 S_{i j} S_{i j}}$ is the strain rate magnitude.
$μ_{t}$ is dynamic eddy viscosity.

Turbulent Dissipation Rate ε

P_{ε} = C_{ε 1} \frac{ε}{k} μ_{t} S^{2} = C_{ε 1} \frac{ε}{k} P_{k}

Dissipation Modeling

Turbulent Kinetic Energy k

D_{k} = - ρ ε

Turbulent Dissipation Rate ε

D_{ε} = - C_{ε 2} ρ \frac{ε^{2}}{k}

Modeling of Turbulent Viscosity $μ_{t}$

μ_{t} = C_{μ} \frac{k^{2}}{ε}

Model Coefficients

$C_{ε 1}$ =1.44, $C_{ε 2}$ =1.92, $C_{μ}$ =0.09, $σ_{k}$ =1.0, $σ_{ε}$ =1.3.

Renormalization Group (RNG) k-ε Model

The RNG k-ε turbulence model (Yakhot and Orszag, 1986) deduces the behavior of large scale eddies from that of the smaller ones by utilizing the scale similarity properties that are inherent in the energy cascade (Bradshaw, 1997).

This model employs a modified coefficient in the dissipation rate equation to account for the interaction between the turbulent dissipation and mean shear. It results in better prediction of flows containing high streamline curvature, flows over a backward facing step (Yakhot et al, 1992), and flows in an expansion duct than the standard k-ε turbulence model. However, its performance worsens when predicting flows in a contraction duct (Hanjalic, 2004).

Transport Equations

Turbulent Kinetic Energy k

\frac{\partial (ρ k)}{\partial t} + \frac{\partial (ρ \bar{u_{j}} k)}{\partial x_{j}} = \frac{\partial}{\partial x_{j}} [(μ + \frac{μ_{t}}{σ_{k}}) \frac{\partial k}{\partial x_{j}}] + P_{k} + D_{k}

Turbulent Dissipation Rate ε

\frac{\partial (ρ ε)}{\partial t} + \frac{\partial (ρ \bar{u_{j}} ε)}{\partial x_{j}} = \frac{\partial}{\partial x_{j}} [(μ + \frac{μ_{t}}{σ_{ε}}) \frac{\partial ε}{\partial x_{j}}] + P_{ε} + D_{ε}

Production Modeling

Turbulent Kinetic Energy k

P_{k} = μ_{t} S^{2}

Turbulent Dissipation Rate ε

P_{ε} = C_{ε 1} \frac{ε}{k} μ_{t} S^{2} = C_{ε 1} \frac{ε}{k} P_{k}

Dissipation Modeling

Turbulent Kinetic Energy k

D_{k} = - ρ ε

Turbulent Dissipation Rate ε

D_{ε} = - C'_{ε 2} ρ \frac{ε^{2}}{k}

Modeling of Turbulent Viscosity $μ_{t}$

μ_{t} = C_{μ} \frac{k^{2}}{ε}

Model Coefficients

$C_{ε 1}$ =1.44, $C_{ε 2}$ =1.92, $C_{μ}$ =0.09, $σ_{k}$ =1.0, $σ_{ε}$ =1.3, $C'_{ε 2} = \tilde{C_{ε 2}} + \frac{C_{μ} λ^{3} (1 - λ / λ_{0})}{1 + β λ^{3}}$ , $\tilde{C_{ε 2}}$ = 1.68, $C_{μ}$ = 0.085, $λ$ = $\frac{k}{ε} S$ , $λ_{0}$ = 4.38, $β$ = 0.012, $σ_{k}$ = 0.72, $σ_{ε}$ = 0.72.

Realizable k-ε Model

The standard k-ε turbulence model and RNG k-ε turbulence model do not satisfy mathematical constraints on the Reynolds stresses for the consistency with physics of turbulence.

These constraints include the positivity of normal stresses ( $\bar{u_{α}^{'} u_{α}^{'}} \geq 0$ ) and Schwarz's inequality ( $\bar{u_{α}^{'} u_{α}^{'}} \bar{u_{β}^{'} u_{β}^{'}} \geq \bar{u_{α}^{'} u_{β}^{'}} \bar{u_{α}^{'} u_{β}^{'}}$ ). The realizable k-ε model proposed by Shih et al. (1995) employs a formulation of turbulent viscosity with a variation of $C_{μ}$ to satisfy these constraints. This is why this model is referred to as realizable. The second improvement made with the realizable model is to have a source term in the dissipation rate transport equation, which is derived from the transport equation of the mean square vorticity fluctuation. Compared to the standard k-ε turbulence model, it performs better when simulating planar jet flows, round jet flows and flows under adverse pressure gradients. Although the realizable model has significant performance improvement over the standard model, it has a problem in predicting impinging flow, swirling flow and secondary flows in a square duct since this model is based on the Boussinesq assumption. In addition, this model needs additional functions to simulate the near wall effects.

Transport Equations

Turbulent Kinetic Energy k

\frac{\partial (ρ k)}{\partial t} + \frac{\partial (ρ \bar{u_{j}} k)}{\partial x_{j}} = \frac{\partial}{\partial x_{j}} [(μ + \frac{μ_{t}}{σ_{k}}) \frac{\partial k}{\partial x_{j}}] + P_{k} + D_{k}

Turbulent Dissipation Rate ε

\frac{\partial (ρ ε)}{\partial t} + \frac{\partial (ρ \bar{u_{j}} ε)}{\partial x_{j}} = \frac{\partial}{\partial x_{j}} [(μ + \frac{μ_{t}}{σ_{ε}}) \frac{\partial ε}{\partial x_{j}}] + P_{ε} + D_{ε} + ρ C_{1} S ε

Production Modeling

Turbulent Kinetic Energy k

P_{k} = μ_{t} S^{2}

Turbulent Dissipation Rate ε

P_{ε} = C_{ε 1} \frac{ε}{k} P_{k}

Dissipation Modeling

Turbulent Kinetic Energy k

D_{k} = - ρ ε

Turbulent Dissipation Rate ε

D_{ε} = - ρ C_{ε 2} \frac{ε^{2}}{k + \sqrt{ν ε}}

Modeling of Turbulent Viscosity $μ_{t}$

μ_{t} = C_{μ} \frac{k^{2}}{ε}

Model Coefficients

$C_{ε 1}$ = 1.44, $C_{ε 2}$ = 1.9, $σ_{k}$ = 1.0, $σ_{ε}$ = 1.22.

where $U^{*}$ = $\sqrt{S_{i j} S_{i j} + \tilde{Ω_{i j}} \tilde{Ω_{i j}}}$ , $\tilde{Ω_{i j}} = Ω_{i j} - 2 ε_{i j k} ω_{k}$ , $Ω_{i j} = \bar{Ω_{i j}} - 2 ε_{i j k} ω_{k}$ , $A_{0} = 4.04$ , $A_{s} = \sqrt{6} c o s ϕ$ , $ϕ = \frac{1}{3} c o s^{- 1} (\sqrt{6} W)$ , $W = \frac{S_{i j} S_{j k} S_{k j}}{\sqrt{S_{i j} S_{i j}}}$ ,

$C_{1} = m a x [0.43, \frac{η}{η + 5}]$

where $η = S \frac{k}{ε}$ , $S = \sqrt{2 S_{i j} S_{i j}}$ is the strain rate magnitude.

Wilcox k-ω Model

Since all three k-ε turbulence models cannot be integrated all the way to walls, wall damping wall functions must be employed to provide correct near wall behavior. It is also known that the standard k-ε turbulence model fails to predict the flow separation under adverse pressure gradients.

Wilcox proposed a turbulence model similar to the standard k-ε turbulence model but replaced the dissipation rate (ε) equation with the eddy frequency (ω) equation (Wilcox, 2006; Wilcox, 2008). The eddy frequency (ω) is often referred to the specific dissipation rate and is defined as $ω = ε / k$ . The Wilcox k-ω turbulence model has an advantage over the k-ε turbulence model as the k-ω model does not require any wall functions for the calculation of the velocity distribution near walls. As a result, the k-ω turbulence model has better performance for flows with adverse pressure gradient when compared to the k-ε turbulence models. However, the k-ω model exhibits a strong sensitivity to the freestream boundary condition (Wilcox, 2006) for external flow applications.

Transport Equations

Turbulent Kinetic Energy k

\frac{\partial (ρ k)}{\partial t} + \frac{\partial (ρ \bar{u_{j}} k)}{\partial x_{j}} = \frac{\partial}{\partial x_{j}} [(μ + σ_{k} μ_{t}) \frac{\partial k}{\partial x_{j}}] + P_{k} + D_{k}

Eddy Frequency (Specific Dissipation Rate) ω

\frac{\partial (ρ k)}{\partial t} + \frac{\partial (ρ \bar{u_{j}} k)}{\partial x_{j}} = \frac{\partial}{\partial x_{j}} [(μ + σ_{k} μ_{t}) \frac{\partial k}{\partial x_{j}}] + P_{k} + D_{k}

Production Modeling

Turbulent Kinetic Energy k

P_{k} = μ_{t} S^{2}

Eddy Frequency ω

P_{ω} = \frac{γ ω}{k} μ_{t} S^{2}

where $γ = \frac{β_{0}}{β^{*}} - \frac{σ_{ω} κ^{2}}{\sqrt{β^{*}}}$ , $β = β_{0} f_{β}$ , $f_{β} = \frac{1 + 85 χ_{ω}}{1 + 100 χ_{ω}}$ , $χ_{ω} = | \frac{Ω_{i j} Ω_{j k} {\hat{S}}_{k i}}{{(β^{*} ω)}^{3}} |$ , ${\hat{S}}_{k i} = S_{k i} - \frac{1}{2} \frac{\partial \bar{u_{m}}}{\partial x_{m}} δ_{k i}$ , $S_{i j} = \frac{1}{2} (\frac{\partial \bar{u_{i}}}{\partial x_{j}} + \frac{\partial \bar{u_{j}}}{\partial x_{i}})$ , $Ω_{i j} = \frac{1}{2} (\frac{\partial \bar{u_{i}}}{\partial x_{j}} - \frac{\partial \bar{u_{j}}}{\partial x_{i}})$

Dissipation Modeling

Turbulent Kinetic Energy (k)

D_{k} = - ρ β^{*} k ω

Eddy Frequency (ω)

D_{ω} = - ρ β ω^{2}

Modeling of Turbulent Viscosity $μ_{t}$

μ_{t} = \frac{k}{\overset{´}{ω}}

where $\overset{´}{ω} = m a x [ω, C_{l i m} \sqrt{\frac{2 {\bar{S}}_{i j} {\bar{S}}_{i j}}{β^{*}}}]$ , ${\bar{S}}_{i j} = S_{i j} - \frac{1}{3} \frac{\partial \bar{u_{k}}}{\partial x_{k}} δ_{i j}$ , $C_{l i m} = \frac{7}{8}$ ,

Model Coefficients

$σ_{k}$ = 0.6, $σ_{ω}$ = 0.5, $β^{*}$ = 0.09, $β_{0}$ = 0.0708, $κ$ = 0.4, $σ_{d} = {\begin{matrix} 0.0 f o r \frac{\partial k}{\partial x_{j}} \frac{\partial ω}{\partial x_{j}} \leq 0 \\ \frac{1}{8} f o r \frac{\partial k}{\partial x_{j}} \frac{\partial ω}{\partial x_{j}} > 0 \end{matrix}$ .

Menter Shear Stress Transport (SST) k-ω Model

Menter (1994) suggested the SST model to overcome the freestream value sensitivity of the standard k-ω turbulence model by transforming the k-ε model into the k-ω model in the near-wall region, and by utilizing the k-ε model in the turbulent region far from the wall.

The SST model employs a modified source term in the eddy frequency equation. Among the different SST versions, the SST 2003 model (Menter, 2003) will be focused on in this section as it has been shown to be more accurate when predicting flows near stagnation zones by introducing a production limiter to constrain the kinetic energy production. This version also employs the strain rate magnitude in the definition of eddy viscosity rather than the vorticity magnitude employed in the standard SST model (1994). When compared to the k-ε model, the SST 2003 model achieves better accuracy for attached boundary layers and flow separation. The SST 2003 model also overcomes the freestream sensitivity of the k-ω turbulence model. However, this model shares a similar range of weakness with the k-ε equation models for the predictions of flows with strong streamline curvature and/or rotation, as well as flows under extra strains and body forces.

Transport Equations

Turbulent Kinetic Energy k

\frac{\partial (ρ k)}{\partial t} + \frac{\partial (ρ \bar{u_{j}} k)}{\partial x_{j}} = \frac{\partial}{\partial x_{j}} [(μ + σ_{k} μ_{t}) \frac{\partial k}{\partial x_{j}}] + {\tilde{P}}_{k} + D_{k}

Eddy Frequency (Specific Dissipation Rate) ω

\frac{\partial (ρ ω)}{\partial t} + \frac{\partial (ρ \bar{u_{j}} ω)}{\partial x_{j}} = \frac{\partial}{\partial x_{j}} [(μ + σ_{ω} μ_{t}) \frac{\partial ω}{\partial x_{j}}] + P_{ω} + D_{ω} + 2 (1 - F_{1}) \frac{ρ σ_{ω 2}}{ω} \frac{\partial k}{\partial x_{j}} \frac{\partial ω}{\partial x_{j}}

where

The blending function $F_{1} = t a n h {{(m i n [m a x \frac{\sqrt{k}}{β^{*} ω d}, \frac{500 υ}{d^{2} ω}, \frac{4 ρ σ_{ω 2} k}{C D_{k ω} d^{2}}])}^{4}}$ ,
$F_{1} {\begin{matrix} 0 : a c t i v a t i o n o f k - ε m o d e l f o r t u r b u l e n t c o r e f l o w s \\ 1 : a c t i v a i t i o n o f k - ω f o r f l o w s n e a r w a l l s \end{matrix}$
$C D_{k ω} = m a x (2 ρ σ_{ω 2} \frac{1}{ω} \frac{\partial k}{\partial x_{i}} \frac{\partial ω}{\partial x_{i}}, 10^{- 10})$ ,
$d$ is the distance to the nearest wall.

Production Modeling

Turbulent Kinetic Energy k

{\tilde{P}}_{k} = \min (P_{k}, 10 ρ β^{*} k ω)

where

Production $P_{k} = μ_{t} S^{2}$
Strain rate magnitude $S = \sqrt{2 S_{i j} S_{i j}}$
Strain rate tensor $S_{i j} = \frac{1}{2} (\frac{\partial \bar{u_{i}}}{\partial x_{j}} + \frac{\partial \bar{u_{j}}}{\partial x_{i}})$
Constant $β^{*}$

Eddy Frequency (ω)

P_{ω} = α ρ S^{2}

where

Blending function for specifying constants $α = α_{1} F_{1} + α_{2} (1 - F_{1})$
$α_{1} = \frac{5}{9}$ for the inner layer
$α_{2}$ = 0.44 for the outer layer.

Dissipation Modeling

Turbulent Kinetic Energy k

D_{k} = - ρ β^{*} k ω

Eddy Frequency ω

D_{ω} = - ρ β ω^{2}

Modeling of Turbulent Viscosity $μ_{t}$

μ_{t} = \frac{ρ a_{1} k}{max (a_{1} ω, S F_{2})}

where

The second blending function $F_{2} = t a n h {(m a x (2 \frac{\sqrt{k}}{β^{*} ω d}, \frac{500 ν}{d^{2} ω}))}^{2}$ ,
Fluid kinematic viscosity $ν$ ,
Constant $β^{*}$ .

Model Coefficients

$β^{*}$ = 0.09.

Following constants for SST are computed by a blending function $ϕ = ϕ_{1} F_{1} + ϕ_{2} (1 - F_{1})$ : $σ_{ω 1} = 0.5$ , $σ_{ω 2} = 0.856$ , $σ_{k 1} = 0.85$ , $σ_{k 2} = 1.00$ , $β_{1} = \frac{3}{40}$ , $β_{2} = 0.0828$ .

Shear Stress Transport (SST) Model with Rotation/Curvature Correction

Since the SST model relies on the Boussineq approximation, it also has poor performance for the prediction of flows with streamline curvature and system rotation because the Reynolds stress tensor is aligned to the mean strain rate tensor.

In order to remedy this problem, the rotation curvature correction for the SST 2003 model modifies the production terms in both the kinetic energy equation and eddy frequency equation (Smirnov and Menter, 2009). The corrected model yields better performance over the SST base model when predicting flows with strong streamline curvature and rotating flows caused by strong swirl or geometric constraints.

Production Modeling

Turbulent Kinetic Energy k

{\tilde{P}}_{k} = \min (f_{r 1} P_{k}, 10 f_{r 1} ρ β^{*} k ω)

where

$f_{r 1}$ = $m a x (m i n [f_{r o t a t i o n}, 1.25], 0.0)$
$f_{r o t a t i o n}$ = $(1 + C_{r 1}) \frac{2 r^{*}}{1 + r^{*}} (1 - C_{r 3} t a n^{- 1} [C_{r 2} \hat{r}]) - C_{r 1}$
$\hat{r} = \frac{2 Ω_{i j} S_{i j}}{\sqrt{Ω_{i j} Ω_{i j}} D} (\frac{D S_{i j}}{D t} + [ε_{i m n} S_{j n} + ε_{j m n} S_{i n}] Ω_{m}^{r o t})$
$Ω_{i j} = \frac{1}{2} [(\frac{\partial \bar{u_{i}}}{\partial x_{j}} - \frac{\partial \bar{u_{j}}}{\partial x_{i}}) + 2 ε_{m j i} Ω_{m}^{r o t}]$
$S_{i j} = \frac{1}{2} (\frac{\partial \bar{u_{i}}}{\partial x_{j}} + \frac{\partial \bar{u_{j}}}{\partial x_{i}})$
$D = \sqrt{m a x (S^{2}, 0.09 ω^{2})}$
$r^{*} = \frac{S}{W}$
$S = \sqrt{2 S_{i j} S_{i j}}$
$W = \sqrt{2 W_{i j} W_{i j}}$

Eddy Frequency ω

P_{ω} = f_{r 1} α ρ S^{2}

Model Coefficients

$β^{*}$ = 0.09, $C_{r 1}$ = 1.0, $C_{r 2}$ = 2.0, $C_{r 3}$ = 1.0.

Turbulent Transition Models

All of the previously described models are incapable of predicting boundary layer transition. To include the effects of transition additional equations are necessary.

Transitional flows can be found in many industrial application cases including gas turbine blades, airplane wings and wind turbines. It is known that conventional turbulence models over predict the wall shear stress for transitional flows. Thus, transition models can be used to improve the accuracy of CFD solutions when flows encounter turbulent transition of the boundary layer.

γ-Reθ Model

Langtry and Menter (2009) proposed one of the most commonly used transition models in industrial CFD applications. The $γ - R e_{θ}$ model is a correlation based intermittency model that predicts natural, bypass and separation induced transition mechanisms.

The $γ - R e_{θ}$ model is coupled with the Shear Stress Transport (SST) turbulence model and requires two additional transport equations for the turbulence intermittency ( $γ$ ) and transition momentum thickness Reynolds number ( $R e_{θ}$ ). The turbulence intermittency is a measure of the flow regime and is defined as $γ = \frac{t_{turbulent}}{t_{laminar} + t_{turbulent}}$ where $t_{turbulent}$ represents the time that the flow is turbulent at a given location, while $t_{laminar}$ represents the time that the flow is laminar.

For example, while the intermittency value is zero, the flow is considered to be laminar. A value of one is for fully turbulent flow. The transition momentum thickness Reynolds number is responsible for capturing the non local effects of turbulence intensity and is defined as $R e_{θ} = \frac{ρ \bar{u} θ}{μ}$ where the momentum thickness is $θ = \int_{0}^{δ} \frac{u}{\bar{u}} (1 - \frac{u}{\bar{u}}) d y$ .

To couple with the SST model, Langtry and Menter (2009) modified the production term and dissipation term of the turbulent kinetic energy to account for the changes in the flow intermittency.

Transport Equations

Turbulent Kinetic Energy k

\frac{\partial (ρ k)}{\partial t} + \frac{\partial (ρ \bar{u_{j}} k)}{\partial x_{j}} = \frac{\partial}{\partial x_{j}} [(μ + σ_{k} μ_{t}) \frac{\partial k}{\partial x_{j}}] + {\tilde{P}}_{k} + D_{k}

Eddy Frequency (Specific Dissipation Rate) ω

\frac{\partial (ρ ω)}{\partial t} + \frac{\partial (ρ \bar{u_{j}} ω)}{\partial x_{j}} = \frac{\partial}{\partial x_{j}} [(μ + σ_{ω} μ_{t}) \frac{\partial ω}{\partial x_{j}}] + P_{ω} + D_{ω} + 2 (1 - F_{1}) \frac{ρ σ_{ω 2}}{ω} \frac{\partial k}{\partial x_{j}} \frac{\partial ω}{\partial x_{j}}

Turbulent Intermittency $γ$

\frac{\partial (ρ γ)}{\partial t} + \frac{\partial (ρ \bar{u_{j}} γ)}{\partial x_{j}} = \frac{\partial}{\partial x_{j}} [(μ + \frac{μ_{t}}{σ_{f}}) \frac{\partial γ}{\partial x_{j}}] + P_{γ} + D_{γ}

Transition Momentum Thickness Reynolds Number ${Re}_{θ}$

\frac{\partial (ρ {Re}_{θ})}{\partial t} + \frac{\partial (ρ \bar{u_{j}} {Re}_{θ})}{\partial x_{j}} = \frac{\partial}{\partial x_{j}} [σ_{f} (μ + μ_{t}) \frac{\partial {Re}_{θ}}{\partial x_{j}}] + P_{θ}

Production Modeling

Turbulent Kinetic Energy k

{\tilde{P}}_{k} = γ_{e f f} min (P_{k}, 10 ρ β^{*} k ω) \frac{\partial \bar{u_{i}}}{\partial x_{j}}

where

Production $P_{k} = μ_{t} S^{2}$
Strain rate magnitude $S = \sqrt{2 S_{i j} S_{i j}}$
Strain rate tensor $S_{i j} = \frac{1}{2} (\frac{\partial \bar{u_{i}}}{\partial x_{j}} + \frac{\partial \bar{u_{j}}}{\partial x_{i}})$
Constant $β^{*}$

Turbulent Dissipation Rate ε

P_{ε} = α ρ S^{2}

where

Blending function for specifying constants $α = α_{1} F_{1} + α_{2} (1 - F_{1})$
$α_{1} = \frac{5}{9}$ for the inner layer
$α_{2}$ = 0.44 for the outer layer.
Turbulent Intermittency $γ$

P_{γ} = F_{l e n g t h} C_{a 1} S \sqrt{(γ F_{o n s e t})} (1 - γ C C_{e 1})

where $F_{o n s e t} = m a x {m i n (m a x [\frac{R e_{ν}}{2.193 R e_{θ 0}}, {(\frac{R e_{ν}}{2.193 R e_{θ 0}})}^{4}], 2) - m a x (1 - {(\frac{R_{T}}{2.5})}^{3}, 0)}$ $R e_{ν} = \frac{S d^{2}}{ν}$ , $R_{T} = \frac{k ν}{ω}$

Transition Momentum Thickness Reynolds Number ${Re}_{θ}$

P_{θ} = C_{θ t} \frac{1}{t} ({Re}_{θ} - \bar{{Re}_{θ}}) (1 - F_{θ})

where $F_{θ} = \min (\max [F_{w a k e} e^{- {(\frac{d}{δ})}^{4}}, 1 - {(\frac{γ - 1 / C_{e 2}}{1 - 1 / C_{e 2}})}^{2}], 1)$ , $F_{w a k e} = e^{- {(\frac{R e_{ω}}{1 E + 5})}^{2}}$ , $R e_{ω} = \frac{ω d^{2}}{ν}$ , $δ = \frac{50 Ω d}{\bar{u}} δ_{B L}$ , $δ_{B L} = 7.5 θ_{B L}$ , $θ_{B L} = \frac{\bar{R e_{θ}} ν}{\bar{u}}$

Dissipation Modeling

Turbulent Kinetic Energy k

D_{k} = - \min (\max [γ_{e f f}, 0.1], 1.0) ρ β^{*} k ω

Turbulent Dissipation Rate ω

D_{ω} = - ρ β ω^{2}

Turbulent Intermittency $γ$

D_{γ} = - C_{a 2} F_{turbulent} Ω γ (C_{e 2} γ - 1)

where $F_{turbulent} = e^{{(- \frac{R_{T}}{4})}^{4}}$

Modeling of Turbulent Viscosity $μ_{t}$

μ_{t} = \frac{ρ a_{1} k}{max (a_{1} ω, S F_{2})}

where

The second blending function $F_{2} = t a n h {(m a x (2 \frac{\sqrt{k}}{β^{*} ω d}, \frac{500 ν}{d^{2} ω}))}^{2}$ ,
Fluid kinematic viscosity $ν$ ,
Constant $β^{*}$ .

Correlations

R e_{θ} = (1173.51 - 589.428 T u + \frac{0.2196}{T u^{2}}) F (λ_{θ})

R e_{θ t} = {\begin{matrix} (1173.51 - 589.428 T u + \frac{0.2196}{T u^{2}}) F (λ_{θ}) f o r T u \leq 1.3 \\ 331.5 {(T u - 0.5658)}^{- 0.671} F (λ_{θ}), f o r T u > 1.3 \end{matrix}

F (λ_{θ}) {\begin{matrix} 1 - (- 12.986 λ_{θ} - 123.66 λ_{θ}^{2} - 405.689 λ_{θ}^{3}) e^{- {(\frac{T u}{1.5})}^{5}} f o r λ_{θ} \leq 0 \\ 1 + 0.275 {(1 - e^{- 35 λ_{θ}})}^{- 0.671} e^{- (\frac{T u}{1.5})}, f o r λ_{θ} > 0 \end{matrix}

where $λ_{θ} = \frac{θ^{2}}{ν} \frac{d \bar{u}}{d s}$ is the pressure gradient.

Model Coefficients

$β^{*}$ = 0.09.

The following constants for SST are computed by a blending function $ϕ = ϕ_{1} F_{1} + ϕ_{2} (1 - F_{1})$ : $σ_{ω 1}$ = 0.5, $σ_{ω 2}$ = 0.856, $σ_{k 1}$ = 0.85, $σ_{k 2}$ = 1.00, $β_{1} = \frac{3}{40}$ , $β_{2}$ = 0.0828. $C_{e 1}$ = 1.0, $C_{e 2}$ = 50, $C_{a 1}$ = 2.0, $C_{a 2}$ = 0.06, $C_{θ 1}$ = 1.0, $C_{θ 2}$ = 2.0, $σ_{f}$ = 1.0.

Large Eddy Simulation (LES)

LES is often regarded as an impractical tool for industrial CFD applications as it requires large computational resources.

LES is associated with the severe requirements of small mesh size and time step size so that the turbulent length scale and time scale in the inertial range is properly resolved. These simulations also have an issue with resolving flows close to walls where energetic vortical structures become very small with the increase of Reynolds number. However, there are cases where LES provides significant benefits, including combustion chemical reaction problems and acoustic problems.

Filtered Navier-Stokes Equations

While Reynolds-averaged Navier-Stokes (RANS) resolves the mean flow and requires a turbulence model to account for the effect of turbulence on the mean flow, Large Eddy Simulation (LES) computes both the mean flow and the large energy containing eddies.

A subgrid model is used to capture the effects of small scale turbulent structures. The spatial filtering process is used to filter out the turbulent structures from the instantaneous flow field which are smaller than a given filter size. This filtering process is based on the decomposition of instantaneous variables (velocity, pressure) into filtered (resolved) and sub filter (unresolved or residual) variables. Here, velocity is used as an example.

u_{i} = \tilde{u_{i}} + u "_{i}

where

$u_{i}$ : instantaneous velocity,
$\tilde{u_{i}}$ : filtered velocity,
$u "_{i}$ : sub filtered (unresolved) velocity.

The filtered velocity field is obtained from a low pass filtering operation due to a weighted filter G, defined as $\tilde{u_{i}} = \int \int \int G (x, x^{'}, Δ) u_{i} (x^{'}, t) d x'_{1} d x'_{2} d x'_{3}$ .

The weighted filter includes

$G (x, x^{'}) = \frac{1}{π} (\frac{sin {(\frac{π (x - x^{'})}{Δ})}^{2}}{(x - x')})$ : a cut-off filter,
$G (x, x^{'}) = \sqrt{\frac{6}{π Δ^{2}}} e x p (\frac{- 6 {(x - x')}^{2}}{Δ^{2}})$ : a Gaussian filter,
$G (x, x^{'}) = {\begin{matrix} \frac{1}{Δ} \forall | x - x^{'} | \leq \frac{1}{2} Δ \\ 0 \forall | x - x^{'} | > \frac{1}{2} Δ \end{matrix}$ : a top-hat filter,

where $Δ = \sqrt[3]{Δ x Δ y Δ z}$ is a cutoff width, representing a spatial averaging over a grid element. Among those, a top-hat filter (or similar one) is a popular choice for commercial CFD codes where unstructured meshes are usually adopted.

Although the LES decomposition method resembles the Reynolds method, they have an important difference due to the following. $\tilde{(\tilde{u_{i}})} \neq \tilde{u_{i}}$ , $\tilde{u "_{i}} \neq 0$

Once this concept is substituted into the instantaneous Navier-Stokes equations, and then the spatial-averaging (or filtering) is made, the filtered Navier-Stokes equations are obtained. The filtered Navier-Stokes equations include the equations for the filtered continuity and filtered momentum equations, which are given below.

\frac{\partial (ρ \tilde{u_{i}})}{\partial x_{i}} = 0

\frac{\partial (ρ \tilde{u_{i}})}{\partial t} + \frac{\partial (ρ \tilde{u_{i} u_{j}})}{\partial x_{j}} = - \frac{\partial (\tilde{p} δ_{i j})}{\partial x_{j}} + 2 μ \frac{\partial \tilde{S_{i j}}}{\partial x_{j}}

where

$\tilde{u_{i} u_{j}} = \tilde{(\tilde{u_{i}} \tilde{u_{j}})} + \tilde{(u "_{i} \tilde{u_{j}}) +} \tilde{(\tilde{u_{i}} u "_{j})} + \tilde{(u "_{i} u "_{j})}$ is due to the decomposition of the nonlinear convective term in the momentum equation.
$\tilde{S_{i j}} = \frac{1}{2} (\frac{\partial \tilde{u_{i}}}{\partial x_{j}} + \frac{\partial \tilde{u_{j}}}{\partial x_{i}})$ is the filtered strain rate tensor.

The filtered momentum equation is rearranged as

\frac{\partial (ρ \tilde{u_{i}})}{\partial t} + \frac{\partial (ρ \tilde{\tilde{u_{i}} \tilde{u_{j}}})}{\partial x_{j}} = - \frac{\partial \tilde{p}}{\partial x_{j}} + 2 μ \frac{\partial \tilde{S_{i j}}}{\partial x_{j}} + \frac{\partial τ_{i j}^{*}}{\partial x_{j}}

where

$τ_{i j}^{*} = - ρ (C_{i j} + R_{i j})$ is a double decomposition stress tensor.
$C_{i j} = \tilde{(u "_{i} \tilde{u_{j}}) +} \tilde{(\tilde{u_{i}} u "_{j})}$ is the cross stress tensor, representing the interactions between large and small turbulence eddies.
$R_{i j} = \tilde{(u "_{i} u "_{j})}$ is the Reynolds subgrid stress tensor.

Since $\tilde{(\tilde{u_{i}} \tilde{u_{j}})}$ in the convection term of the filtered momentum equation needs a secondary filtering process, it is rewritten as shown below (Leonard, 1974),

\tilde{(\tilde{u_{i}} \tilde{u_{j}})} = L_{i j} + \tilde{u_{i}} \tilde{u_{j}}

where $L_{i j} = \tilde{(\tilde{u_{i}} \tilde{u_{j}})} - \tilde{u_{i}} \tilde{u_{j}}$ is the Leonard stress tensor, representing interactions among large turbulent eddies.

Utilizing the decomposition process shown above, the filtered momentum equations can be rewritten as

\frac{\partial (ρ \tilde{u_{i}})}{\partial t} + \frac{\partial (ρ \tilde{u_{i}} \tilde{u_{j}})}{\partial x_{j}} = - \frac{\partial \tilde{p}}{\partial x_{j}} + 2 μ \frac{\partial \tilde{S_{i j}}}{\partial x_{j}} + \frac{\partial τ_{i j}^{'}}{\partial x_{j}}

where

$\tilde{S_{i j}} = \frac{1}{2} (\frac{\partial \tilde{u_{i}}}{\partial x_{j}} + \frac{\partial \tilde{u_{j}}}{\partial x_{i}})$ is the filtered strain rate tensor,
$τ_{i j}^{'} = - ρ (C_{i j} + R_{i j} + L_{i j}) = ρ \tilde{u_{i} u_{j}} - ρ \tilde{u_{i}} \tilde{u_{j}}$ is the subgrid stress tensor.

Similar to RANS, the subgrid stress tensor is an unknown and must be computed by a subgrid model. However, the subgrid stress tensor differs from the Reynolds stress tensor. The table below summarizes the differences between time averaging for RANS and spatial averaging for LES. Common subgrid scale (SGS) models available to solve the filtered Navier-Stokes equations include:

Smagorinsky-Lilly SGS model
Germano Dynamic Smagorinsky-Lilly model
Wall-Adapting Local Eddy-Viscosity (WALE)

Table 1. Time-Averaging (RANS) vs. Spatial-Averaging (LES)
	RANS	LES
Double averaging	$\bar{(\bar{u_{i}})} = \bar{u_{i}}$	$\tilde{(\tilde{u_{i}})} \neq \tilde{u_{i}}$
Turbulence averaging	$\bar{u'_{i}} = 0$	$\tilde{u "_{i}} \neq 0$
Averaging (or filtering) of convective term of Navier Stokes	$\bar{u_{i} u_{j}} = (\bar{\bar{u_{i}} \bar{u_{j}}}) + (\bar{u'_{i} \bar{u_{j}}}) +$ $(\bar{\bar{u_{i}} u'_{j}}) + \bar{u'_{i} u'_{j}}$	$\tilde{u_{i} u_{j}} = \tilde{(\tilde{u_{i}} \tilde{u_{j}})} + \tilde{(u "_{i} \tilde{u_{j}}) +}$ $\tilde{(\tilde{u_{i}} u "_{j})} + \tilde{(u "_{i} u "_{j})}$
Turbulent stress	$τ_{i j}^{R} = - ρ \bar{u_{i}^{'} u_{j}^{'}}$	$τ_{i j}^{'} = ρ \tilde{u_{i} u_{j}} - ρ \tilde{u_{i}} \tilde{u_{j}}$

Smagorinsky-Lilly Subgrid Scale Model

The Smagorinsky-Lilly SGS model is based on the Prandtl’s mixing length model and assumes that a kinematic SGS viscosity can be expressed in terms of the length scale and the strain rate magnitude of the resolved flow.

τ_{i j}^{'} = ρ \tilde{u_{i} u_{j}} - ρ \tilde{u_{i}} \tilde{u_{j}} = - 2 μ_{s} \tilde{S_{i j}}

where

$μ_{s} = ρ {(C_{S} Δ)}^{2} \sqrt{2 \tilde{S_{i j}} \tilde{S_{i j}}}$ is the subgrid turbulent viscosity,
Model constant 0.17 < $C_{s}$ < 0.21 (Lilly, 1996)
$\tilde{S_{i j}} = \frac{1}{2} (\frac{\partial \tilde{u_{i}}}{\partial x_{j}} + \frac{\partial \tilde{u_{j}}}{\partial x_{i}})$ is the filtered (resolved) strain rate tensor

A major drawback of this model is that the model constant ( $C_{s}$ ) does not vary in space and time. Furthermore, this model has no correct wall behavior and it is too dissipative for laminar turbulent transition cases. These limitations led to the development of a dynamic model for which the model constant is allowed to vary depending on the grid resolution and flow regime.

Dynamic Subgrid Scale Model

Recognizing $C_{s}$ variations in space and time, Germano et al. (1991) proposed the dynamic model to compute the value of $C_{s}$ rather than specifying it explicitly.

It is implemented by utilizing two filters: a cutoff filter $Δ$ and a test (coarse) cutoff filter $\tilde{Δ}$ .

The subgrid stress tensor $τ_{i j}^{'}$ with the cutoff filter ( $Δ$ ) is: $τ_{i j}^{'} = ρ \tilde{u_{i} u_{j}} - ρ \tilde{u_{i}} \tilde{u_{j}} = - 2 ρ {(C_{S} Δ)}^{2} | \tilde{S} | \tilde{S_{i j}}$ . Where $| \tilde{S} | = \sqrt{2 \tilde{S_{i j}} \tilde{S_{i j}}}$ is the strain rate magnitude.

Figure 3 shows a resolved turbulence region utilizing Large Eddy Simulation (LES) and a modeled region assuming the subgrid tensor

τ_{i j}^{'}

Figure 3. Energy Spectrum for LES Using the Cutoff Filter Width ( $Δ$ )

The test subgrid stress tensor $T_{i j}$ with the coarse filter ( $\tilde{Δ}$ ) can be written as

T_{i j} = ρ \tilde{\tilde{u_{i} u_{j}}} - ρ \tilde{\tilde{u_{i}}} \tilde{\tilde{u_{j}}} = - 2 ρ {(C_{S} \tilde{Δ})}^{2} | \tilde{\tilde{S}} | \tilde{\tilde{S_{i j}}}

where

$| \tilde{\tilde{S}} | = \sqrt{2 \tilde{\tilde{S_{i j}}} \tilde{\tilde{S_{i j}}}}$ is the coarse filtered strain rate magnitude.
$\hat{\tilde{S_{i j}}} = \frac{1}{2} (\frac{\partial \hat{\tilde{u_{i}}}}{\partial x_{j}} + \frac{\partial \hat{\tilde{u_{j}}}}{\partial x_{i}})$ is the filtered strain rate tensor, using the coarse cutoff filter.

Figure 4 shows a resolved LES region and a corresponding subgrid modelled region (T) when the coarse filter is employed.

Figure 4. Energy Spectrum for LES Using the Test Cutoff Filter Width ( $\tilde{Δ}$ )

Because of the coarse filtering, the test (coarse) subgrid stress tensor $T_{i j}$ should be a summation of the coarse filtered subgrid stress tensor $\tilde{τ_{i j}^{'}}$ and the Leonard stress tensor $L_{i j}$ .

L_{i j} = T_{i j} - \tilde{τ_{i j}^{'}} = ρ \tilde{\tilde{u_{i} u_{j}}} - ρ \tilde{\tilde{u_{i}}} \tilde{\tilde{u_{j}}} - ρ \tilde{\tilde{u_{i} u_{j}}} + ρ \tilde{\tilde{u_{i}} \tilde{u_{j}}} = ρ \tilde{\tilde{u_{i}} \tilde{u_{j}}} - ρ \tilde{\tilde{u_{i}}} \tilde{\tilde{u_{j}}}

where

$\tilde{τ_{i j}^{'}}$ is the subgrid tensor for the cutoff filter (or grid filtered), then test filtered.

\tilde{τ_{i j}^{'}} = ρ \tilde{\tilde{u_{i} u_{j}}} - ρ \tilde{\tilde{u_{i}} \tilde{u_{j}}} = - 2 ρ {(C_{S} Δ)}^{2} \tilde{| \tilde{S} | \tilde{S_{i j}}}

$L_{i j}$ is the Leonard subgrid stress tensor, representing the contribution to the subgrid stresses by turbulence length scales smaller than the test filter but larger than the cutoff filter.

The Leonard subgrid stress tensor can be arranged as

L_{i j} = - 2 ρ {(C_{S} \hat{Δ})}^{2} | \tilde{\tilde{S}} | \tilde{\tilde{S_{i j}}} + 2 ρ {(C_{S} Δ)}^{2} \tilde{| \tilde{S} | \tilde{S_{i j}}} = 2 ρ C_{S}^{2} Δ^{2} (\tilde{| \tilde{S} | \tilde{S_{i j}}} - α^{2} | \tilde{\tilde{S}} | \tilde{\tilde{S_{i j}}})

where $α = \hat{Δ} / Δ$ .

The Leonard subgrid stress tensor can be rewritten as

L_{i j} = 2 ρ C_{S}^{2} Δ^{2} (\tilde{| \tilde{S} | \tilde{S_{i j}}} - α^{2} | \tilde{\tilde{S}} | \tilde{\tilde{S_{i j}}}) = C_{S}^{2} M_{i j}

where $M_{i j} = 2 ρ Δ^{2} (\tilde{| \tilde{S} | \tilde{S_{i j}}} - α^{2} | \tilde{\tilde{S}} | \tilde{\tilde{S_{i j}}})$ .

Since the above equation is overdetermined a minimum least square error method is used to determine the coefficient $C_{s}$ .

C_{s}^{2} = \frac{L_{i j} M_{i j}}{M_{i j} M_{i j}}

In order to avoid numerical instabilities associated with the above equation, as the numerator could become negative, averaging of the error in the minimization is employed.

C_{s}^{2} = \frac{L_{i j} M_{i j}}{M_{i j} M_{i j}}

Hybrid Simulations

In recent years, hybrid methods have increasingly been employed for the simulation of unsteady turbulent flows.

These models form a bridge between Large Eddy Simulation (LES) and Reynolds-averaged Navier-Stokes (RANS) by utilizing RANS for attached boundary layer calculations and LES for the separated regions. In general, hybrid simulations need spatial filtering processes to determine the local sub grid turbulent viscosity. Compared to LES and Direct Numerical Simulation (DNS), hybrid simulations are extremely robust since the numerical requirement is less severe than the two other approaches. Hybrid simulations include Detached Eddy Simulation (DES), Delayed Detached Eddy Simulation (DDES) and Improved Delayed Detached Eddy Simulation (IDDES).

Detached Eddy Simulations

Detached Eddy Simulation (DES) uses one of the one or two equation Reynolds-averaged Navier-Stokes (RANS) turbulence models to define the turbulence length scale or distance from the wall. The distance is then used to determine the region for which the RANS or Large Eddy Simulation (LES) models will be used.

DES using the Spalart-Allmaras (SA) model employs a modified distance function to replace the distance (d) in the SA model’s dissipation term (Spalart et al., 1997).

\tilde{d} = \min (d, C_{D E S} Δ)

where d is the distance to the closest wall, $C_{D E S}$ is a constant and $Δ$ is a metric of the local element size, often chosen as the largest grid spacing in all three directions ( $Δ = max [Δ x, Δ y, Δ z]$ ).

The SA model is activated in a boundary layer region where d < ∆, while a Smagorinsky-like LES model is used in regions where ∆ < d. Similarly, Strelets (2001) adopted the SST model as the baseline model for DES, which computes a turbulence length scale and compares it with a grid length size for a switch between LES and RANS. However, both approaches need careful determinations of local grid sizes because these grid sizes act as the switch for the activation of LES and RANS.

Delayed Detached Eddy Simulations

Delayed Detached Eddy Simulation is an improved version of DES that avoids the undesired activation of Large Eddy Simulation (LES) in boundary layer regions when the maximum grid size is less than the distance to the closest wall (∆ < d).

In this model (Spalart et al., 2006), the length scale is redefined as:

\tilde{d} = d - f_{d} \max (0, d - C_{D E S} Δ)

where $f_{d} = 1 - t a n h ({(8 h)}^{3})$ , $h = \frac{\tilde{ν}}{\sqrt{{\bar{u}}_{i j} {\bar{u}}_{i j}} κ^{2} d^{2}}$ , $κ = 0.41$ .

Improved Delayed Detached Eddy Simulations (DDES)

Shur et al. (2008) proposed an improved DDES to ensure that most of the turbulence is resolved when the model is operating as a wall modeled Large Eddy Simulation (LES).

This is achieved by introducing a blending function of length scales defined as:

\tilde{d} = \tilde{f_{d}} (1 + f_{e}) d + (1 - \tilde{f_{d}}) C_{D E S} Δ

where

$\tilde{f_{d}} = \max ([1 - f_{d}], f_{b})$ ,
$f_{b} = \min (2 e^{- 9 α^{2}}, 1.0)$ ,
$α = 0.25 - \frac{d}{Δ}$ ,
$f_{e} = \max ([f_{e 1} - 1], 1.0) f_{e 2}$ ,
$f_{e 1} = {\begin{matrix} 2 e^{- 11.09 α^{2}} \frac{Δ}{d}, i f α \geq 0 \\ 2 e^{- 9 α^{2}} \frac{Δ}{d}, i f α < 0 \end{matrix}$ ,
$f_{e 2} = 1.0 - \max (\tanh {((C_{t}^{2} h_{t})}^{3}), \tanh {((C_{l}^{2} h_{l})}^{10}))$ ,
$h_{t} = \frac{ν_{t}}{\sqrt{{\bar{u}}_{i j} {\bar{u}}_{i j}} κ^{2} d^{2}}$ ,
$h_{l} = \frac{ν}{\sqrt{{\bar{u}}_{i j} {\bar{u}}_{i j}} κ^{2} d^{2}}$ .

Near-Wall Modeling

For internal wall bounded flows, proper mesh resolution is required in order to calculate the steep gradients of the velocity components, turbulent kinetic energy, dissipation, as well as the temperature.

Furthermore, the first grid locations from the wall and the stretching ratio between subsequent points are also determinant factors that affect solution accuracy. Since boundary layer thickness is reduced as the Reynolds number increases, the computational cost increases for higher Reynolds number flows due to the dense grid requirements near the walls. Additionally, it is hard to determine the first grid locations for complex three dimensional industrial problems without adequate testing and simulation.

For high Reynolds number flows, it is numerically efficient if flows within boundary layers are modeled rather than resolved down to the wall. This suggests that a coarse mesh is used at the expense of numerical accuracy when compared to fully wall resolved approaches.

Wall Function

Figure 5 shows velocity profiles over a flat plate under the zero pressure gradient.

Velocity profile and wall distances are scaled by the friction velocity (

u_{τ} = \sqrt{\frac{τ_{w}}{ρ}}

where

τ_{w}

is the wall shear stress and

ρ

is the fluid density). This is referred to as the log-law velocity profile. As shown in Figure 5, the velocity profiles are divided into three distinguished regions, the viscous sublayer, the buffer layer and the logarithmic layer.

The wall function utilizes the universality of the log-law velocity profile to obtain the wall shear stress. This wall function approach allows the reduction of mesh requirement near the wall as the first mesh location is placed outside the viscous sublayer.

Overall, the log-law based wall model is economical and reasonably accurate in most flow conditions, especially for high Reynolds number flows. However, it tends to display poor performance in situations with low Reynolds number flows, strong body forces (rotational effect, buoyancy effect), and massively separated flows with adverse pressure gradients.

Velocity Profile in the Viscous Sublayer (y+ < 5)

In the viscous sublayer, the normalized velocity profile ( $U^{+}$ ) has a linear relationship with the normalized wall distance.

U^{+} = y^{+}

where $U^{+} = \frac{\bar{u}}{u_{τ}}$ is the velocity ( $\bar{u}$ ) parallel to the wall, normalized by the friction velocity. The friction velocity is defined as $u_{τ} = \sqrt{\frac{τ_{w}}{ρ}}$ , $τ_{w}$ is the wall shear stress and $ρ$ is the fluid density. $y^{+} = \frac{ρ u_{τ} y}{μ}$ is the normalized wall distance (or wall unit). The distance from the wall is $y$ and $μ$ is the fluid dynamic viscosity.

Velocity Profile in the Logarithmic Layer (30 < y+ < 500)

In the logarithmic layer the velocity profile can be given by a logarithmic function

U^{+} = \frac{1}{κ} \log (y^{+}) + B

where $κ$ = 0.4 is the Von Kármán constant and B = 5.5 is a constant.

The wall shear stress can be estimated from the above equation via the iterative solution procedure.

Velocity Profile in the Buffer Layer (5 < y+ < 30)

Since two equations mentioned before are not valid in the buffer layer a special function is needed to bridge the viscous sublayer and the logarithmic layer. Details are not covered here.

Kinematic Eddy Viscosity

The kinematic eddy viscosity can be obtained as $ν_{t} = κ y u_{τ}$ .

Two Layer Wall Model

The two equation turbulence models based on eddy frequency (ω) and the Spalart-Allmaras (SA) model do not require special wall treatments to solve the boundary layer as these models are valid through the viscous sub layer.

However, the two equation turbulence models based on turbulence dissipation rate (ε) need additional functions to simulate the near wall effects. The two layer wall model is one of them. In the two layer model, turbulent kinetic energy is determined from the turbulent kinetic energy transport equation, while dissipation rate is resolved with a one equation turbulence model (Wolfstein, 1969) in the near wall, viscous-affected regions where $R e_{y} = \frac{ρ y \sqrt{k}}{μ} < 200$ .

ε = \frac{k^{3 / 2}}{l_{ε}}

where

$l_{ε} = C_{l} y (1 - \exp [- \frac{R e_{y}}{A_{ε}}])$ ,
$C_{l} = κ C_{μ}^{- 3 / 4}$ ,
$κ$ is the Von Kármán constant,
$A_{ε}$ = 5.08 (Chen and Patel, 1988).

The eddy viscosity can be written as

ν_{t} = C_{μ} l_{μ} \sqrt{k}

where

$l_{μ} = C_{l} y (1 - \exp [- \frac{R e_{y}}{A_{μ}}])$ ,
$A_{μ}$ = 70.

Inlet Turbulence Parameters

CFD simulations require specification of turbulence variables at inlet boundaries.

When measured turbulence data is available you can explicitly specify turbulence variables. For example, eddy viscosity for the Spalart-Allmaras (SA) model, turbulent kinetic energy and eddy frequency/dissipation rate for k-ε and k-ω based models. When measured data is not available, there are estimations for the turbulence values that are based on the turbulence intensity and turbulence length scale or eddy viscosity ratio.

Turbulence Intensity

The turbulence intensity (I) is defined as the ratio of the root-mean-square of the turbulent velocity fluctuations ( $u'$ ) and the mean velocity ( $\bar{u}$ ):

I \equiv \frac{u'}{\bar{u}}

where $u^{'} = \sqrt{\frac{(\bar{u_{x}^{'}^{2}} + \bar{u_{y}^{'}^{2}} + \bar{u_{z}^{'}^{2}})}{3}} = \sqrt{\frac{2 k}{3}}$ and $\bar{u} = \sqrt{({\bar{u}}_{x}^{2} + {\bar{u}}_{y}^{2} + {\bar{u}}_{z}^{2})}$ .

For fully developed internal flows the turbulence intensity can be estimated as

I = 0.16 R e_{h}^{- 1 / 8}

where $R e_{h} = \frac{ρ \bar{u} D_{h}}{μ}$ is the Reynolds number and $D_{h}$ is a hydraulic diameter.

The following turbulence intensity values can be assumed

5 percent < $I$ < 20 percent for rotating machineries, for example, turbines and compressors.
1 percent < $I$ < 5 percent for internal flows.
$I$ ~0.05 percent for external flows.

Turbulence Length Scale

The turbulence length scale represents the characteristic size of the turbulent eddies within a flow field. This parameter is often used to characterize the nature of the turbulence and appears in some form in nearly every turbulence model. RANS turbulence models have different definitions of the turbulence length scale based on the turbulence model and the type of application.

Estimated turbulence length scales based on the flow conditions:

$l = 0.038 D_{h}$ for fully developed internal flows.
$l = 0.22 δ$ for developing flows.

where $δ \approx 0.382 \frac{x}{R e_{x}^{1 / 5}}$ is the turbulent boundary layer thickness over a flat plate and $R e_{x} = \frac{ρ \bar{u} x}{μ}$ is the Reynolds number and x is the distance from the start of the boundary layer.

Eddy Viscosity Ratio

The eddy viscosity ratio ( $ν_{r}$ ) is the ratio between the eddy viscosity ( $ν_{t}$ ) and fluid kinematic viscosity ( $ν$ ). The eddy viscosity ratio is set between 1 and 10 for internal flows, while it should be between 0.2 and 1.3 for external flows.

Inlet Turbulence Specification for the SA Model

The SA model has three options to estimate the turbulence variable at the inlet.

Option 1: Eddy viscosity $ν_{t}$
Option 2: Turbulence intensity $I$ and length scale $l$ . For the eddy viscosity, $ν_{t} = \sqrt{\frac{3}{2}} \bar{u} I l$ .
Option 3: Viscosity ratio $ν_{r}$ . For the eddy viscosity, $ν_{t} = ν ν_{r}$

Inlet Turbulence Specification for the k-ε Models

Option 1: Kinetic energy $k$ and dissipation rate $ε$
Option 2: Turbulence intensity $I$ and length scale $l$ . For the kinetic energy, $k = \frac{3}{2} {(\bar{u} I)}^{2}$ . For the dissipation rate, $ε = C_{μ} \frac{k^{3 / 2}}{l}$ , where $C_{μ}$ = 0.09.
Option 3: Turbulence intensity $I$ and viscosity ratio $ν_{r}$ . For the kinetic energy, $k = \frac{3}{2} {(\bar{u} I)}^{2}$ . For the dissipation rate, $ε = \frac{C_{μ} k^{2}}{ν ν_{r}}$ .

Inlet Turbulence Specification for the k-ω Models

Option 1: Kinetic energy $k$ and eddy frequency $ω$
Option 2: Turbulence intensity $I$ and length scale $l$ . For the kinetic energy, $k = \frac{3}{2} {(\bar{u} I)}^{2}$ . For the eddy frequency, $ω = \frac{\sqrt{k}}{l C_{μ}^{1 / 4}}$ .
Option 3: Turbulence intensity $I$ and viscosity ratio $ν_{r}$ . For the kinetic energy, $k = \frac{3}{2} {(\bar{u} I)}^{2}$ . For the dissipation rate, $ω = \frac{k}{ν ν_{r}}$ .

Turbulence Wall Function

The first option for computing turbulent boundary layers is to fully resolve them. When computing the near wall gradients explicitly, AcuSolve integrates the governing equations directly to the wall. As a result, this option is more accurate, provided that sufficient mesh density is used. $y^{+}$ (defined below) of the first mesh point must be less than 10 (preferably 5); otherwise, gross errors in traction, heat flux, and mass transfer may result. This option is typically used for applications where the near wall flow profile plays an important role in the physics of the simulation, that is, cases having adverse pressure gradients, flow separation, and so on. This option is activated by specifying the turbulence_wall_type = low_reynolds_number (or low_re) parameter in the SIMPLE_BOUNDARY_CONDITION command.

y^{+} = \frac{ρ y u^{*}}{μ}

where $μ$ is the viscosity, $u^{*} = \sqrt{\frac{τ_{w}}{ρ}}$ is the turbulent friction velocity, $τ_{w}$ is the wall shear stress, $ρ$ is the density.

The second type of wall treatment for turbulent boundary layers allows you to approximate the near wall flow field, without using fine near wall mesh, by employing wall functions. This approach can greatly reduce the size of a mesh by eliminating the need for fine mesh spacing normal to no-slip walls. When this approach is applied, AcuSolve assumes a shape for the near wall flow field. This assumed profile is based on the “Law of the Wall” for turbulent boundary layers. The “Law of the Wall” is a relation that is based on theoretical and experimental arguments and relates the stream wise velocity profile with the normal distance from the wall. This relation was formulated for fully developed boundary layers with favorable pressure gradients. This option is activated by specifying the turbulence_wall_type = wall_function (or func) parameter in the SIMPLE_BOUNDARY_CONDITION command. $y^{+}$ of the first mesh point may be as large as 300.

The third type of wall model that is offered by AcuSolve is the running average wall function. When this model is employed, the wall function is evaluated using the running average velocity field and not the instantaneous field. This approach is typically used with LES and DES models, but may also be used with RANS if appropriate.

Note: This requires the Running Average field to be turned on in the simulation. This option is activated by specifying the turbulence_wall_type = running_average_wall_function parameter in the SIMPLE_BOUNDARY_CONDITION command.

Turbulence Modeling

Reynolds-Averaged Navier-Stokes (RANS) Simulations

General Form of Turbulence Models

One Equation Eddy Viscosity Models

Spalart-Allmaras (SA) Model

Transport Equations

Production of υ ^

Destruction of υ ^

Modeling of Turbulent Viscosity μ t

Model Coefficients

Spalart-Allmaras (SA) Model with Rotation/Curvature Correction

Modified Production of υ ^

Model Coefficients

Two Equation Eddy Viscosity Models

Standard k-ε Model

Transport Equations

Production Modeling

Dissipation Modeling

Modeling of Turbulent Viscosity μ t

Model Coefficients

Renormalization Group (RNG) k-ε Model

Transport Equations

Production Modeling

Dissipation Modeling

Modeling of Turbulent Viscosity μ t

Model Coefficients

Realizable k-ε Model

Transport Equations

Production Modeling

Dissipation Modeling

Modeling of Turbulent Viscosity μ t

Model Coefficients

Wilcox k-ω Model

Transport Equations

Production Modeling

Dissipation Modeling

Modeling of Turbulent Viscosity μ t

Model Coefficients

Menter Shear Stress Transport (SST) k-ω Model

Transport Equations

Production Modeling

Dissipation Modeling

Modeling of Turbulent Viscosity μ t

Model Coefficients

Shear Stress Transport (SST) Model with Rotation/Curvature Correction

Production Modeling

Model Coefficients

Turbulent Transition Models

γ-Reθ Model

Transport Equations

Production Modeling

Dissipation Modeling

Modeling of Turbulent Viscosity μ t

Correlations

Model Coefficients

Large Eddy Simulation (LES)

Filtered Navier-Stokes Equations

Smagorinsky-Lilly Subgrid Scale Model

Dynamic Subgrid Scale Model

Hybrid Simulations

Detached Eddy Simulations

Delayed Detached Eddy Simulations

Improved Delayed Detached Eddy Simulations (DDES)

Near-Wall Modeling

Wall Function

Velocity Profile in the Viscous Sublayer (y+ < 5)

Velocity Profile in the Logarithmic Layer (30 < y+ < 500)

Velocity Profile in the Buffer Layer (5 < y+ < 30)

Kinematic Eddy Viscosity

Two Layer Wall Model

Inlet Turbulence Parameters

Turbulence Intensity

Turbulence Length Scale

Eddy Viscosity Ratio

Inlet Turbulence Specification for the SA Model

Inlet Turbulence Specification for the k-ε Models

Inlet Turbulence Specification for the k-ω Models

Turbulence Wall Function

Production of $\hat{υ}$

Destruction of $\hat{υ}$

Modeling of Turbulent Viscosity $μ_{t}$

Modified Production of $\hat{υ}$

Modeling of Turbulent Viscosity $μ_{t}$

Modeling of Turbulent Viscosity $μ_{t}$

Modeling of Turbulent Viscosity $μ_{t}$

Modeling of Turbulent Viscosity $μ_{t}$

Modeling of Turbulent Viscosity $μ_{t}$

Modeling of Turbulent Viscosity $μ_{t}$