# 3.2.6. Turbulence Closure Models

The Favre-filtered turbulent flow equations of the previous section have been modeled in terms of , the turbulent eddy viscosity for RANS simulations and the subgrid turbulent eddy viscosity for LES. Evaluation of this eddy viscosity is dependent upon the closure model selected. All models supported by SIERRA/Fuego are described below.

## 3.2.6.1. Standard RANS Model

The standard closure model is a two-equation type of model, where transport equations for the turbulent kinetic energy and the turbulent dissipation rate are solved to obtain length-scale and time-scale estimates for the local turbulence field, to be used for modeling the turbulent eddy viscosity . The turbulent kinetic energy, , and the dissipation rate of turbulent kinetic energy, , are given by (Gran et al. [16])

(3.118)

(3.119)

respectively, where the turbulence production rate, , is defined as

(3.120)

and is modeled using the same Boussinesq approximation as in (3.85),

(3.121)

The turbulent eddy viscosity is then given by the Prandtl-Kolmogorov relationship,

(3.122)

where . The filter time, is provided by the temporally filtered Navier Stokes model (Tieszen et al. [17]). The parameters , , , and are adjustable constants.

Frequently, although not formally justified in high Reynolds flows, the diffusion coefficient for the turbulent kinetic energy and turbulence dissipation, (3.118) and (3.119), may include the molecular viscosity. This option is supported within Fuego by entering the following command line in the Fuego region block, `include molecular viscosity in k-e diffusion term`.

## 3.2.6.2. Low Reynolds Number RANS Model

In the case of the low Reynolds number turbulent flows, the standard transport equations can be modified to contain additional damping functions to improve their accuracy. The low Reynolds number model of Launder and Sharma [18] is used, along with recent modification proposed by Mathur and He [19]. In this formulation, (3.118), includes an additional right-hand-side source term,

(3.123)

Moreover, the dissipation rate equation, which is now solved for , i.e., the transformed dissipation rate that is numerically zero at the wall, also includes a new source term,

(3.124)

where the strain-rate magnitude is defined by .

The constants in the dissipation rate equation are modified by damping coefficients, and , where is unity and .

In this formulation, the eddy viscosity relationship defined by,

(3.125)

where the dampening function is and is now a function of . The constant is 3.4 with .

The physical turbulent dissipation rate is related to the numerical turbulent dissipation rate by,

(3.126)

Wall functions for momentum and other turbulence quantities and transported variables are not used with this model.

## 3.2.6.3. RNG RANS Model

The RNG model was derived using a rigorous statistical decomposition of the velocity field called renormalization group (RNG) theory. This model has several significant benefits over the standard model, including improved accuracy for rapidly strained flows, swirling flows, and low Reynolds number flows, without additional modifications. Additionally, values for the model constants are derived analytically rather than being evaluated empirically. Papageorgakis and Assanis [20] describe the version of the RNG model as implemented here.

The same turbulent kinetic energy equation as in the standard model, (3.118), is used for the RNG equation. The turbulent kinetic energy dissipation rate equation is the same as (3.119), with the addition of a single source term on the right-hand-side of the equation,

(3.127)

where , , and are model constants, and

(3.128)

As with the standard model, the turbulent eddy viscosity is then given by the Prandtl-Kolmogorov relationship,

(3.129)

## 3.2.6.4. RANS Model

Durbin [21] introduced a method for handling the wall region without using either wall functions or damping functions. In his method a fine grid is required near the wall (e.g., the first grid point is typically within one dimensionless unit of distance from the wall where the coordinate normal to the wall is nondimensionalized with the inner scale for a turbulent boundary layer, at the first grid point, where is the friction velocity, ). The model employs two transport equations in addition to slightly modified and equations to account for the nonhomogeneous region near the wall. The eddy viscosity is formulated using the component of turbulent kinetic energy normal to the wall for velocity scaling (instead of using as in the standard model).

The turbulent kinetic energy, , is given by (3.118) while the dissipation rate of turbulent kinetic energy, , is given by

