3.2.5. Turbulent Flow Equations, Favre-Averaged

The Favre-averaged turbulent transport equations are derived from the laminar equations of Laminar Flow Equations by passing the equations through either the RANS temporal filter of (3.68) or the LES spatial filter of (3.72). The mathematical form of the equations are essentially identical between the two filtering methods, so only a single set of equations will be presented. Care should be taken to interpret the filters as either temporal or spatial, depending on the closure models selected. While it is the Favre-averaged form of the equations that are solved, a comparison of the simple Reynolds-averaged and the Favre-averaged form is given at the end of this page for reference.

The approach most commonly used in turbulence modeling is called the Boussinesq eddy viscosity approximation, which relates the turbulent stress tensor to the filtered strain rate tensor through a modeled turbulent eddy viscosity. This general modeling approach has shown remarkable success for a broad range of problems (Wilcox [13]), and is the approach used in SIERRA/Fuego. A similar approach is used for scalar transport, where the scalar flux vector is related to scalar gradients through a modeled diffusion coefficient.

The following subsections describe the turbulent transport equations expressed in terms of a turbulent eddy viscosity or turbulent diffusion coefficient through the Boussinesq approximation. The treatment of these coefficients is dependent upon which of the many closure models are selected, and will be described in Turbulence Closure Models.

3.2.5.1. Conservation of Mass

The integral form of the Favre-filtered continuity equation used for turbulent transport is

(3.77)\int {{\partial \bar{\rho}} \over {\partial t}} {\rm d}V
  + \int \bar{\rho} \tilde{u}_j n_j {\rm d}S = 0.

This equation is in closed form, and no additional modeling is required.

3.2.5.2. Conservation of Momentum

The integral form of the Favre-filtered momentum equations used for turbulent transport are

(3.78)\int {{\partial \bar{\rho} \tilde{u}_i} \over {\partial t}} {\rm d}V
  +  \int \bar{\rho} \tilde{u}_i \tilde{u}_j n_j {\rm d}S
  +  \int \bar{p} n_i {\rm d}S  =
  \int \bar{\tau}_{ij} n_j {\rm d}S
  + \int \tau_{u_i u_j} n_j {\rm d}S
  + \int \left(\bar{\rho} - \rho_{\circ} \right) g_i {\rm d}V,

where the turbulent stress \tau_{u_i u_j} is defined as

(3.79)\tau_{u_i u_j} \equiv - \bar{\rho} ( \widetilde{u_i u_j} -
  \tilde{u}_i \tilde{u}_j ).

3.2.5.2.1. RANS Modeling

For RANS simulations, \tau_{u_i u_j} represents the Reynolds stress tensor and can be reduced to the form \tau_{u_i u_j} = -\overline{\rho u_i'' u_j''} by substitution of the Favre decomposition u_i \equiv \tilde{u}_i + u_i'' of each variable and simplifying. The deviatoric (trace-free) part of the stress tensor is defined as

(3.80)\tau_{u_i u_j}^D  \equiv  \tau_{u_i u_j} - \frac{1}{3} \tau_{u_k u_k} \delta_{ij}

(3.81)\tau_{u_i u_j}^D =  \tau_{u_i u_j} + \frac{2}{3} \bar{\rho}\tilde{k} \delta_{ij}

where the turbulent kinetic energy is defined as \tilde{k} \equiv \frac{1}{2} \widetilde{u_k'' u_k''}. The deviatoric part of the Reynolds stress tensor is modeled by the Boussinesq approximation which relates the Reynolds stresses to the filtered strain rate tensor through a modeled turbulent viscosity \mu_t, resulting in

(3.82)\tau_{u_i u_j}^D  =  \mu_t \left( \frac{\partial \tilde{u}_i}{\partial x_j}
  + {{\partial \tilde{u}_j} \over {\partial x_i}} \right) - {2 \over 3}
  \mu_t \frac{\partial \tilde{u}_k}{\partial x_k} \delta_{ij}

(3.83)\tau_{u_i u_j}^D  =  2 \mu_t \left( \tilde{S}_{ij} - \frac{1}{3} \tilde{S}_{kk}
  \delta_{ij} \right),

where the filtered strain rate tensor is defined by

(3.84)\tilde{S}_{ij} \equiv \frac{1}{2} \left( \frac{ \partial \tilde{u}_i }{
  \partial x_j } + \frac{ \partial \tilde{u}_j }{ \partial x_i } \right).

