3.2.6. Turbulence Closure Models

The Favre-filtered turbulent flow equations of the previous section have been modeled in terms of \mu_t, 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 k - \epsilon RANS Model

The standard k - \epsilon 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 \mu_t. The turbulent kinetic energy, k, and the dissipation rate of turbulent kinetic energy, \epsilon, are given by (Gran et al. [16])

(3.118)\int {{\partial \bar{\rho}{k}} \over {\partial t}} {\rm d}V + \int \bar{\rho}{k} \tilde{u}_j n_j {\rm d}S = \int {{\mu_t} \over {\sigma_k}} {{\partial {k}} \over {\partial x_j}} n_j {\rm d}S + \int \left(P_k - \bar{\rho} {\epsilon}\right) {\rm d}V

(3.119)\int {{\partial \bar{\rho} {\epsilon}} \over {\partial t}} {\rm d}V + \int \bar{\rho} {\epsilon} \tilde{u}_j n_j {\rm d}S = \int {{\mu_t} \over {\sigma_\epsilon}} {{\partial {\epsilon}} \over {\partial x_j}} n_j {\rm d}S + \int {\epsilon \over {k}} \left(C_{\epsilon 1} P_k - C_{\epsilon 2} \bar{\rho} {\epsilon}\right) {\rm d}V ,

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

(3.120)P_k \equiv -\overline{\rho u_i'' u_j''} {{\partial \tilde{u}_i} \over {\partial x_j}},

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

(3.121)P_k &= \mu_t \left( {{\partial \tilde{u}_i} \over {\partial x_j}} + {{\partial \tilde{u}_j} \over {\partial x_i}} \right) {{\partial \tilde{u}_i} \over {\partial x_j}} - {2 \over 3} \left( \bar{\rho} {k} + \mu_t{{\partial \tilde{u}_k} \over {\partial x_k}} \right) {{\partial \tilde{u}_m} \over {\partial x_m}} &= \left[ 2 \mu_t \left( \tilde{S}_{ij} - \frac{1}{3} \tilde{S}_{kk} \delta_{ij} \right) - \frac{2}{3} \bar{\rho} \tilde{k} \delta_{ij} \right] \frac{\partial \tilde{u}_i}{\partial \tilde{x}_j}.

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

(3.122)\mu_t = {\rm C}_{\mu} \bar{\rho} {{k} {\tau}}.

where \tau = min( {k \over \epsilon}, \Delta t_f). The filter time, \Delta t_f is provided by the temporally filtered Navier Stokes model (Tieszen et al. [17]). The parameters C_{\epsilon 1}, C_{\epsilon 2}, \sigma_k, and \sigma_\epsilon 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 k - \epsilon RANS Model

In the case of the low Reynolds number turbulent flows, the standard k - \epsilon 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)S_k^{lr} = -2{\mu}\left( \partial \sqrt{k} \over {\partial x_j} \right)^2.

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

(3.124)S_{\tilde \epsilon}^{lr} = 2\frac{\mu \mu^T}{\rho} \left(\frac{\partial S} {\partial x_k } \right)^2,

where the strain-rate magnitude is defined by S = \sqrt{2 S_{ij}S_{ij}}.

The constants in the dissipation rate equation are modified by damping coefficients, C_{\epsilon_1} = f_1 C_{\epsilon_1} and C_{\epsilon_2} = f_2 C_{\epsilon_2}, where f_1 is unity and f_2 = 1-0.3 e^{-Re_t^2}.

In this formulation, the eddy viscosity relationship defined by,

(3.125)\mu_t = f_{\mu} {\rm C}_{\mu} \bar{\rho} k \tau,

where the dampening function is f_\mu = e^{-\frac{A_\mu}{(1+Re_t/50)^2}} and \tau is now a function of \tilde \epsilon. The constant A_\mu is 3.4 with Re_t = \frac{\rho k^2}{\mu \tilde{\epsilon}}.

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

(3.126)\epsilon = \tilde{\epsilon} + 2\frac{\mu}{\rho}\left( \partial \sqrt{k} \over {\partial x_j} \right)^2.

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

3.2.6.3. RNG k - \epsilon RANS Model

The RNG k - \epsilon 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 k - \epsilon 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 k - \epsilon model as implemented here.