(3.130)

The time scale, , is the usual time scale , away from the wall region; however, near the wall, if becomes smaller than the Kolmogorov time scale , then the latter is used for . This is formally stated by

(3.131)

(3.132)

where

(3.133)

and the modified constant, , is given by

(3.134)

The model includes a transport equation for ,

(3.135)

An elliptic relaxation model equation is formulated to solve for the variable in the above equation. The purpose of the elliptic relaxation model is to account for nonlocal effects such as wall blocking; the equation is given by

(3.136)

Finally, the turbulent eddy viscosity is given by

(3.137)

## 3.2.6.5. RANS Model

The turbulence model and its variants are similar in structure to the models. However, instead of computing the turbulent dissipation rate directly, the model models the transport the reciprocal of a turbulent timescale referred to as the turbulent frequency. This quantity, , can be related to the turbulent dissipation by

(3.138)

The the transport equations are given by the 2006 model, (Wilcox [22]),

(3.139)

(3.140)

The user is to note the above standard for writing the effective diffusive flux coefficient. The model also has a number of adjustable parameters: , , , , , and . The constant is given by,

(3.141)

where

(3.142)

The value of is as follows:

(3.143)

The production term is the same as in . Typically limiters are used to prevent it from exceeding the dissipation rate by too large an amount. Although the 2006 description does not speak of production limiters, other sources that use the 2006 model do, i.e.

(3.144)

The value of is expected to be a user specified quantity (see input file manual for more details). In general, this term is defaulted to a very high number.

The eddy viscosity is

(3.145)

where is,

(3.146)

## 3.2.6.6. Shear Stress Transport (SST)

It has been observed that standard 1998 models display a strong sensitivity to the free stream value of . To remedy, this, an alternative set of transport equations have been used that are based on smoothly blending the model near a wall with away from the wall (see Mentor [23]). Because of the relationship between and , the transport equations for turbulent kinetic energy and dissipation can be transformed into equations involving and . Aside from constants, the transport equation for is unchanged. However, an additional cross-diffusion term is present in the equation. Blending is introduced by using smoothing which is a function of the distance from the wall, . The transport equations for the Mentor 2003 model ( [23]) are provided by the following:

(3.147)

(3.148)

The model coefficients, , , and must also be blended, which is represented by

(3.149)

where , , , , , , and .

The blending function is given by

(3.150)

where

(3.151)

The final parameter is

(3.152)

In the 2003 SST model description, the production term is expected to be limited:

(3.153)

The value of is expected to be a user specified quantity (see input file manual for more details). In general, this term is defaulted to a very high number.

An important component of the SST model is the different expression used for the eddy viscosity,

(3.154)

where is another blending function given by

(3.155)

The final parameter is

(3.156)

## 3.2.6.7. Standard Smagorinsky LES Model

The standard Smagorinsky LES closure model approximates the subgrid turbulent eddy viscosity using a mixing length-type model, where the LES grid filter size provides a natural length scale. The subgrid eddy viscosity is modeled simply as (Smagorinsky [24])

(3.157)

where the strain rate tensor magnitude is defined as . The constant coefficient typically varies between and and should be carefully tuned to match the problem being solved (Rogallo and Moin [25]). It is assigned a value of here.

Although this model is desirable due to its simplicity and efficiency, care should be taken in its application. It is known to predict subgrid turbulent eddy viscosity proportional to the shear rate in the flow, independent of the local turbulence intensity. Non-zero subgrid turbulent eddy viscosity is even predicted in completely laminar regions of the flow, sometimes even preventing a natural transition to turbulence. Therefore, this model should only be used when this behavior will not adversely affect results.

## 3.2.6.8. Dynamic Smagorinsky LES Model

As mentioned in the previous section, the standard Smagorinsky model requires careful tuning of the constant model coefficient for the particular problem being simulated, and it is often overly-dissipative due to its inability to adapt to the local turbulent environment. Germano et al. [26] developed an improvement over the standard Smagorinsky model, where the coefficient is dynamically calculated based on the local turbulence field. A generalization of this method for variable-density flow is used here (Moin et al. [14]).