Substituting this into (3.81) yields the modeled form of the full Reynolds stress tensor (Kuo [12], p. 445)

(3.85)\tau_{u_i u_j} = 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}.

The Favre-filtered momentum equations then become

(3.86)\lefteqn{ \int {{\partial \bar{\rho} \tilde{u}_i} \over {\partial t}}
  {\rm d}V + \int \bar{\rho} \tilde{u}_i \tilde{u}_j n_j {\rm d}S
  + \int \left(\bar{p} + {2\over3}\bar{\rho}\tilde{k} \right) n_i
  {\rm d}S = }

(3.87)\int 2 (\mu + \mu_t) \left( \tilde{S}_{ij} - \frac{1}{3}
  \tilde{S}_{kk} \delta_{ij} \right) n_j {\rm d}S
  + \int \left(\bar{\rho} - \rho_{\circ} \right) g_i {\rm d}V,

where RANS closure models for the turbulent viscosity \mu_t are presented in Turbulence Closure Models.

3.2.5.2.2. LES Modeling

For LES, \tau_{u_i u_j} in (3.78) represents the subgrid stress tensor. The deviatoric part of the subgrid stress tensor is defined as

(3.88)\tau_{u_i u_j}^D \equiv  \tau_{u_i u_j} - \frac{1}{3} \tau_{u_k u_k}
  \delta_{ij}

(3.89)\tau_{u_i u_j}^D =  \tau_{u_i u_j} + \frac{2}{3} \bar{\rho} q^2 \delta_{ij},

where the subgrid turbulent kinetic energy is defined as q^2 \equiv \frac{1}{2} (\widetilde{u_k u_k} - \tilde{u_k} \tilde{u_k} ). The deviatoric part of the subgrid stress tensor is then modeled similar to RANS closures as (Moin, et al. [14])

(3.90)\tau_{u_i u_j}^D = 2 \mu_t \left( \tilde{S}_{ij} - \frac{1}{3} \tilde{S}_{kk}
  \delta_{ij} \right).

Substituting this into (3.89) yields the modeled form of the full subgrid stress tensor

(3.91)\tau_{u_i u_j} = 2 \mu_t \left( \tilde{S}_{ij} - \frac{1}{3} \tilde{S}_{kk}
  \delta_{ij} \right) - \frac{2}{3} \bar{\rho} q^2 \delta_{ij}.

For low Mach-number flows, a vast majority of the turbulent kinetic energy is contained at resolved scales (Erlebacher, et al. [15]). For this reason, the subgrid turbulent kinetic energy q^2 will not be directly treated and will instead be included in the pressure as an additional normal stress. The Favre-filtered momentum equations then become

(3.92)\lefteqn{ \int {{\partial \bar{\rho} \tilde{u}_i} \over {\partial t}}
  {\rm d}V + \int \bar{\rho} \tilde{u}_i \tilde{u}_j n_j {\rm d}S
  + \int \left( \bar{p} + \frac{2}{3} \bar{\rho} q^2 \right)
  n_i {\rm d}S = }

(3.93)\int 2 (\mu + \mu_t) \left( \tilde{S}_{ij} - \frac{1}{3}
  \tilde{S}_{kk} \delta_{ij} \right) n_j {\rm d}S
  + \int \left(\bar{\rho} - \rho_{\circ} \right) g_i {\rm d}V,

where LES closure models for the subgrid turbulent eddy viscosity \mu_t are presented in Turbulence Closure Models.

3.2.5.3. Conservation of Energy

The integral form of the Favre-filtered energy equation used for turbulent transport is

(3.94)\int {{\partial \bar{\rho} \tilde{h}} \over {\partial t}} {\rm d}V
  + \int \bar{\rho} \tilde{h} \tilde{u}_j n_j {\rm d}S
  = - \int \bar{q}_j n_j {\rm d}S
  - \int \tau_{h u_j} n_j {\rm d}S
  - \int {{\partial \bar{q}_i^r} \over {\partial x_i}} {\rm d}V

(3.95)+ \int \left({{\partial P} \over {\partial t}} + \tilde{u}_j {{\partial P} \over {\partial x_j}}\right){\rm d}V
  + \int \overline{\tau_{ij} {{\partial u_i} \over {\partial x_j }}} {\rm d}V.