The same turbulent kinetic energy equation as in the standard k - \epsilon model, (3.118), is used for the RNG k - \epsilon 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)S_\epsilon^\mathrm{RNG} = - \frac{ C_\mu \eta^3 (1 - \eta/\eta_o) }{ 1 + \beta \eta^3 } \frac{\epsilon^2}{k},

where C_\mu, \beta, and \eta_o are model constants, and

(3.128)\eta = (2 \tilde{S}_{ij} \tilde{S}_{ij})^\frac{1}{2} \frac{k}{\epsilon}.

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

(3.129)\mu_t = C_{\mu} \bar{\rho} {{k} {\tau}}.

3.2.6.4. v^2-f 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, y^{+} = y u_\tau / \nu < 1 at the first grid point, where u_\tau is the friction velocity, \sqrt{\tau_w / \bar{\rho}}). The model employs two transport equations in addition to slightly modified k and \epsilon 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 \sqrt{k} as in the standard k - \epsilon model).

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

(3.130)\int {{\partial \bar{\rho} {\epsilon}} \over {\partial t}} {\rm d}V + \int \bar{\rho} {\epsilon} \tilde{u}_j n_j {\rm d}S = \int {{\mu_t} \over {\sigma_\epsilon}} {{\partial {\epsilon}} \over {\partial x_j}} n_j {\rm d}S + \int {1 \over {T}} \left(C_{\epsilon 1}' P_k - C_{\epsilon 2} \bar{\rho} {\epsilon}\right) {\rm d}V.

The time scale, T, is the usual time scale k/\epsilon, away from the wall region; however, near the wall, if k/\epsilon becomes smaller than the Kolmogorov time scale \sqrt{\nu/\epsilon}, then the latter is used for T. This is formally stated by

(3.131)T = \min \left[ T_1, {\alpha \over {2 \sqrt{3}}} {k \over {\overline{v^2} C_\mu \sqrt{\tilde{S}^2}}} \right]

(3.132)T_1 = \max \left[ {k \over \epsilon}, 6 \sqrt{\nu \over \epsilon} \right],

where

(3.133)\tilde{S}^2 = \tilde{S}_{ij} \tilde{S}_{ij} = {1 \over 4} \left( {{\partial \tilde{u}_i} \over {\partial x_j}} + {{\partial \tilde{u}_j} \over {\partial x_i}} \right) \left( {{\partial \tilde{u}_i} \over {\partial x_j}} + {{\partial \tilde{u}_j} \over {\partial x_i}} \right)

and the modified constant, C_{\epsilon_1}', is given by

(3.134)C_{\epsilon_1}'= C_{\epsilon_1} \left( 1 + 0.045\sqrt{k/\overline{v^2}} \right).

The model includes a transport equation for \overline{v^2},

(3.135){{\partial \bar{\rho} \overline{v^2}} \over {\partial t}} + {{\partial \bar{\rho} \tilde{u}_j \overline{v^2}} \over {\partial x_j}} = {{\partial} \over {\partial x_j}} \left[ \left( \mu + \mu_t \right) {{\partial \overline{v^2}} \over {\partial x_j}} \right] + \bar{\rho} k f - { \bar{\rho} N \overline{v^2} \over T_1 } .

An elliptic relaxation model equation is formulated to solve for the variable f 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)f - L^2 {{\partial} \over {\partial x_j}} \left( {{\partial f} \over {\partial x_j}} \right) = C_1 {{\left(2/3 - \overline{v^2} / k \right)} \over {T_1}} + C_2 2 \nu_t {{\tilde{S}^2} \over k} + \left( N - 1 \right) {{\overline{v^2} / k} \over {T_1}} .

Finally, the turbulent eddy viscosity is given by

(3.137)\mu_t = C_\mu \bar{\rho} \overline{v^2} {\tau} .

3.2.6.5. k - \omega RANS Model

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

(3.138)\epsilon = \beta^* k \omega.

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

(3.139)\int{\partial \bar{\rho} k \over \partial t}\text{d}V + \int \bar{\rho} k\tilde{u}_{j} n_{j} \text{d} S = \int ({\mu+\sigma_k\frac{\bar{\rho} k}{\omega}}){\partial k \over \partial x_{j}} n_{j} \text{d} V + \int \left(P_{k}^{\omega} - \beta^* \bar{\rho} k \omega\right) \text{d} V,