Similar to the standard Smagorinsky LES closure model, the subgrid eddy viscosity is modeled by the mixing length approximation

(3.158)

where the strain rate tensor magnitude is defined as . The coefficient is dynamically evaluated by taking advantage of scale similarity in the inertial range of the turbulence spectrum, near the minimum resolved scales. This is done by introducing a “test filter” which is identical to the grid filter defined in (3.72) except for having a larger filter size denoted by . The test filter of variable is denoted by .

The previously-defined subgrid stress tensor can be rewritten as

(3.159)

and an analogous larger-scale “subtest” stress can be analogously defined as

(3.160)

where the notation denotes resolved quantities that have been passed through the test filter. These two stresses can be related to each other through the algebraic identity of Germano [27],

(3.161)

(3.162)

Note that the right-hand side of (3.162) is completely computable in terms of resolved quantities.

By modeling the two stresses in (3.161) and equating them to (3.162), the model coefficient can be dynamically evaluated. The subtest stress is modeled analogously to the subgrid stress, as

(3.163)

(3.164)

where is assumed to be the same at both scales. The test-filtered strain rate tensor is defined similar to as

(3.165)

Notice that when the modeled forms of and are substituted into (3.161), appears inside a test filtering operation. Formally solving this system of equations for requires the expensive proposition of solving an additional set of coupled integro-differential equations (Ghosal et al. [28]). Alternatively, it is common practice to remove from the test filter with the assumption that it is varying slowly over distances on the order of the test filter size. This greatly simplifies calculations, although it yields a system of over-determined equations for this single constant. The square of the error involved in this approximation is , where

(3.166)

(3.167)

Minimizing this error in a least-squares fashion yields an expression for the modeled Smagorinsky coefficient (Lilly [29]),

(3.168)

that can be used directly in (3.158) for the subgrid turbulent eddy viscosity.

Due to the above simplifications, the model constant can sometimes fluctuate wildly to both large positive and negative values. These fluctuations can possibly lead to numerical instability, so they must be controlled. A common solution, and one that is taken here, is to pass the numerator and denominator of (3.168) through the test filter, yielding

(3.169)

This can be crudely justified by recognizing that was already assumed to vary slowly over distances equal to the test filter size, so that this filtering operation is simply enforcing that assumption.

This form of the dynamic Smagorinsky closure model allows energy backscatter, which is an intermittent transfer of turbulent kinetic energy from small scales to larger scales rather than the typical cascade from large to small scales. While backscatter can occur in real turbulent flows, the predicted negative eddy viscosities of the dynamic Smagorinsky model are more often attributable to model errors than to a real physical backscatter process. This can easily destabilize a simulation, so negative eddy viscosity is disallowed in the present formulation.

The only free parameter in the dynamic Smagorinsky closure model is the ratio between the test and grid filter sizes, . Solutions are fairly insensitive to the choice of , although values of around are usually considered optimal (Germano et al. [26]). This ratio is dictated by the box filter formulation used in Fuego and the mesh topology selected by the user. The test filter volume for a particular CVFEM node is defined as the volume of all surrounding finite elements that contain that node. (See Numerics for more information about the CVFEM formulation.) On uniform hexahedral and uniform quadrilateral meshes, the test filter ratio will have a value of . The ratio will be around for uniform tetrahedral meshes and around for uniform triangular meshes, which are still reasonable values.

## 3.2.6.9. Subgrid-Scale Kinetic Energy One-Equation LES Model

The subgrid scale kinetic energy one-equation turbulence model, or model, represents a simple LES closure model. The transport equation for subgrid turbulent kinetic energy is given by

(3.170)

The production of subgrid turbulent kinetic energy, , is modeled by (3.121) while the dissipation of turbulent kinetic energy, , is given by

(3.171)

where the grid filter length, , is given in terms of the grid cell volume by

(3.172)