The simple Fickian diffusion velocity approximation, (3.46), is assumed, so that the mean diffusive heat flux vector \bar{q}_j is

(3.96)\bar{q}_j = - \overline{ \left[ {\mu \over {\rm Pr}}
  {{\partial h} \over {\partial x_j}}
  -  {\mu \over {\rm Pr}}
  \sum_{k=1}^K h_k {{\partial Y_k} \over {\partial x_j}} \right] }
  - \overline{ {\mu \over {\rm Sc}}
  \sum_{k=1}^K h_k {{\partial Y_k} \over {\partial x_j}} }.

If Sc = Pr, i.e., unity Lewis number (Le = 1), then the diffusive heat flux vector simplifies to \bar{q}_j = -\frac{\mu}{\mathrm{Pr}} \frac{\partial \tilde{h}}{\partial x_j}. The viscous dissipation term is closed by

(3.97)\overline{\tau_{ij} {{\partial u_i} \over {\partial x_j }}}
  &= \left(\left(\mu + \mu_t\right) \left( {{\partial \tilde{u}_i}
  \over {\partial x_j}}
  + {{\partial \tilde{u}_j} \over {\partial x_i}} \right)
  - {2 \over 3} \left( \bar{\rho} \tilde{k} +
  \mu_t{{\partial \tilde{u}_k} \over {\partial x_k}} \right)
  \delta_{ij} \right) {{\partial \tilde{u}_i} \over {\partial x_j}}

  &= \left[ 2 \mu \tilde{S}_{ij}
  + 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 x_j}.

The turbulent diffusive flux vector \tau_{h u_j} in (3.95) is defined as

(3.98)\tau_{h u_j} \equiv \bar{\rho} \left( \widetilde{h u_j} -
  \tilde{h} \tilde{u}_j \right).

For RANS simulations, \tau_{h u_j} represents the turbulent energy diffusive flux vector and is simplified to the form \tau_{h u_j} = \overline{\rho h'' u_j''} by substitution of the Favre decomposition of each variable. It is then modeled by

(3.99)\tau_{h u_j} = \overline{\rho h'' u_j''} = - \frac{\mu_t}{\mathrm{Pr}_t} \frac{\partial \tilde{h}}{\partial x_j},

where \mathrm{Pr}_t is the turbulent Prandtl number and \mu_t is the modeled turbulent eddy viscosity from momentum closure. For LES, \tau_{h u_j} represents the subgrid turbulent energy diffusive flux vector, and is modeled in the same way as

(3.100)\tau_{h u_j} = - \frac{\mu_t}{\mathrm{Pr}_t} \frac{\partial \tilde{h}}{\partial x_j},

where \mathrm{Pr}_t is the subgrid turbulent Prandtl number and \mu_t is the modeled subgrid turbulent eddy viscosity from momentum closure.

The resulting filtered and modeled turbulent energy equation for both RANS and LES is given in Libby and Williams [11], p. 25, as

(3.101)\int {{\partial \bar{\rho} \tilde{h}} \over {\partial t}} {\rm d}V
  + \int \bar{\rho} \tilde{h} \tilde{u}_j n_j {\rm d}S
  = \int \left( {\mu \over {\rm Pr}} + {{\mu_t} \over {{\rm Pr}_t}} \right)
  {{\partial \tilde{h}} \over {\partial x_j}}  n_j {\rm d}S
  - \int {{\partial \bar{q}_i^r} \over {\partial x_i}} {\rm d}V
  + \int \left({{\partial P} \over {\partial t}} + \tilde{u}_j
  {{\partial P} \over {\partial x_j}}\right){\rm d}V
  + \int \overline{\tau_{ij} {{\partial u_j} \over {\partial x_j }}} {\rm d}V.

This equation is also given in Gran et al. [16] (without the transient and radiation source terms and the additional term for laminar transport). The turbulent Prandtl number must have the same value as the turbulent Schmidt number for species transport to maintain unity Lewis number.

3.2.5.4. Conservation of Species

The integral form of the Favre-filtered species equation used for turbulent transport is