(3.140)\int {\partial \bar{\rho} \omega \over \partial t}\text{d} V + \int \bar{\rho} \omega \tilde{u}_{j} n_{j} \text{d} S = \int ({\mu+\sigma_{\omega}\frac{\bar{\rho} k}{\omega}} ){\partial \omega \over \partial x_{j}} n_{j} \text{d} V + \int \left(\gamma {\omega \over k} P_{k}^{\omega} - \beta \bar{\rho} \omega^{2} + \frac{\bar{\rho} \sigma_d} {\omega} {\partial k \over \partial x_j} {\partial \omega \over \partial x_j} \right) \text{d} V.

The user is to note the above standard for writing the effective diffusive flux coefficient. The model also has a number of adjustable parameters: \beta_0 = 0.0708, \beta^* = 0.09, \gamma = \frac{13}{25}, C_{lim} = \frac{7}{8}, \sigma_{k} = 0.6, and \sigma_{\omega} = 0.5. The constant \beta is given by,

(3.141)\beta = \beta_0 f_{\beta}

where

(3.142)f_{\beta} = \frac{1+85 \chi_{\omega}}{1+100 \chi_{\omega}}

The value of \chi_{\omega} is as follows:

(3.143)\chi_{\omega} = | \frac{\Omega_{ij}\Omega_{jk} S_{ki}}{(\beta^*\omega)^3} |

The production term is the same as in k-\epsilon. 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)P_{k}^{\omega} = \max\left(P_{k}, 10 \bar{\rho} k \omega \right).

The value of 10 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)\mu_{T} = \bar{\rho}{k \over \hat{\omega}}.

where \hat{\omega} is,

(3.146)\hat{\omega} = \max(\omega, C_{lim} \sqrt\frac{2 S_{ij} S_{ij}}{\beta^*}).

3.2.6.6. Shear Stress Transport (SST)

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

(3.147)\int{\partial \bar{\rho} k \over \partial t}\text{d}V + \int \bar{\rho} k\tilde{u}_{j} n_{j} \text{d} S = \int {(\mu + \hat \sigma_k \mu_{t})} {\partial k \over \partial x_{j}} n_{j} + \int \left(P_{k}^{\omega} - \beta^* \bar{\rho} k \omega\right) \text{d} V,

(3.148)\int {\partial \bar{\rho} \omega \over \partial t}\text{d} V + \int \bar{\rho} \omega \tilde{u}_{j} n_{j} \text{d} S = \int {(\mu + \hat\sigma_{\omega} \mu_{t})} {\partial \omega \over \partial x_{j}} n_{j} + \int {2(1-F) \frac{\bar{\rho}\sigma_{\omega2}} {\omega} {\partial k \over \partial x_j} {\partial \omega \over \partial x_j} } \text{d} V \newline + \int \left(\frac{\hat\gamma}{\nu_t} P_{k}^{\omega} - \hat \beta \bar{\rho} \omega^{2}\right) \text{d} V.

The model coefficients, \hat\sigma_k, \hat\sigma_{\omega}, \hat\gamma and \hat\beta must also be blended, which is represented by

(3.149)\hat \phi = F\phi_1+ (1-F)\phi_2.

where \sigma_{k1} = 0.85, \sigma_{k2} = 1.0, \sigma_{\omega1} = 0.5, \sigma_{\omega2} = 0.856, \gamma_1 = \frac{5}{9}, \gamma_2 = 0.44, \beta_1 = 0.075 and \beta_2 = 0.0828.

The blending function is given by

(3.150)F = \tanh(A_{1}^{4}),

where

(3.151)A_{1} = \min \left( \max \left( {\sqrt{k} \over \beta^* \omega y}, {500 \mu \over \bar{\rho} y^{2} \omega}\right), {4 \bar{\rho} \sigma_{\omega2} k \over CD_{k\omega} y^{2}} \right).

The final parameter is

(3.152)CD_{k\omega} = \max \left( 2 \bar{\rho} \sigma_{\omega2} \frac{1}{\omega} {\partial k \over \partial x_{j}} {\partial \omega \over \partial x_{j}}, 10^{-10} \right).

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