The subgrid turbulent eddy viscosity is then provided by

(3.173)

where the values of and are 0.845 and 0.0856, respectively.

### 3.2.6.9.1. Wall-Resolved LES Approach

The baseline one-equation model is designed to support applications where wall boundary layers are not resolved, i.e, wall-modeled LES (WMLES). For simulations that include wall-resolved boundary layer meshing approaches, i.e., near wall normalized wall unit, less than unity, additional transport terms and blending functions must be added. These additional contributions ensure that quantities such as turbulent viscosity and normalized turbulence kinetic energy scale cubically and quadratically with . The approach followed in Fuego to support the wall-resolved LES (WRLES) use case is to support the model developed by Inagaki [30].

The one-equation form for the Inagaki model is very similar to the baseline model with only a few select changes that allow proper scaling of the subgrid-scale turbulence kinetic energy near a wall boundary. In this model, the dissipation rate, now includes an added dissipation term,

(3.174)

while the turbulent viscosity is modified via a blending coefficient ,

(3.175)

The functional form for provided in Inagaki [30] is,

(3.176)

where the model constants and are 20 and 2, the former of which closely replicates Van Driest scaling while the value of 2 provides the functional form of (here, is the minimum distance to the wall, a value that is internally computed within Fuego).

The functions and are,

(3.177)

and

(3.178)

that include the static model constant , which is specified to be a value of 4. As implied by (3.174), is the modified dissipation rate that includes the augmented dissipation by gradients of ,

(3.179)

The model constants for and in Inagaki are provided as 0.835 and 0.054, respectively. Validation test cases demonstrate a somewhat weak dependence on the precise values of the model constants as with the diffusive flux vector scaling, , which was provided as 0.54. Note that this model is not designed to be activated in concert with the dynamic procedure to follow.

## 3.2.6.10. Dynamic Subgrid-Scale Kinetic Energy One-Equation LES Model

Similar to the dynamic Smagorinsky model in Dynamic Smagorinsky LES Model, a dynamic version is developed for the subgrid kinetic energy model. The standard version with fixed coefficients over-predicts turbulent viscosity while the dynamic version is known to offer a better predictability. In Fuego, and are calculated dynamically which is considered to be a standard approach for the dynamic model [31]. The concept of “test filter” is identical to that of the dynamic Smagorinsky model in Dynamic Smagorinsky LES Model. Subgrid-scale kinetic energy for grid-filter and test-filter levels are

(3.180)

(3.181)

Exact form of the dissipation in (3.170) is

(3.182)

where , , and . Meanwhile, dissipates by both molecular and turbulent viscosities of the grid-filtered level since the quantity is fully resolved in the test-filter level [31].

(3.183)

Using scale similarity, (3.171) applies to the test-filter level as

(3.184)

and therefore, is calculated by

(3.185)

The other coefficient, , is computed similarly to the (3.168) as

(3.186)

where is defined identically to (3.166) and is simplified by

(3.187)

Note that dynamic subgrid kinetic energy model does not require an additional filtering as in (3.169). The method will predict negative values of turbulent viscosity, similar to dynamic smagorinsky. Numerical issues can arise if model predicts an overly negative value of the turbulent viscosity for an overly long timescale. Because of this, clipping is suggested [31]

(3.188)

with a value of of unity, such that the effective viscosity remains non-negative. In Fuego, we allow this value to be user-defined with a default value of zero. A sufficiently large value of will effectively result in using the directly modeled turbulent viscosity.

## 3.2.6.11. Buoyancy Models for the Production Rate

There are three supported models that augment the production of turbulent kinetic energy via buoyancy contributions, buoyant vorticity generation [32], Rodi’s [33], and de Ris’ [34] buoyancy term.

The buoyant vorticity generation model has been developed and validated by Sandia National Laboratories group 9132 for use in large scale buoyant plumes. The model attempts to augment the production of turbulent kinetic energy by adding a source term, to both the turbulent kinetic energy and dissipation rate equation that is related to the baroclinic torque,