(3.102)\int {{\partial \bar{\rho} \tilde{Y}_k} \over {\partial t}} {\rm d}V
  + \int \bar{\rho} \tilde{Y}_k \tilde{u}_j n_j {\rm d}S =
  - \int \tau_{Y_k u_j} n_j {\rm d}S
  - \int \overline{\rho Y_k \hat{u}_{j,k}} n_j {\rm d}S +
  \int \overline{\dot{\omega}_k} {\rm d}V,

where the form of diffusion velocities (see (3.46)) assumes the Fickian approximation with a constant value of diffusion velocity for consistency with the turbulent form of the energy equation, (3.95).

The turbulent diffusive flux vector \tau_{Y_k u_j} is defined as

(3.103)\tau_{Y_k u_j} \equiv \bar{\rho} \left( \widetilde{Y_k u_j} -
  \tilde{Y_k} \tilde{u}_j \right).

For RANS simulations, \tau_{Y_k u_j} represents the turbulent species diffusive flux vector and is simplified to the form \tau_{Y_k u_j} = \overline{\rho Y_k'' u_j''} by substitution of the Favre decomposition of each variable. It is then modeled as

(3.104)\tau_{Y_k u_j} = \overline{\rho Y_k'' u_i''}
  = - \frac{\mu_t}{\mathrm{Sc}_t} \frac{\partial \tilde{Y}_k}{\partial x_i},

where \mathrm{Sc}_t is the turbulent Schmidt number for all species and \mu_t is the modeled turbulent eddy viscosity from momentum closure. For LES, \tau_{Y_k u_j} represents the subgrid turbulent species diffusive flux vector, and is modeled identically as

(3.105)\tau_{Y_k u_j} = - \frac{\mu_t}{\mathrm{Sc}_t} \frac{\partial
  \tilde{Y}_k}{\partial x_i},

where \mathrm{Sc}_t is the subgrid turbulent Schmidt number for all species and \mu_t is the subgrid modeled turbulent eddy viscosity from momentum closure.

The Favre-filtered and modeled turbulent species transport equation for both RANS and LES then becomes (Gran et al. [16])

(3.106)\int {{\partial \bar{\rho} \tilde{Y}_k} \over {\partial t}} {\rm d}V
  + \int \bar{\rho} \tilde{Y}_k \tilde{u}_j n_j {\rm d}S =
  \int \left( {\mu \over {\rm Sc}}
  + {{\mu_t} \over {{\rm Sc}_t}}  \right)
  {{\partial \tilde{Y}_k} \over
  {\partial x_j}} n_j {\rm d}S +
  \int \overline{\dot{\omega}}_k {\rm d}V .

If transporting both energy and species equations, the laminar Prandtl number must be equal to the laminar Schmidt number and the turbulent Prandtl number must be equal to the turbulent Schmidt number to maintain unity Lewis number. Although there is a species conservation equation for each species in a mixture of K species, only K-1 species equations need to be solved since the mass fractions sum to unity and

(3.107)\tilde{Y}_k = 1 - \sum_{j \ne k}^{K} \tilde{Y}_j .

3.2.5.5. Radiation Transport

The Favre-averaged energy equation, (3.101), requires the time-averaged radiative flux divergence. From (3.63), the time-averaged radiative flux divergence is given by

(3.108)\overline{{{\partial q_i^r} \over {\partial x_i}}} =
  4 \sigma \overline{\mu_a T^4} - \overline{\mu_a G} .

For optically thin turbulent eddies, which is the case for many combustion applications, fluctuations in the absorption coefficient and the scalar flux are weakly correlated [7] so (3.108) may be simplified to

(3.109)\overline{{{\partial q_i^r} \over {\partial x_i}}} =
  4 \sigma \overline{\mu_a T^4}
  - \bar{\mu}_a \bar{G} .

The time averaged scalar flux is obtained from the time averaged Boltzmann radiative transport equation

(3.110)s_i {{\partial} \over {\partial x_i}} \overline{I \left(s\right)}
  + \bar{\mu}_a \overline{I\left(s\right)} =
  {\overline{\mu_a \sigma T^4} \over {\pi}} ,

where the correlation between the turbulent fluctuations in the absorption coefficient and the intensity is assumed small to simplify the absorption term.

Both (3.109) and (3.110) include the time averaged emission term, \overline{\alpha T^4}, which may significantly increase the radiative emission from a turbulent flame above what would be estimated from the mean temperature and absorption coefficient values. The details of the closure used for this term are discussed in the turbulent combustion model section.