(3.153)P_{k}^{\omega} = \max\left(P_{k}, 10 \bar{\rho} k \omega \right).

The value of 10 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)\mu_{t} = \frac {a_1 \bar{\rho} k} {\max\left( a_1 \omega, S F_2 \right) },

where F_2 is another blending function given by

(3.155)F_2 = \tanh(A_{2}^{2}).

The final parameter is

(3.156)A_{2} = \max\left({2 \sqrt{k} \over \beta^* \omega y}, {500 \mu \over \bar{\rho} \omega y^{2}} \right).

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 \Delta provides a natural length scale. The subgrid eddy viscosity is modeled simply as (Smagorinsky [24])

(3.157)\mu_t = \rho \left(C_s \Delta \right)^2 |\tilde{S}|,

where the strain rate tensor magnitude is defined as |\tilde{S}| \equiv (2 \tilde{S}_{ij} \tilde{S}_{ij})^{{1} \over {2}}. The constant coefficient C_s typically varies between 0.1 and 0.24 and should be carefully tuned to match the problem being solved (Rogallo and Moin [25]). It is assigned a value of 0.17 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 C_s 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)\mu_t = C_R \bar{\rho} \Delta^2 |\tilde{S}|,

where the strain rate tensor magnitude is defined as |\tilde{S}| \equiv (2 \tilde{S}_{ij} \tilde{S}_{ij})^{{1} \over {2}}. The coefficient C_R 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 \hat{\Delta}. The test filter of variable \phi is denoted by \hat{\phi}.

The previously-defined subgrid stress tensor can be rewritten as

(3.159)\tau_{u_i u_j} &\equiv - ( \bar{\rho} \widetilde{u_i u_j} - \bar{\rho} \tilde{u}_i \tilde{u}_j ) & = - \left( \overline{\rho u_i u_j} - \frac{ \overline{\rho u_i}\:\overline{\rho u_j} }{ \bar{\rho} } \right)

and an analogous larger-scale “subtest” stress T_{u_i u_j} can be analogously defined as

(3.160)T_{u_i u_j} \equiv - \left( \widehat{\overline{\rho u_i u_j}} - \frac{ \widehat{\overline{\rho u_i}}\: \widehat{\overline{\rho u_j}} }{ \hat{\bar{\rho}} } \right),

where the \hat{\bar{()}} 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)L_{u_i u_j} \equiv T_{u_i u_j} - \widehat{\tau_{u_i u_j}}

(3.162)L_{u_i u_j} = - \left( \widehat{\bar{\rho} \tilde{u}_i \tilde{u}_j} - \frac{ \widehat{\bar{\rho} \tilde{u}_i}\: \widehat{\bar{\rho} \tilde{u}_j} }{ \hat{\bar{\rho}} } \right).

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 C_R can be dynamically evaluated. The subtest stress is modeled analogously to the subgrid stress, as

(3.163)\tau_{u_i u_j} \approx 2 C_R \bar{\rho} \Delta^2\, |\tilde{S}| \left( \tilde{S}_{ij} - \frac{1}{3}\tilde{S}_{kk} \delta_{ij} \right)

(3.164)T_{u_i u_j} \approx 2 C_R \hat{\bar{\rho}} \hat{\Delta}^2\, |\hat{\tilde{S}}| \left( \widehat{ \tilde{S}_{ij}} - \frac{1}{3} \widehat{ \tilde{S}_{kk}} \delta_{ij} \right),

where C_R is assumed to be the same at both scales. The test-filtered strain rate tensor is defined similar to |\tilde{S}| as

(3.165)|\hat{\tilde{S}}| \equiv \left(2 \widehat{\tilde{S}_{ij}} \widehat{\tilde{S}_{ij}} \right)^\frac{1}{2} .

Notice that when the modeled forms of \tau_{u_i u_j} and T_{u_i u_j} are substituted into (3.161), C_R appears inside a test filtering operation. Formally solving this system of equations for C_R 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 C_R 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 Q = (L_{ij} - C_R M_{ij})^2, where

(3.166)L_{u_i u_j} = - \left( \widehat{\bar{\rho} \tilde{u}_i \tilde{u}_j} - \frac{ \widehat{\bar{\rho} \tilde{u}_i}\, \widehat{\bar{\rho} \tilde{u}_j} }{ \hat{\bar{\rho}} } \right)