(3.189)

of the model.

The buoyancy model of Rodi is given by

(3.190)

De Ris’ buoyancy model offers two versions - flaming

(3.191)

and non-flaming.

(3.192)

Default value of the user-defined coefficient is 0.01. Note that ambient and flame density, rather than local density, matters on the flaming version.

In each model, derivatives are evaluated at the subcontrol volume center while the property values are lumped.

The right hand side of the turbulent kinetic energy equation for all model is . For the dissipation rate equation, the source term is for the buoyant vorticity generation model while it is otherwise. Recall that the inverse time scale is determined by the turbulence model of choice, i.e., for the standard model and provided in (3.131) for the model.

Note that the use of the buoyancy models has not been evaluated with the model.

## 3.2.6.12. Turbulence closure model constants

For each of the aforementioned turbulence closure models, there are several constant coefficients which may be modified by the user in the input deck. Table 3.1, Table 3.2, Table 3.3, and Table 3.4 list these parameters, their mapping to input deck names, and default values. Each of these default values may be modified by the user by specifying the respective `Turbulence Model Parameter` line in the `Global Constants` block under the Sierra domain.

Table 3.1 Constant parameters for turbulence models.
Turbulence Model | Symbol | User Input Name | Default Value
--- | --- | --- | ---
Standard |  | `Cmu` | 0.09
| `Ceps_1` | 1.44
| `Ceps_2` | 1.92
| `Cchi` | 2.0
| `Sigma_K` | 1.0
| `Sigma_E` | 1.3
Low Reynolds |  | `Cmu` | 0.09
| `Ceps_1` | 1.44
| `Ceps_2` | 1.92
| `Sigma_K` | 1.0
| `Sigma_E` | 1.3
| `Amu` | 3.4
RNG |  | `Cmu` | 0.0837
| `Ceps_1` | 1.42
| `Ceps_2` | 1.68
| `Sigma_K` | 0.7194
| `Sigma_E` | 0.7194
v^2 - f |  | `Cmu` | 0.22
| `Ceps_1` | 1.4
| `Ceps_2` | 1.9
| `Sigma_K` | 1.0
| `Sigma_E` | 1.0
| `CF_1` | 0.4
| `CF_2` | 0.3
| `Alpha` | 0.6
| `Nseg` | 6.0
| `CL` | 0.23
| `Ceta` | 70.0

Table 3.2 Constant parameters for turbulence models.
Turbulence Model | Symbol | User Input Name | Default Value
--- | --- | --- | ---
|  | `Beta_Zero` | 0.0708
| `Beta_Star` | 0.09
| `Sigma_K` |
| `Sigma_W` | 0.5
| `Gamma` |
| `Clim` |
SST |  | `A_One` | 0.31
| `Beta_One` | 0.075
| `Beta_Two` | 0.0828
| `Beta_Star` | 0.09
| `Gamma_One` |
| `Gamma_Two` | 0.44
| `Sigma_K_One` | 0.85
| `Sigma_K_Two` | 1.0
| `Sigma_W_One` | 0.5
| `Sigma_W_Two` | 0.856

Table 3.3 Constant parameters for LES turbulence models.
Turbulence Model | Symbol | User Input Name | Default Value
--- | --- | --- | ---
One-equation |  | `Cv` | 0.5
| `Ceps` | 0.845
| `Cmueps` | 0.0856
Standard Smagorinsky |  | `Cv` | 0.5
| `Cs` | 0.17
Dynamic Smagorinsky |  | `Cs` | 0.17

Table 3.4 Constant parameters for miscellaneous turbulence models. Default values may be changed using the model parameters input.
Model | Symbol | User Input Name | Default Value
--- | --- | --- | ---
Buoyant vorticity generation |  | `Cbvg` | 0.35
| `Ceps_3` | 0.0
Rodi’s source term |  | `C_eps4` | 0.0
EDC laminar limit |  | `Cgammalam` | 2.0
| `Ctaulam` | 0.02
| `Clamtrans` | 40.0