3.2.5.6. Comparison Between Time-Averaging and Favre-Averaging

The time-averaged (Reynolds-averaged) and Favre-averaged transport equations are given in the following section.

3.2.5.6.1. Conservation of Mass

The continuity equation:

  • Time averaged:

    (3.111)\int {{\partial \overline{\rho}} \over {\partial t}} {\rm d}V
  + \int \overline{\rho} \overline{u}_j n_j {\rm d}S
  + \int \overline{\rho' u'_j} n_j {\rm d}S = 0

  • Favre averaged:

    (3.112)\int {{\partial \overline{\rho}} \over {\partial t}} {\rm d}V
  + \int \overline{\rho} \tilde{u}_j n_j {\rm d}S = 0

3.2.5.6.2. Conservation of Momentum

The momentum transport equations:

  • Time averaged:

    (3.113)\int {{\partial \overline{\rho} \overline{u}_i} \over {\partial t}} {\rm d}V
  + \int {{\partial \overline{\rho' u'_i}} \over {\partial t}} {\rm d}V
  +  \int \overline{\rho}\ \overline{u}_i \overline{u}_j n_j {\rm d}S
  +  \int \overline{\rho} \overline{u'_i u'_j} n_j {\rm d}S
  +  \int \overline{u_j} \overline{\rho' u'_i} n_j {\rm d}S
  +   \int \overline{u_i} \overline{\rho' u'_j} n_j {\rm d}S
  +  \int \overline{\rho' u'_i u'_j} n_j {\rm d}S
  +  \int \overline{p} n_i {\rm d}S
  =  \int \overline{\tau}_{ij} n_j {\rm d}S
  + \int \overline{\rho} g_i {\rm d}V

  • Favre averaged:

    (3.114)\int {{\partial \overline{\rho} \tilde{u}_i} \over {\partial t}} {\rm d}V
  +  \int \overline{\rho} \tilde{u}_i \tilde{u}_j n_j {\rm d}S
  +  \int \overline{p} n_i {\rm d}S  =
  \int \overline{\tau}_{ij} n_j {\rm d}S
  - \int \overline{\rho u_i'' u_j''} n_j {\rm d}S
  + \int \overline{\rho} g_i {\rm d}V

3.2.5.6.3. Conservation of Energy

The energy transport equation, assume Lewis number is one:

  • Favre averaged:

    (3.115)\int {{\partial \overline{\rho} \tilde{h}} \over {\partial t}} {\rm d}V
  + \int \overline{\rho} \tilde{h} \tilde{u}_j n_j {\rm d}S
  =  \int \overline{{\kappa \over {C_p}}
  {{\partial h} \over {\partial x_j}}}  n_j {\rm d}S
  - \int \overline{\rho h'' u''_j} n_j {\rm d}S
  -  \int {{\partial \overline{q_i^r}} \over {\partial x_i}} {\rm d}V

3.2.5.6.4. Conservation of Species

The species transport equation:

  • Time averaged:

    (3.116)\int {{\partial \overline{\rho} \overline{Y}_k} \over {\partial t}} {\rm d}V
  +  \int {{\partial \overline{\rho' Y'_k}} \over {\partial t}} {\rm d}V
  + \int \overline{\rho} \overline{Y}_k \overline{u}_j n_j {\rm d}S
  + \int \overline{\rho} \overline{Y'_k u'_j} n_j {\rm d}S
  + \int \overline{u_j} \overline{\rho' Y'_k} n_j {\rm d}S
  +  \int \overline{Y_k} \overline{\rho' u'_j} n_j {\rm d}S
  + \int \overline{\rho' Y'_k u'_j} n_j {\rm d}S
  =  \int \overline{\rho Y_k \hat{u}_{j,k}} n_j {\rm d}S +
  \int \overline{\dot{\omega}_k} {\rm d}V

  • Favre averaged:

    (3.117)\int {{\partial \overline{\rho} \tilde{Y}_k} \over {\partial t}} {\rm d}V
  + \int \overline{\rho} \tilde{Y}_k \tilde{u}_j n_j {\rm d}S =
  - \int \overline{\rho Y_k'' u_j''} n_j {\rm d}S
  + \int \overline{\rho Y_k \hat{u}_{j,k}} n_j {\rm d}S +
  \int \overline{\dot{\omega}_k} {\rm d}V