(3.167)M_{u_i u_j} = 2 \hat{\bar{\rho}} \hat{\Delta}^2 |\hat{\tilde{S}}| \left( \widehat{\tilde{S}_{ij}} - \frac{1}{3} \widehat{ \tilde{S}_{kk}} \delta_{ij} \right) - \widehat{ 2 \bar{\rho} \Delta^2\, |\tilde{S}| \left(\tilde{S}_{ij} - \frac{1}{3}\tilde{S}_{kk} \delta_{ij} \right) }.

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

(3.168)C_R = \frac{ L_{u_i u_j} M_{u_i u_j} }{ M_{u_i u_j} M_{u_i u_j} },

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

Due to the above simplifications, the model constant C_R 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)C_R = \frac{ \widehat{L_{u_i u_j} M_{u_i u_j}} }{ \widehat{M_{u_i u_j} M_{u_i u_j}} }.

This can be crudely justified by recognizing that C_R 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, \alpha = \hat{\Delta}/\Delta. Solutions are fairly insensitive to the choice of \alpha, although values of around \alpha = 2 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 2.0. The ratio will be around 1.59 for uniform tetrahedral meshes and around 1.73 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 Ksgs model, represents a simple LES closure model. The transport equation for subgrid turbulent kinetic energy is given by

(3.170)\int {{\partial \bar{\rho}{k^\mathrm{sgs}}} \over {\partial t}} {\rm d}V + \int \bar{\rho}{k^\mathrm{sgs}} \tilde{u}_j n_j {\rm d}S = \int {{\mu_t} \over {\sigma_k}} {{\partial {k^\mathrm{sgs}}} \over {\partial x_j}} n_j {\rm d}S + \int \left(P_k^\mathrm{sgs} - D_k^\mathrm{sgs}\right) {\rm d}V.

The production of subgrid turbulent kinetic energy, P_k^\mathrm{sgs}, is modeled by (3.121) while the dissipation of turbulent kinetic energy, D_k^\mathrm{sgs}, is given by

(3.171)D_k^\mathrm{sgs} = \bar{\rho} C_{\epsilon} { { {k^\mathrm{sgs}}^{{3}\over {2}} } \over { \Delta} },

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

(3.172)\Delta = V^{{1}\over{3}}.

The subgrid turbulent eddy viscosity is then provided by

(3.173)\mu_t = \bar{\rho} C_{\mu_{\epsilon}} \Delta {k^\mathrm{sgs}}^{{1} \over {2}},

where the values of C_{\epsilon} and C_{\mu_{\epsilon}} 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, y^+ 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 y^+. 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, D_k^\mathrm{sgs} now includes an added dissipation term,

(3.174)D_k^\mathrm{sgs} = \bar{\rho} \epsilon_{sgs} = \bar{\rho} C_{\epsilon} { { {k_{sgs}}^{{3}\over {2}} } \over { \Delta} } + 2 \mu \frac{\partial \sqrt{k_{sgs}}}{\partial x_j}\frac{\partial \sqrt{k_{sgs}}}{\partial x_j},

while the turbulent viscosity is modified via a blending coefficient f_\mu,

(3.175)\mu_t = f_\mu \bar{\rho} C_{\mu_{\epsilon}} \Delta {k^\mathrm{sgs}}^{{1} \over {2}}.

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

(3.176)f_\mu = 1 - exp \left[ - \left( \left(\frac{y'_\epsilon}{A_o} \right)^{\frac{2}{3}} \right)^{B_o} \right]

where the model constants A_o and B_o are 20 and 2, the former of which closely replicates Van Driest scaling while the B_o value of 2 provides the functional form of f_{\mu} \propto y^2 (here, y is the minimum distance to the wall, a value that is internally computed within Fuego).

The functions y'_\epsilon and u'_{\epsilon} are,

(3.177)y'_\epsilon = \frac{y u'_{\epsilon}}{\nu}

and

(3.178)u'_{\epsilon} = \left( \nu \epsilon_{sgs} \right)^{\frac{1}{4}} \sqrt {C_l \frac{y}{\Delta}},

that include the static model constant C_l, which is specified to be a value of 4. As implied by (3.174), \epsilon_{sgs} is the modified dissipation rate that includes the augmented dissipation by gradients of \sqrt {k_{sgs}},

(3.179)\epsilon_{sgs} = C_{\epsilon} { { {k_{sgs}}^{{3}\over {2}} } \over { \Delta} } + 2 \nu \frac{\partial \sqrt{k_{sgs}}}{\partial x_j}\frac{\partial \sqrt{k_{sgs}}}{\partial x_j}.

The model constants for C_{\epsilon} and C_{\mu_{\epsilon}} 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, \sigma^k, 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, C_{\epsilon} and C_{\mu_{\epsilon}} are calculated dynamically which is considered to be a standard approach for the dynamic Ksgs 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)k^\mathrm{sgs} = \frac{1}{2} \left( \widetilde{u_k u_k} - \tilde{u_k} \tilde{u_k} \right) ,

(3.181)k^\mathrm{test} = \frac{1}{2} \left( \widehat{\tilde{u_k} \tilde{u_k}} - \hat{\tilde{u_k}} \hat{\tilde{u_k}} \right) .

Exact form of the dissipation D_k^\mathrm{sgs} in (3.170) is

(3.182)D_k^\mathrm{sgs} = 2\bar{\mu} \left[ \widetilde{S_{ij}^* D_{ij}^*} - \tilde{S_{ij}^*} \tilde{D_{ij}^*} \right]

where S_{ij}^* = {S}_{ij} - \frac{1}{3}{S}_{kk} \delta_{ij}, D_{ij}^* = {D}_{ij} - \frac{1}{3}{D}_{kk} \delta_{ij}, and D_{ij}=\partial{u_i} / \partial{x_j}. Meanwhile, k^\mathrm{test} 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)D_k^\mathrm{test} = 2(\bar{\mu}+\mu_t) \left[ \widehat{ \tilde{S_{ij}^*} \tilde{D_{ij}^*} } - \widehat{\tilde{S_{ij}^*}} \widehat{\tilde{D_{ij}^*}} \right]

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

(3.184)D_k^\mathrm{test} = C_{\epsilon} \hat{\bar{\rho}} { { {k^\mathrm{test}}^{{3}\over {2}} } \over { \hat{\Delta}} },

and therefore, C_\epsilon is calculated by

(3.185)C_\epsilon = \frac{ 2(\bar{\mu}+\mu_t) \left[ \widehat{ \tilde{S_{ij}^*} \tilde{D_{ij}^*} } - \widehat{\tilde{S_{ij}^*}} \widehat{\tilde{D_{ij}^*}} \right] } { \widehat{\bar{\rho}} {{k^\mathrm{test}}^{{3}\over {2}}}/{\hat{\Delta}} }

The other coefficient, C_{\mu_{\epsilon}}, is computed similarly to the (3.168) as

(3.186)C_{\mu_{\epsilon}} = \frac{ L_{u_i u_j} M_{u_i u_j} }{ M_{u_i u_j} M_{u_i u_j} },

where L_{u_i u_j} is defined identically to (3.166) and M_{u_i u_j} is simplified by

(3.187)M_{u_i u_j} = 2 \hat{\bar{\rho}} \hat{\Delta} {{k^\mathrm{test}}^{{1}\over {2}}} \widehat{\tilde{S}_{ij}^*} .

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)\mu_{t, {\rm clip}} = \operatorname{max}(\mu_{t, {\rm model}}, - \tau_{\rm clip} \bar{\mu})

with a value of \tau_{\rm clip} 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 \tau_{\rm clip} 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, G_B to both the turbulent kinetic energy and dissipation rate equation that is related to the baroclinic torque,

(3.189)G_B = {C_{bvg} (\mu+ \mu_t) \left\| {\partial {\rho} \over {\partial x_j}} \times {\partial {P} \over {\partial x_j} } \right\| \over \rho^2}.

of the model.

The buoyancy model of Rodi is given by

(3.190)G_B = \beta {\mu_t \over Pr_t} { \partial T \over \partial x_j } g_j.

De Ris’ buoyancy model offers two versions - flaming

(3.191)C_{deris} (\rho_\infty-\rho_f) \vec{g} k^{0.5}

and non-flaming.

(3.192)C_{deris} \Delta k^{0.5} (|\nabla\rho\times \vec{g}| - \nabla\rho\cdot \vec{g})

Default value of the user-defined coefficient C_{deris} 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 \rm{rhs} += \int G_B dV. For the dissipation rate equation, the source term is \rm{rhs} += \int C_{\epsilon3} {1 \over T} G_B dV for the buoyant vorticity generation model while it is \rm{rhs} += \int C_{\epsilon1}' C_{\epsilon4} {1 \over T} G_B dV otherwise. Recall that the inverse time scale is determined by the turbulence model of choice, i.e., \epsilon \over k for the standard k-\epsilon model and provided in (3.131) for the v^2-f model.

Note that the use of the buoyancy models has not been evaluated with the v^2-f 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 k - \epsilon turbulence models.

Turbulence Model

Symbol

User Input Name

Default Value

Standard k - \epsilon

C_{\mu}

Cmu

0.09

C_{\epsilon 1}

Ceps_1

1.44

C_{\epsilon 2}

Ceps_2

1.92

C_{\chi}

Cchi

2.0

\sigma_{k}

Sigma_K

1.0

\sigma_{\epsilon}

Sigma_E

1.3

Low Reynolds k - \epsilon

C_{\mu}

Cmu

0.09

C_{\epsilon 1}

Ceps_1

1.44

C_{\epsilon 2}

Ceps_2

1.92

\sigma_{k}

Sigma_K

1.0

\sigma_{\epsilon}

Sigma_E

1.3

A_{\mu}

Amu

3.4

RNG k - \epsilon

C_{\mu}

Cmu

0.0837

C_{\epsilon 1}

Ceps_1

1.42

C_{\epsilon 2}

Ceps_2

1.68

\sigma_{k}

Sigma_K

0.7194

\sigma_{\epsilon}

Sigma_E

0.7194

v^2 - f

C_{\mu}

Cmu

0.22

C_{\epsilon 1}

Ceps_1

1.4

C_{\epsilon 2}

Ceps_2

1.9

\sigma_{k}

Sigma_K

1.0

\sigma_{\epsilon}

Sigma_E

1.0

C_1

CF_1

0.4

C_2

CF_2

0.3

\alpha

Alpha

0.6

C_{T}

Nseg

6.0

C_{L}

CL

0.23

C_{\eta}

Ceta

70.0
Table 3.2 Constant parameters for k - \omega turbulence models.

Turbulence Model

Symbol

User Input Name

Default Value

k - \omega

\beta_{0}

Beta_Zero

0.0708

\beta^{*}

Beta_Star

0.09

\sigma_{k}

Sigma_K

3 / 5

\sigma_{\omega}

Sigma_W

0.5

\gamma

Gamma

13 / 25

C_{lim}

Clim

7 / 8

SST

A_{1}

A_One

0.31

\beta_{1}

Beta_One

0.075

\beta_{2}

Beta_Two

0.0828

\beta^{*}

Beta_Star

0.09

\gamma_{1}

Gamma_One

5 / 9

\gamma_{2}

Gamma_Two

0.44

\sigma_{k 1}

Sigma_K_One

0.85

\sigma_{k 2}

Sigma_K_Two

1.0

\sigma_{\omega 1}

Sigma_W_One

0.5

\sigma_{\omega 2}

Sigma_W_Two

0.856
Table 3.3 Constant parameters for LES turbulence models.

Turbulence Model

Symbol

User Input Name

Default Value

One-equation

C_{v}

Cv

0.5

C_{\epsilon}

Ceps

0.845

C_{\mu_{\epsilon}}

Cmueps

0.0856

Standard Smagorinsky

C_{v}

Cv

0.5

C_{s}

Cs

0.17

Dynamic Smagorinsky

C_{s}

Cs

0.17
Table 3.4 Constant parameters for miscellaneous turbulence models. Default values may be changed using the k - \epsilon model parameters input.

Model

Symbol

User Input Name

Default Value

Buoyant vorticity generation

C_{BVG}

Cbvg

0.35

C_{\epsilon 3}

Ceps_3

0.0

Rodi’s source term

C_{\epsilon 4}

C_eps4

0.0

EDC laminar limit

C_{\gamma,lam}

Cgammalam

2.0

C_{\tau,lam}

Ctaulam

0.02

C_{lam,trans}

Clamtrans

40.0