3.2.11. Laminar Flamelet Turbulent Combustion Model

Laminar flamelet models for non-premixed turbulent combustion treat turbulent flames as an ensemble of laminar diffusion flames. [53] Nonequilibrium chemistry effects may be included in the model by accounting for localized fluid strain, resulting in what is classically called the Strained Laminar flamelet Model (SLFM). Nonadiabatic effects may also be included by accounting for losses to the surroundings in the ensemble of flamelets.

The fundamental assumption is that the chemical time scales of the important reactions are fast enough that they occur only in a thin layer around stoichiometry, thinner (ideally) than the smallest scales of the turbulence. Defining a small quantity $\epsilon = \ell_{\rm reaction \, zone} / \ell_{\rm mixing \, layer} \ll 1$ , we can examine the governing equations in that thin region using a multiscale asymptotic expansion as

(3.339) $Y_i = Y_i(\zeta, \tau, x, t)+\epsilon Y_i^1(\zeta,\tau, x, t) + \ldots, \qquad \zeta = \frac{Z(x,t)-Z_{\rm st}}{\epsilon} \mbox { and } \tau = t / \epsilon^2.$

Collecting the dominant terms, making some standard simplifications, and assuming that the chemical reaction scales as $\epsilon^{-2}$ , the state of the gas depends on the flow scale $Z$ and $\chi = 2 D |\nabla Z|^2$ :

(3.340) $\rho \pd{Y_i}{t} - \frac{\rho \chi}{2} \frac{1}{{\rm Le}_i}\pd{^2 Y_i}{Z^2} - \dot{\omega_i}(\vec{\Phi}) = 0$

(3.341) $\rho \pd{T}{t} - \frac{\rho \chi}{2} \left( \pd{^2 T}{Z^2} + \frac{1}{c_p} \pd{c_p}{Z} \pd{T}{Z}\right) - \dot{\omega}_T(\vec{\Phi}) = 0$

(3.342) $T(Z=0,t) = T_{\rm ox.} \mbox{, } T(Z=1,t) = T_{\rm fuel} \mbox{, } Y_i(Z=0,t) = Y_{\rm i, ox.} \mbox{, } Y_i(Z=1,t) = Y_{\rm i, fuel},$

(3.343) $\rho = \rho(\vec{\Phi}) \mbox{ and } c_p = c_p(\vec{\Phi}),$

where $\vec{\Phi}$ is the state vector $\vec{\Phi} = (P_{\rm{th}}, T, Y_0, Y_1, \ldots, Y_N)$ . The approximation allows us to resolve the chemical scales in the phase space of the mixture fraction instead of on a three-dimensional grid, granting dramatic computational savings. If we make the additional assumption that the chemistry is quasi-steady on the scale of the flow, then the chemical structure in mixture fraction space can be pre-computed offline from the simulation for a range of flow parameters $\chi$ and tabulated (using fuego_tabular_props). During the flow simulation, the solution of the flamelet simulation can be queried to determine required flow properties, e.g. $\rho = \rho(Z,\chi)$ . Note that the flamelet formulation in Eq. Laminar Flamelet Turbulent Combustion Model is specifically for a “two-stream” problem, with constant Lewis numbers, where the boundary and initial conditions of the simulation can be completely described by linear combinations of two constant state vectors. Additional “streams” and boundary heat losses will require additional transport equations to be solved.

This section summarizes the basic formulation and implementation details of both the adiabatic and nonadiabatic flamelet model and SLF model, including both the property table generation procedure in fuego_tabular_props and the usage of the property table in fuego to evaluate turbulent filtered quantities of interest for both adiabatic and nonadiabatic configurations.

3.2.11.1. Adiabatic Property Table Generation

3.2.11.1.1. Laminar Flamelet Generation

Unstrained flamelet libraries, where nonequilibrium chemistry effects may be neglected with respect to fluid strain rates, can be generated directly with the fuego_tabular_props application. These libraries should be used either in laminar flow or in turbulent flow where the turbulence/chemistry interactions may be neglected.

Equilibrium chemistry, Burke-Schumann chemistry, or nonreacting flow scenarios are supported in configurations where there are two or more streams that may be mixed and potentially reacted. The stream composition is parameterized by the mixture fraction vector $Z_m$ , where each of the $M$ component represents the fraction of mass that originated in that stream, where there are $N$ streams and $M=N-1$ . The mixture fraction for the final stream may be evaluated as $Z_N = 1 - \sum_{m=1}^M Z_m$ .

The resulting flamelet data can then be assembled into a sequence of multi-dimensional tables of dependent variable $\phi$ as a function of the mixture fraction vector, $\phi(Z_m)$ , and can be used directly for laminar simulations. Adding turbulence interactions, nonequilibrium effects, and nonadiabatic effects will increase the dimensionality of this lookup table and require additional processing. See the following sections for more information.

3.2.11.1.2. Strained Laminar Flamelet Libraries

Strained laminar flamelet data may be generated for use in Fuego with the Spitfire code. This data is two-dimensional in nature, determined by the mixture fraction and a reference scalar dissipation rate $\chi_o$ at a reference mixture fraction $Z_o$ . The instantaneous laminar scalar dissipation rate is defined as

(3.344) $\chi = 2 D \frac{\partial Z}{\partial x_i} \frac{\partial Z}{\partial x_i},$

with D being the molecular mass diffusion coefficient. The reference value $\chi_o$ is arbitrary, although typical choices include the stoichiometric value $\chi_\mathrm{st} = \chi(Z_\mathrm{st})$ or the maximum value $\chi_\mathrm{max}=\chi(Z=0.5)$ . Stoichiometric values are used in Fuego. Spitfire will then assemble it into a sequence of multi-dimensional tables of dependent variable $\phi$ as a function of the mixture fraction vector and reference scalar dissipation rate, $\phi(Z_m, \chi_o)$ . These tables may be then be used to evaluate properties in Fuego.

3.2.11.1.3. Aksit-Moss SNL soot model

A modified Aksit-Moss soot model [54] is available, solving for soot moles per mass ( $N_s$ ) and soot mass fraction ( $Y_s$ )

(3.345) $\int \frac{\partial \bar{\rho} \widetilde{Z}}{\partial t} \, {\rm d} V + \int \bar{\rho} \widetilde{u}_j \widetilde{Z} \, n_j {\rm d} S = \int \left(\frac{\bar{\mu}}{\rm \, Sc} \frac{\partial \widetilde{Z}}{\partial x_j} - \tau^{sgs}_{Z,j} \right) \, n_j {\rm d} S + \int\dot{S}_Z \, {\rm d} V$

(3.346) $\int \frac{\partial \bar{\rho} \widetilde{h}}{\partial t} \, {\rm d} V + \int \bar{\rho} \widetilde{u}_j \widetilde{h} \, n_j {\rm d} S = \int \left(\frac{\bar{\mu}}{\rm \, Pr}\frac{\partial \widetilde{h}}{\partial x_j} - \tau^{sgs}_{h,j} \right) \, n_j {\rm d} S - \int \dot{S}_h \, {\rm d} V$

(3.347) $\int \frac{\partial \bar{\rho} \widetilde{h}_c}{\partial t} \, {\rm d} V + \int \bar{\rho} \widetilde{u}_j \widetilde{h_c} \, n_j {\rm d} S = \int \left(\frac{\bar{\mu}}{\rm \, Pr} \frac{\partial \widetilde{h}_c}{\partial x_j} - \tau^{sgs}_{h_c,j} \right) \, n_j {\rm d} S$

(3.348) $\int \frac{\partial \bar{\rho} \widetilde{N_s}}{\partial t} \, {\rm d} V + \int \bar{\rho} \widetilde{u}_j \widetilde{N}_s \, n_j {\rm d} S = \int \left(\frac{\bar{\mu}}{{\rm \, Sc}_s} \frac{\partial \widetilde{N_s}}{\partial x_j} - \tau^{sgs}_{N_s,j} \right) \, n_j {\rm d} S + \int \dot{S}_{N_s} \, {\rm d} V$

(3.349) $\int \frac{\partial \bar{\rho} \widetilde{Y_s}}{\partial t} \, {\rm d} V + \int \bar{\rho} \widetilde{u}_j \widetilde{Y}_s \, n_j {\rm d} S = \int \left(\frac{\bar{\mu}}{{\rm \, Sc}_s} \frac{\partial \widetilde{Y_s}}{\partial x_j} - \tau^{sgs}_{Y_s,j} \right) \, n_j {\rm d} S + \int \dot{S}_{Y_s} \, {\rm d} V$

with $W_s$ denoting the molecular weight of soot and $\gamma_L$ a leaning coefficient representing the rate at which mixture fraction is leaned per soot mass fraction rate. Source terms are defined as

(3.350) $\dot{S}_{N_s} = \dot{S}_{\rm nucleation} - \dot{S}_{\rm coagulation} \widetilde{Y_s}^{1/6} \widetilde{N}_s^{11/6}$

(3.351) $\dot{S}_{Y_s} = W_s \dot{S}_{\rm nucleation} + \left[\dot{S}_{\rm surface \; growth} - \dot{S}_{\rm O_2 \; oxid.} - \dot{S}_{\rm OH \; oxid.} \right] \widetilde{Y_s}^{2/3} \widetilde{N}_s^{1/3}$

(3.352) $\dot{S}_Z = -\gamma_L \dot{S}_{Y_s}$

(3.353) $\dot{S}_h = 4 \sigma \left( \overline{q^r_{\rm s}} Y_s + \overline{q^r_{\rm gas}} \right) - \left(\overline{\alpha_{soot}} Y_s + \overline{\alpha_{\rm gas}} \right) \overline{G}.$

Each of the rates, such as $\dot{S}_{\rm nucleation}$ , must be provided as a tabulated source term to the model, e.g. from activating the “Aksit-Moss SNL” soot model in SpitFire.

3.2.11.1.4. Turbulent Averaging

In turbulent simulations, a filtered form of the governing equations are solved to reduce the resolution requirements to an affordable level. Temporal filtering is used in Reynolds Averaged Navier-Stokes (RANS) models and spatial filtering is used in Large Eddy Simulation (LES) models. Both types of filtering are represented with the notation $\bar{\phi}$ , and are handled similarly in the present work. Density-weighted, or Favre filtering greatly simplifies the treatment of variable-density flow. A Favre-filtered quantity is represented by $\tilde{\phi} \equiv \overline{\rho \phi} / \bar{\rho}$ .

For use in turbulent simulations, a Favre-filtered version of the variables in the property table must be calculated. This is performed by convoluting the property variable with the joint PDF of the independent variable sub-filter fluctuations, and is mathematically expressed as

(3.354) $\tilde{\phi}(\tilde{Z}_m,\widetilde{Z''^2},\tilde{\chi}) = \int\limits_0^\infty \int\limits_0^1 \phi(Z_m, \chi_o) \; p_{Z \chi}(Z_i,\chi;\tilde{Z}_m,\widetilde{Z''^2},\tilde{\chi}) \, \mathrm{d}Z_i\, \mathrm{d} \chi,$

where $p_{Z \chi}(Z_i,\chi;\tilde{Z}_m,\widetilde{Z''^2}, \tilde{\chi})$ is the joint PDF of sub-filter fluctuations of the dependent variable $\phi$ in $Z_i$ - $\chi$ space, parameterized by the filtered mixture fractions $\tilde{Z}_m$ and the variance $\widetilde{Z''^2}$ of mixture fraction component $Z_i$ , and the filtered scalar dissipation rate $\tilde{\chi}$ . The reference scalar dissipation rate has the functionality $\chi_o(\tilde{Z}_i,\widetilde{Z''^2},\tilde{\chi})$ , which will be discussed in the following section. Variance of only a single component of mixture fraction, $Z_i$ , is considered at present for simplicity, although extensions to include additional components are possible. Statistical independence will be assumed between $Z_i$ and $\chi$ fluctuations, so that

(3.355) $\tilde{\phi}(\tilde{Z}_m,\widetilde{Z''^2},\tilde{\chi}) = \int\limits_0^\infty \int\limits_0^1 \phi(Z_m, \chi_o) \; p_Z(Z_i;\tilde{Z}_m,\widetilde{Z''^2}) \, p_\chi(\chi;\tilde{\chi}) \, \mathrm{d}Z_i\, \mathrm{d} \chi.$

For the present work, $p_Z(Z_i;\tilde{Z}_m,\widetilde{Z''^2})$ will be modeled as either a beta PDF or a clipped Gaussian PDF and $p_\chi(\chi;\tilde{\chi})$ will be modeled as the delta function $\delta(\chi - \tilde{\chi})$ .

3.2.11.1.5. Property Table Implementation

The convolution integral in (3.355) would be prohibitively expensive to evaluate each time a value for $\tilde{\phi}$ is needed by a turbulent reacting simulation. Therefore, this integral will be pre-calculated so that each property table query will only involve an interpolation from a table of values.

Storing the final $\tilde{\phi}(\tilde{Z}_m,\widetilde{Z''^2}, \tilde{\chi})$ values directly is undesirable since the range of possible $\tilde{\chi}$ values for each flamelet is different, resulting in a non-orthogonal table. Instead, the values $\tilde{\phi}_T(\tilde{Z}_m,\widetilde{Z''^2}, \chi_o)$ are stored in an orthogonal table that is indexed by $\tilde{Z}_m$ , $\widetilde{Z''^2}$ , and $\chi_o(\tilde{Z}_i,\widetilde{Z''^2},\tilde{\chi})$ . These tabulated values are calculated by

(3.356) $\tilde{\phi}_T(\tilde{Z}_m,\widetilde{Z''^2},\chi_o) = \int_0^1 \phi(Z_m, \chi_o) \; p_Z(Z_i;\tilde{Z}_m, \widetilde{Z''^2}) \, \mathrm{d}Z_i.$

The reference scalar dissipation rate $\chi_o$ needed for lookup in the table for $\tilde{\phi}_T(\tilde{Z}_m,\widetilde{Z''^2}, \chi_o)$ can be evaluated from the local filtered scalar dissipation rate $\tilde{\chi}$ through laminar flamelet theory. The instantaneous scalar dissipation rate $\chi$ can be approximated by

(3.357) $\chi &= \chi_\mathrm{max} \exp\left( -2 [\mathrm{erfc}^{-1}(2Z)]^2 \right) &= \chi_\mathrm{max} F_\chi(Z),$

where $\chi_\mathrm{max}$ is the maximum scalar dissipation rate found in the counterflow diffusion flame, which occurs at the stagnation point where $Z = 0.5$ . (Note that this expression has not yet been extended to multiple mixture fractions, so that this treatment is only applicable for two-stream problems.) The value of $\chi$ at any reference location in the flamelet can be similarly approximated, so that $\chi_o = \chi_\mathrm{max} F_\chi(Z_o)$ . Combining these models by equating the unknown $\chi_\mathrm{max}$ yields a closed-form expression linking the scalar dissipation rate at any location to the reference value on the flamelet with the same characteristic $\chi_\mathrm{max}$ ,

(3.358) $\chi = \chi_o \frac{ F_\chi(Z) }{ F_\chi(Z_o) }.$

Applying the filtering operation in (3.355) to both sides of (3.358) for a single-mixture fraction configuration yields

(3.359) $\tilde{\chi} & = \int\limits_0^\infty \int\limits_0^1 \chi_o \frac{ F_\chi(Z) }{ F_\chi(Z_o) } \; p_Z(Z;\tilde{Z},\widetilde{Z''^2}) \, p_\chi(\chi;\tilde{\chi}) \, \mathrm{d}Z\, \mathrm{d} \chi & = \frac{\chi_o}{F_\chi(Z_o)} \int_0^\infty p_\chi(\chi;\tilde{\chi}) \, \mathrm{d}\chi \int_0^1 F_\chi(Z) \, p_Z(Z;\tilde{Z},\widetilde{Z''^2}) \, \mathrm{d}Z & = \frac{ \chi_o }{ F_\chi(Z_o) } \int_0^1 F_\chi(Z) \, p_Z(Z;\tilde{Z},\widetilde{Z''^2}) \, \mathrm{d}Z,$

so that the filtered reference scalar dissipation rate can be calculated from the filtered quantities provided by the turbulent flame simulation as

(3.360) $\chi_o(\tilde{Z},\widetilde{Z''^2},\tilde{\chi}) = \frac{ \tilde{\chi} \, F_\chi(Z_o) }{ \int_0^1 F_\chi(Z) \, p_Z(Z;\tilde{Z},\widetilde{Z''^2}) \, \mathrm{d}Z, }.$

To decrease computational cost, the integral in the denominator can be interpolated from pre-calculated values in a two-dimensional table as a function of $\tilde{Z}$ and $\tilde{Z''^2}$ .

To summarize, the turbulent reacting simulation will query the property table for the variable $\tilde{\phi}(\tilde{Z}_m,\widetilde{Z''^2},\tilde{\chi})$ . Internally, (3.360) will be used to calculate $\chi_o$ as a function of the provided filtered independent variables. This value will then be used along with the provided independent variables to interpolate a value for $\tilde{\phi}_T(\tilde{Z}_m,\widetilde{Z''^2}, \chi_o)$ from the stored table that was pre-calculated with (3.356). This interpolated value will then be returned to the main simulation as the requested value for $\tilde{\phi}(\tilde{Z}_m,\widetilde{Z''^2},\tilde{\chi})$ .

If turbulence/chemistry interactions are to be neglected in the simulation, the delta function $\delta(Z - \tilde{Z})$ may be used for $p_Z(Z;\tilde{Z})$ in (3.360) so that the reference scalar dissipation rate can be computed simply as

(3.361) $\chi_o(\tilde{Z},\tilde{\chi}) = \frac{ \tilde{\chi}\, F(Z_o) }{ F(\tilde{Z}) }.$

Once the multidimensional property table has been generated, it can be imported into fuego and queried for the dependent variables as a function of the independent variables $\tilde{Z}_m$ , $\widetilde{Z''^2}$ , and $\tilde{\chi}$ . Models are required for each of these independent variables used by the flamelet property table. Filtered Scalar Dissipation Rate to Filtered Scalar Variance present models for each of these quantities for each of the supported turbulence closure models.

3.2.11.2. Nonadiabatic Property Table Generation

When including the effects of radiative or convective heat losses in a flamelet simulation, additional parameterizations beyond those in the previous section are required. These are the “conserved enthalpy”, $h^*$ and heat loss parameter $\gamma$ , where the heat loss parameter is defined as $\gamma = h - h^*$ . The conserved enthalpy is identical to the traditional enthalpy except that its transport equation omits all source terms (typically due to radiative losses).

This formulation is used as a way to parameterize losses in a manner that is consistent with the opposed diffusion flame burner simulations used to generate the flamelet libraries. In these burner simulations, the inflowing pure stream states are fixed and cannot experience any heat losses; Losses only occur in the interior of the burner, and are represented by $\gamma$ variation. A range of inflowing pure stream states may also be computed, and are parameterized through $h^*$ variation. In this way, the full range of possible states may be tabulated and retrieved in a fire simulation through values of $h$ and $h^*$ , which are both straightforward to compute.

For turbulent simulations, the Favre-filtered property variable $\tilde{\phi}$ is evaluated as

(3.362) $\tilde{\phi}(\tilde{Z}_m,\widetilde{Z''^2},\tilde{\chi}, \tilde{\gamma}, \tilde{h}^*) = \int\limits_{-\infty}^{\infty} \int\limits_{-\infty}^{\infty} \int\limits_0^\infty \int\limits_0^1 \phi(Z_m, \chi_o, \gamma_o, h_o^*) \; p_{Z \chi \gamma h^*}(Z_m,\chi,\gamma,h^*;$

(3.363) $\tilde{Z}_m, \widetilde{Z''^2},\tilde{\chi},\tilde{\gamma}, \tilde{h}^*) \, \mathrm{d}Z_m\, \mathrm{d} \chi\, \mathrm{d} \gamma\, \mathrm{d} h^*,$

where $\gamma_o$ and $h_o^*$ are reference values of the heat loss parameter and the conserved enthalpy, respectively, to be defined in the following sections. Statistical independence will be assumed between fluctuations of each $Z_m$ component, $\chi$ , $\gamma$ , and $h^*$ , so that

(3.364) $\tilde{\phi}(\tilde{Z}_m,\widetilde{Z''^2},\tilde{\chi},\tilde{\gamma}, \tilde{h}^*) = \int\limits_{-\infty}^{\infty} \int\limits_{-\infty}^{\infty} \int\limits_0^\infty \int\limits_0^1 \int\limits_0^1 \phi(Z_m, \chi_o, \gamma_o, h_o^*) \; p_{Z_i}(Z_i;\tilde{Z}_i,\widetilde{Z''^2}) \, p_{Z_m}(Z_m;\tilde{Z}_m)$

(3.365) $p_\chi(\chi;\tilde{\chi}) \, p_\gamma(\gamma;\tilde{\gamma}) \, p_{h^*}(h^*;\tilde{h}^*) \, \mathrm{d}Z_i \, \mathrm{d}Z_{m \neq i} \, \mathrm{d}\chi \, \mathrm{d}\gamma \, \mathrm{d}h^*.$

For the present work, $p_{Z_i}(Z_i;\tilde{Z}_i,\widetilde{Z''^2})$ will be modeled as either a beta PDF or a clipped Gaussian PDF, and $p_{Z_m}(Z_m;\tilde{Z}_m)$ , $p_\chi(\chi;\tilde{\chi})$ , $p_\gamma(\gamma;\tilde{\gamma})$ , and $p_{h^*}(h^*;\tilde{h}^*)$ will be modeled as the delta functions $\delta(Z_m - \tilde{Z}_m)$ , $\delta(\gamma - \tilde{\gamma})$ , and $\delta(\chi - \tilde{\chi})$ , $\delta(\gamma - \tilde{\gamma})$ , and $\delta(h^* - \tilde{h}^*)$ , respectively.

The convolution integral in (3.365) would be prohibitively expensive to evaluate each time a value for $\tilde{\phi}$ is needed by a turbulent reacting simulation. Therefore, this integral will be pre-calculated so that each property table query will only involve an interpolation from a table of values.

Storing the final $\tilde{\phi}(\tilde{Z}_m,\widetilde{Z''^2}, \tilde{\chi}, \tilde{\gamma}, \tilde{h^*})$ values directly is undesirable since the range of possible $\tilde{\chi}$ , $\tilde{\gamma}$ , and $\tilde{h}^*$ values for each flamelet is different, resulting in a non-orthogonal table. Instead, the values $\tilde{\phi}_T(\tilde{Z}_m,\widetilde{Z''^2}, \chi_o, \gamma_o, h_o^*)$ are stored in an orthogonal table that is indexed by $\tilde{Z}_m$ , $\widetilde{Z''^2}$ , $\chi_o(\tilde{Z}_i,\widetilde{Z''^2},\tilde{\chi})$ , $\gamma_o(\tilde{Z}_i,\widetilde{Z''^2},\tilde{\gamma})$ , and $h_o^*(\tilde{Z}_i,\tilde{h}^*)$ . These tabulated values are calculated by

(3.366) $\tilde{\phi}_T(\tilde{Z}_m,\widetilde{Z''^2},\chi_o,\gamma_o,h_o^*) = \int\limits_0^1 \int\limits_0^1 \phi(Z_m, \chi_o, \gamma_o, h_o^*) \; p_{Z_i}(Z_i;\tilde{Z}_i, \widetilde{Z''^2}) \, p_{Z_m}(Z_m;\tilde{Z}_m) \, \mathrm{d}Z_i \, \mathrm{d}Z_{m \neq i}.$

The required reference values of $\gamma_o$ and $h_o^*$ are described in the following sections.

3.2.11.2.1. Boundary heat loss

The addition of a temperature boundary condition on the wall requires a modification of the flamelet formulation of Eq. Laminar Flamelet Turbulent Combustion Model. The equation for a normalized temperature variable, $\theta = T/T_{\rm ox.} - 1$ , is

(3.367) $\partial_t \left(\rho \theta \right) + \nabla \cdot \left( \rho \mathbf{u} \theta \right) - \nabla \cdot \left(\lambda \nabla \theta \right) = L_D$

(3.368) $\theta = \dot{\omega}_T(\theta, \vec{Y}) \mbox{ in } \Omega$

(3.369) $\theta(x \in \partial \Omega_{\rm fuel}) = \theta_{\rm fuel},$

(3.370) $\theta(x \in \partial \Omega_{\rm ox}) = 0,$

(3.371) $\theta(x \in \partial \Omega_{\rm wall}) = \theta_{\rm wall},$

which now has an extra boundary term $\theta(x \in \partial \Omega_{\rm wall}) = \theta_{\rm wall}$ . The extra boundary condition remains when we apply the flamelet approximation, leaving $\theta(Z=Z_{\rm wall}) = \theta_{\rm wall}$ . The value of the mixture fraction at the wall, however, is undecided: we only know that $\nabla Z \cdot n = 0$ . During the simulation, the value of mixture fraction directly evaluated at the wall can be determined dynamically and the value of temperature can be computed. Away from the wall, however, one in principle would need to follow the $\zeta$ coordinate from the flamelet transformation until it intersects the wall. However, given that $\nabla Z \cdot n = 0$ , the gradient trajectory in principle is tangential to the wall. Although the flamelet equation itself is well-posed, the asymptotic derivation of the flamelet model in the very near region to a nonadiabatic wall and the equations need to modified in some fashion to account.

Flamelets can readily be described when they are adiabatic; in the limit of unity Lewis numbers and adiabatic systems the enthalpy is a linear function of the mixture fraction. The existence of radiative transport and wall heat transfer introduces deviations from this linear relationship between $h$ and $Z$ . Heat losses at the predominant boundary temperature are a common scenario. Defining a reference boundary temperature at $T_{\rm fr}(Z) = (1-Z) T_{\rm ox} + Z T_{\rm fuel}$ , then this case a simplified flamelet temperature equation with heat losses could be written as

(3.372) $\rho \partial_t T - \frac{1}{2}\rho \chi \partial_Z^2 T = \dot{\omega}_T(\vec{\Phi}) - H_T (T - T_{\rm fr})$

where $H_T$ represents a heat transfer coefficient that will be further discussed below. This gives a heat loss term that is linear in $T$ . Alternately, the heat loss can be written specifically for radiative-style losses, $\dot{q}_{\rm losses} = \sigma(\vec{\Phi}) (T^4-T_{\rm fr}^4)$ . Regardless, with heat loss expressed in terms of $T_{\rm fr}$ the flamelet enthalpy is no longer linear in $Z$ but instead takes on a roughly inverted triangular form with an extrema at the peak temperature, roughly $Z_{st}$ . This has led us to express the difference between the adiabatic enthalpy, defined as $h_c = h_{\rm ox} (1-Z) + h_{\rm fuel} Z$ , and the actual flamelet computed enthalpy, $h$ , as $\gamma = h - h_c$ . The introduction of $\gamma$ is done strictly as an expedient for generation of flamelet libraries. By assuming a triangular form (or any particular assumed form) we can stretch the table entries into a square format by tabulating as a function of the stoichiometric value of $\gamma_{st}$ . This does require the use of the assumed form for $\gamma$ for converting from the local $Z$ value of $\gamma = h(Z) - h_c(Z)$ to $\gamma_{\rm st}$ . Comparison with DNS and unsteady flamelets for laminar flames shows good agreement with this type of enthalpy defect model for radiation in unity Lewis number flames (which is an appropriate assumption for turbulent, hydrocarbon fires) [55]. We make an assumption of path independence for the solution of at particular integrated heat loss, but in reality the solution will depend somewhat on the value of $H_T$ and the form of the added heat loss term.

The compensation for boundary heat loss can be extended to a full range of temperature (below $T_{\rm fr.}$ ) in the flamelet libraries by not only including an integrated heat loss rate from the flamelets, but also a translation of the flamelet in enthalpy space. This translation is simply denoted $h^\star$ , where the conserved enthalpy line is shifted as $h_c = h_{\rm ox} (1-Z) + h_{\rm fuel} Z + h^\star$ . This allows a full description of wall boundary heat loss. Having two heat loss parameterizations, however, makes the lookup procedure non-unique, requiring a method for deciding which point on $(\gamma, h^\star)$ to use for the flamelet lookup. We prefer $\gamma$ and use $\gamma$ as much as possible. When $\gamma$ is insufficient, which would is the case for overly cold or hot walls (wall temperatures outside of the range of temperatures spanned by the solutions of Eq. (3.372)). At the wall boundaries, the conserved enthalpy is defined to not be affected by heat loss while the true enthalpy is, providing $\gamma$ at the wall.

3.2.11.2.2. Property Table Heat Loss Parameterization

For nonadiabatic flamelet library generation and tabulation, a functional form for the heat loss parameter $\gamma$ in terms of reference quantities is required, similar in concept to the form of $\chi$ in (3.357). The value of $\gamma$ must be zero in each of the pure streams, and should have a maximum value near the stoichiometric flame sheet since this quantity typically represents radiative losses to the environment. A piecewise linear functional form is selected for simplicity. For a single mixture fraction, this form is simply

(3.373) $\gamma = \gamma_o \, F_\gamma(Z,Z_o),$

where $\gamma_o$ is a reference heat loss at reference state $Z_o$ (selected to be the stoichiometric condition $Z_\mathrm{st}$ ) and the nondimensional function $F_\gamma(Z,Z_o)$ is defined as

(3.374) $F_\gamma(Z,Z_o) = \left\{ \begin{array}{l@{\quad:\quad}l} \frac{Z}{Z_o} Z \le Z_o \\ \frac{1-Z}{1-Z_o} Z > Z_o \\ \end{array} \right. .$

For multiple mixture fractions, $\gamma$ is calculated by

(3.375) $\gamma = \gamma_o \, F_\gamma(Z_m,Z_{o,km},\gamma_{o,k}^\mathrm{max}),$

where $\gamma_o$ is the maximum-magnitude reference heat loss in the vector $\gamma_{o,k}^\mathrm{max}$ , which contains the reference heat loss parameters corresponding to maximum thermal losses for the $K$ stoichiometric mixture fractions that can be defined between stream pairs, $Z_{o,km}$ . The multiple stoichiometric mixture fractions are necessary because a single unique stoichiometric mixture fraction does not exist when using multiple mixture fractions.

The functional form for $F_\gamma$ is quite complex for multiple mixture fractions, and will only be described briefly here. In general, for a three-stream problem, there are two independent mixture fractions and the realizable mixture fraction space is the triangle where the two mixture fractions sum to a value less than or equal to unity. The value of $F_\gamma$ must be zero at the “corners” of this space, where the coordinates are $(0,0)$ , $(0,1)$ , and $(1,0)$ . The multiple stoichiometric mixture fractions between stream pairs will define points along the boundaries of this realizable mixture fraction space that represent local maxima in the heat loss distribution along that boundary. Straight lines may be used to connect these points in mixture fraction space, forming a “ridge” in the multidimensional $F_\gamma$ distribution. When definable, a linear fit is used between this ridge and a corner where $\gamma$ is zero. When not uniquely definable, linear fits are used between the ridge and the adjacent boundary value along rays extended from the opposite corner of the state space. Note that the values $\gamma_{o,k}$ are required for the calculation of $F_\gamma$ so that the final function may be normalized to a unity maximum value with appropriate relative scaling between the boundary heat loss values. Note that no more than a three-stream configuration is currently supported by fuego_tabular_props.

Applying the filtering operation in (3.365) to both sides of (3.375) yields

(3.376) $\tilde{\gamma} & = \int\limits_{-\infty}^{\infty} \int\limits_{-\infty}^{\infty} \int\limits_0^\infty \int\limits_0^1 \int\limits_0^1 \gamma_o \, F_\gamma(Z_m,Z_{o,km},\gamma_{o,k}^\mathrm{max}) \, p_{Z_i}(Z_i;\tilde{Z}_i,\widetilde{Z''^2}) \, p_{Z_m}(Z_m;\tilde{Z}_m) p_\chi(\chi;\tilde{\chi}) \, p_\gamma(\gamma;\tilde{\gamma}) \, p_{h^*}(h^*;\tilde{h}^*) \, \mathrm{d}Z_i \, \mathrm{d}Z_{m \neq i} \, \mathrm{d} \chi \, \mathrm{d} \gamma \, \mathrm{d} h^* & = \gamma_o \int_{-\infty}^{\infty} p_{h^*}(h^*;\tilde{h}^*) \; \mathrm{d}h^* \int_{-\infty}^{\infty} p_\gamma(\gamma;\tilde{\gamma}) \; \mathrm{d}\gamma \int_0^\infty p_\chi(\chi;\tilde{\chi}) \; \mathrm{d}\chi \int\limits_0^1 \int\limits_0^1 F_\gamma(Z_m,Z_{o,km},\gamma_{o,k}^\mathrm{max}) \, p_{Z_i}(Z_i;\tilde{Z}_i,\widetilde{Z''^2}) \, p_{Z_m}(Z_m;\tilde{Z}_m) \; \mathrm{d}Z_i \, \mathrm{d}Z_{m \neq i} & = \gamma_o \int\limits_0^1 \int\limits_0^1 F_\gamma(Z_m,Z_{o,km},\gamma_{o,k}^\mathrm{max}) \, p_{Z_i}(Z_i;\tilde{Z}_m,\widetilde{Z''^2}) \, p_{Z_m}(Z_m;\tilde{Z}_m) \; \mathrm{d}Z_i \, \mathrm{d}Z_{m \neq i},$

so that the filtered heat loss parameter can be calculated from the filtered quantities provided by the turbulent flame simulation as

(3.377) $\gamma_o(\tilde{Z}_m,\widetilde{Z''^2},\tilde{\gamma}) = \frac{ \tilde{\gamma} }{ \int\limits_0^1 \int\limits_0^1 F_\gamma(Z_m,Z_{o,km},\gamma_{o,k}^\mathrm{max}) \, p_{Z_i}(Z_i;\tilde{Z}_i,\widetilde{Z''^2}) \, p_{Z_m}(Z_m;\tilde{Z}_m) \; \mathrm{d}Z_i \, \mathrm{d}Z_{m \neq i} }.$

To decrease computational cost, the integral in the denominator can be interpolated from pre-calculated values in a multi-dimensional table as a function of $\tilde{Z}_m$ and $\widetilde{Z''^2}$ . (3.377) can be used during a simulation to convert filtered independent variables to the reference heat loss parameter required to perform table lookups to retrieve $\tilde{\phi}_T$ .

If turbulence/chemistry interactions are to be neglected in the simulation, the delta function $\delta(Z_i - \tilde{Z}_i)$ may be used for $p_{Z_i}(Z_i;\tilde{Z}_i,\widetilde{Z''^2})$ in (3.377) so that the reference heat loss parameter can be computed simply as

(3.378) $\gamma_o(\tilde{Z}_m,\tilde{\gamma}) = \frac{ \tilde{\gamma} }{ F_\gamma(\tilde{Z}_m,Z_{o,km}, \gamma_{o,k}^\mathrm{max}) }.$

3.2.11.2.3. Property Table Conserved Enthalpy Parameterization

For nonadiabatic flamelet library generation and tabulation, a functional form for the conserved enthalpy $h^*$ in terms of reference quantities is required. The value of $h^*$ should vary linearly within the range provided for each of the pure streams as a function of a reference heat loss parameter $h_o^*$ , with an appropriate stream-weighted blending for all other compositions.

The stream-weighted mixture properties are computed with an augmented mixture fraction vector $Z'_n$ in terms of $Z_m$ ,

(3.379) $Z'_n = \left[ Z_1, Z_2, \ldots, Z_M, 1-\sum_{m=1}^M Z_m \right],$

where the last component is simply the last implied mixture fraction to recover a unity sum. A reference augmented mixture fraction is defined as the centroid of the realizable mixture fraction space with each component being identical and equal to

(3.380) $Z'_{o,n} = \frac{1}{N}.$

From these definitions, minimum and maximum reference conserved enthalpy values may be computed as

(3.381) $h_{o,\mathrm{min}}^* = \sum_{n=1}^N h_{\mathrm{stream},\mathrm{min},n}^* \, Z'_{o,n}$

(3.382) $h_{o,\mathrm{max}}^* = \sum_{n=1}^N h_{\mathrm{stream},\mathrm{max},n}^* \, Z'_{o,n},$

where $h_{\mathrm{stream},\mathrm{min},n}^*$ and $h_{\mathrm{stream},\mathrm{max},n}^*$ are vectors of the minimum and maximum conserved enthalpy in pure stream $n$ , respectively. The conserved enthalpy can then be modeled as

(3.383) $h^* = h_{\mathrm{min},Z}^* + \left( h_o^* - h_{o,\mathrm{min}}^* \right) a_Z,$

where the mixture-weighted minimum conserved enthalpy is

(3.384) $h_{\mathrm{min},Z}^* = \sum_{n=1}^N h_{\mathrm{stream},\mathrm{min},n}^* \, Z'_n$

and the mixture-weighted stream variation proportionality constant is

(3.385) $a_Z = \sum_{n=1}^N \left( \frac{ h_{\mathrm{stream},\mathrm{max},n}^* - h_{\mathrm{stream},\mathrm{min},n}^* }{ h_{o,\mathrm{max}}^* - h_{o,\mathrm{min}}^* } \right) Z'_n.$

Applying the filtering operation in (3.365) to both sides of (3.383) yields

(3.386) $\tilde{h}^* = \tilde{h}_{\mathrm{min},Z}^* + \left( h_o^* - h_{o,\mathrm{min}}^* \right) \tilde{a}_Z,$

where the two mixture-weighted quantities are now expressed in terms of the augmented filtered mixture fraction as

(3.387) $\tilde{h}_{\mathrm{min},Z}^* = \sum_{n=1}^N h_{\mathrm{stream},\mathrm{min},n}^* \, \tilde{Z}'_n$

and

(3.388) $\tilde{a}_Z = \sum_{n=1}^N \frac{ h_{\mathrm{stream},\mathrm{max},n}^* - h_{\mathrm{stream},\mathrm{min},n}^* }{ h_{o,\mathrm{max}}^* - h_{o,\mathrm{min}}^* } \, \tilde{Z}'_n.$

This allows the reference conserved enthalpy to be expressed in terms of the filtered quantities provided by the turbulent flame simulation as

(3.389) $h_o^*(\tilde{Z}_m,\tilde{h}^*) = h_{o,\mathrm{min}}^* + \frac{ \tilde{h}^* - \tilde{h}_{\mathrm{min},Z}^* }{ \tilde{a}_Z }.$

3.2.11.3. Filtered Scalar Dissipation Rate

3.2.11.3.1. RANS Model

For RANS turbulence closure models the instantaneous laminar scalar dissipation rate given in (3.344) can be Favre-filtered and expanded to the form

(3.390) $\bar{\rho} \tilde{\chi} & = 2 \overline{ \rho D \frac{\partial{Z}}{\partial x_i} \frac{\partial{Z}}{\partial x_i} } & = 2 \overline{ \rho D } \frac{\partial{\tilde{Z}}}{\partial x_i} \frac{\partial{\tilde{Z}}}{\partial x_i} + 4 \overline{\rho D \frac{\partial{Z''}}{\partial x_i} } \frac{\partial{\tilde{Z}}}{\partial x_i} + 2 \overline{\rho D \frac{\partial{Z''}}{\partial x_i} \frac{\partial{Z''}}{\partial x_i} }.$

The middle term on the RHS is neglected for constant density flow [56]. The first term is referred to as the mean scalar dissipation rate

(3.391) $\bar{\rho} \tilde{\chi}_m = 2 \overline{ \rho D } \frac{\partial{\tilde{Z}}}{\partial x_i} \frac{\partial{\tilde{Z}}}{\partial x_i}$

and the third term is the perturbation scalar dissipation rate $\bar{\rho}\tilde{\chi}_p$ . This term can be modeled as

(3.392) $\bar{\rho} \tilde{\chi}_p = 2 \overline{\rho D \frac{\partial{Z''}}{\partial x_i} \frac{\partial{Z''}}{\partial x_i} }$

(3.393) $\approx C_\chi \bar{\rho}\, \frac{\epsilon}{k} \widetilde{Z''^2}$

for RANS-based turbulence closures where $\frac{\epsilon}{k}$ provides an inverse turbulence time scale, $\widetilde{Z''^2}$ is the scalar variance that will be modeled in Filtered Scalar Variance, and $C_\chi$ is a model constant that typically has a value of $2.0$ . [53]

Expressing the molecular mass diffusivity as $\bar{\rho} D = \mu/\mathrm{Sc}$ , where $\mu$ is the molecular viscosity and $\mathrm{Sc}$ is the Schmidt number, the modeled total filtered scalar dissipation rate for RANS closures is

(3.394) $\tilde{\chi} & = \tilde{\chi}_m + \tilde{\chi}_p & \approx \frac{2}{\bar{\rho}} \frac{\mu}{\mathrm{Sc}} \frac{\partial \tilde{Z}}{ \partial x_i} \frac{\partial \tilde{Z}}{\partial x_i} + C_\chi \frac{\epsilon}{k} \widetilde{Z''^2}.$

3.2.11.3.2. LES Model

For LES closures (3.394) also applies, so that the total filtered scalar dissipation rate is the sum of the mean and the perturbation scalar dissipation rates. The mean scalar dissipation rate is expressed identically to RANS closures as

(3.395) $\tilde{\chi}_m = \frac{2}{\bar{\rho}} \frac{\mu}{\mathrm{Sc}} \frac{\partial \tilde{Z}}{\partial x_i} \frac{\partial \tilde{Z}}{\partial x_i}.$

The perturbation scalar dissipation rate $\tilde{\chi}_p$ represents the sub-filter dissipation of scalar variance, and can be modeled by assuming that sub-filter dissipation is in local equilibrium with sub-filter production, and that the sub-filter production can be modeled with a gradient transport assumption as [57]

(3.396) $\bar{\rho} \tilde{\chi}_p = 2 \overline{ \rho D \frac{\partial Z''}{\partial x_i} \frac{\partial Z''}{\partial x_i} } = - 2 \overline{ \rho u_i'' Z'' } \frac{\partial \tilde{Z}}{ \partial x_i} \bar{\rho} \tilde{\chi}_p \approx 2 \frac{\mu_t}{\mathrm{Sc}_t} \frac{\partial \tilde{Z}}{\partial x_i} \frac{\partial \tilde{Z}}{\partial x_i},$

where $\mu_t$ is the modeled turbulent eddy viscosity and $\mathrm{Sc}_t$ is the turbulent Schmidt number.

This results in the final modeled form for the filtered total scalar dissipation rate for LES closures,

(3.397) $\tilde{\chi} &= \tilde{\chi}_m + \tilde{\chi}_p & \approx \frac{2}{\bar{\rho}} \left( \frac{\mu}{\mathrm{Sc}} + \frac{\mu_t}{\mathrm{Sc}_t} \right) \frac{\partial \tilde{Z}}{\partial x_i} \frac{\partial \tilde{Z}}{\partial x_i}.$

3.2.11.4. Filtered Mixture Fraction

The primary quantity used to identify the chemical state in Flamelet closure models is the mixture fraction, $Z$ . While there are many different definitions of the mixture fraction that have subtle variations that attempt to capture effects like differential diffusion, they can all be interpreted as a local mass fraction of the chemical elements that originated in the fuel stream. [58] The mixture fraction is a conserved scalar that varies between $0$ in the oxidizer stream and $1$ in the fuel stream, and is transported in laminar flow by the equation

(3.398) $\frac{\partial \rho Z}{\partial t} + \frac{ \partial \rho u_i Z }{ \partial x_i} = \frac{\partial}{\partial x_i} \left( \rho D \frac{\partial Z}{\partial x_i} \right),$

where $D$ is an effective molecular mass diffusivity.

Applying either temporal Favre filtering for RANS-based treatments or spatial Favre filtering for LES-based treatments yields

(3.399) $\frac{\partial \bar{\rho} \tilde{Z}}{\partial t} + \frac{ \partial \bar{\rho} \tilde{u}_i \tilde{Z} }{ \partial x_i} = - \tau_{Zu_j} + \frac{\partial}{\partial x_i} \left( \bar{\rho} D \frac{\partial \tilde{Z}}{\partial x_i} \right),$

where sub-filter correlations have been neglected in the molecular diffusive flux vector [59] and the turbulent diffusive flux vector is defined as

(3.400) $\tau_{Z u_j} \equiv \bar{\rho} \left( \widetilde{Z u_i} - \tilde{Z} \tilde{u}_i \right).$

Similar to species transport, this sub-filter correlation is modeled in both RANS and LES closures with the gradient transport approximation

(3.401) $\tau_{Z u_j} \approx - \bar{\rho} D_t \frac{\partial Z}{\partial x_i},$

where $D_t$ is the turbulent mass diffusivity, modeled as $\bar{\rho} D_t = \mu_t / \mathrm{Sc}_t$ where $\mu_t$ is the modeled turbulent viscosity from momentum transport and $\mathrm{Sc}_t$ is the turbulent Schmidt number. Please see the Fuego theory manual for further details. The molecular mass diffusivity is then expressed similarly as $\bar{\rho} D = \mu / \mathrm{Sc}$ so that the final modeled form of the filtered mixture fraction transport equation is

(3.402) $\frac{\partial \bar{\rho} \tilde{Z}}{\partial t} + \frac{ \partial \bar{\rho} \tilde{u}_i \tilde{Z} }{ \partial x_i} = \frac{\partial}{\partial x_i} \left[ \left( \frac{\mu}{\mathrm{Sc}} + \frac{\mu_t}{\mathrm{Sc}_t} \right) \frac{\partial \tilde{Z}}{\partial x_i} \right].$

In integral form as used in Fuego, the mixture fraction transport equation is

(3.403) $\int \frac{\partial \bar{\rho} \tilde{Z}}{\partial t}\, dV + \int \bar{\rho} \tilde{u}_i \tilde{Z} n_i\, dS = \int \left( \frac{\mu}{\mathrm{Sc}} + \frac{\mu_t}{\mathrm{Sc}_t} \right) \frac{\partial \tilde{Z}}{\partial x_i} n_i\, dS.$

3.2.11.5. Filtered Scalar Variance

3.2.11.5.1. RANS Model

For RANS-based turbulence closures, a transport equation is solved for the filtered scalar variance, $\widetilde{Z''^2}$ . This equation can be derived by subtracting (3.399) multiplied by $\tilde{Z}$ from the filter of the multiple of (3.398) and $Z$ , yielding

(3.404) $\frac{\partial \bar{\rho} \widetilde{Z''^2}}{\partial t} + \frac{\partial}{\partial x_i} \left( \bar{\rho} \tilde{u}_i \widetilde{Z''^2} \right) = - \frac{\partial}{\partial x_i} \left( \overline{\rho u_i'' Z''^2} \right) + \frac{\partial}{\partial x_i} \left( \overline{\rho D \frac{\partial Z''^2}{\partial x_i} } \right) + 2 \overline{ Z''^2 \frac{\partial}{\partial x_i} \left( \rho D \frac{\partial \tilde{Z}}{\partial x_i} \right) } - 2 \overline{ \rho u''_i Z'' } \frac{\partial \tilde{Z}}{\partial x_i} - 2 \overline{ \rho D \frac{\partial Z''}{\partial x_i} \frac{\partial Z''}{\partial x_i} },$

where the filtered mixture fraction variance is defined as $\widetilde{Z''^2} \equiv \widetilde{Z^2} - \tilde{Z}^2$ .

All five terms on the RHS of (3.404) require closure models. The first term represents turbulent transport of mixture fraction variance, and is modeled by a gradient-transport assumption as

(3.405) $- \overline{\rho u_i'' Z''^2} \approx \frac{\mu_t}{\mathrm{Sc}_t} \frac{\partial \widetilde{Z''^2}}{ \partial x_i}.$

The second and third terms on the RHS of (3.404) taken together represent molecular diffusion of mixture fraction variance, and is typically neglected with respect to turbulent transport for sufficiently high Reynolds numbers. Its effects are included here with another gradient-transport assumption of the form

(3.406) $\frac{\partial}{\partial x_i} \left( \overline{\rho D \frac{\partial Z''^2}{\partial x_i} } \right) + 2 \overline{ Z''^2 \frac{\partial}{\partial x_i} \left( \rho D \frac{\partial \tilde{Z}}{\partial x_i} \right) } \approx \frac{\partial}{\partial x_i} \left( \frac{\mu}{\mathrm{Sc}} \frac{\partial \widetilde{Z''^2}}{\partial x_i} \right).$

The fourth and fifth terms on the RHS of (3.404) represent production and dissipation of mixture fraction variance, respectively. The production term is similarly modeled with a gradient transport assumption as

(3.407) $- 2 \overline{ \rho u''_i Z'' } \frac{\partial \tilde{Z}}{\partial x_i} \approx 2 \frac{\mu_t}{\mathrm{Sc}_t} \frac{\partial \tilde{Z}}{\partial x_i} \frac{\partial \tilde{Z}}{\partial x_i}.$

The mixture fraction variance dissipation rate term is equal to the perturbation scalar dissipation rate,

(3.408) $2 \overline{ \rho D \frac{\partial Z''}{\partial x_i} \frac{\partial Z''}{\partial x_i} } = \bar{\rho} \tilde{\chi}_p,$

previously defined in (3.392) and modeled in (3.393). An identical treatment of this term is used here.

The final modeled form of the filtered scalar variance transport equation for RANS turbulence closure models is

(3.409) $\frac{\partial \bar{\rho} \widetilde{Z''^2}}{\partial t} + \frac{\partial}{\partial x_i} \left( \bar{\rho} \tilde{u}_i \widetilde{Z''^2} \right) = \frac{\partial}{\partial x_i} \left[ \left( \frac{\mu}{\mathrm{Sc}} + \frac{\mu_t}{\mathrm{Sc}_t} \right) \frac{\partial \widetilde{Z''^2}}{ \partial x_i} \right] + 2 \frac{\mu_t}{\mathrm{Sc}_t} \frac{\partial \tilde{Z}}{\partial x_i} \frac{\partial \tilde{Z}}{\partial x_i} - \bar{\rho} \tilde{\chi}_p.$

3.2.11.5.2. LES Model

For LES turbulence closures, the filtered scalar variance $\widetilde{Z''^2}$ can be modeled with the scaling law [57]

(3.410) $\bar{\rho} \widetilde{Z''^2} \approx C_V \bar{\rho} \Delta^2 \frac{\partial \tilde{Z}}{\partial x_i} \frac{\partial \tilde{Z}}{\partial x_i},$

where $\Delta$ is a length scale corresponding to the grid filter size and $C_V$ is a model coefficient. For the $k^\mathrm{sgs}$ closure and the non-dynamic Smagorinsky closure, $C_V$ has a fixed value of $0.5$ . For the dynamic Smagorinsky LES closure, $C_V$ can be dynamically calculated based on the local instantaneous flow-field.

To dynamically evaluate the filtered scalar variance model coefficient, begin by defining the grid filter-scale correlation

(3.411) $\tau_{Z''^2} & \equiv \bar{\rho} \widetilde{Z''^2} & = \bar{\rho} \widetilde{Z^2} - \bar{\rho} \tilde{Z}^2 & = \overline{\rho Z^2} - \frac{ \left(\overline{\rho Z}\right)^2 }{ \bar{\rho} }.$

Similarly, define an equivalent correlation at a larger test-filter scale

(3.412) $T_{Z''^2} \equiv \widehat{\overline{\rho Z^2}} - \frac{ \left( \widehat{\overline{\rho Z}} \right)^2 }{ \hat{\bar{\rho}} }.$

Now, define the quantity $L_{Z''^2}$ as a combination of these two correlations which reduces to an expression that can be evaluated in closed form,

(3.413) $L_{Z''^2} \equiv T_{Z''^2} - \widehat{\tau_{Z''^2}}$

(3.414) $L_{Z''^2} = \widehat{ \bar{\rho} \tilde{Z}^2 } - \frac{ \left( \widehat{ \bar{\rho} \tilde{Z} } \right)^2 }{ \hat{\bar{\rho}} }.$

By modeling the two correlations in (3.413) and equating them to (3.414), the model coefficient $C_V$ can be dynamically evaluated. The correlations at the two filter scales are modeled analogously as

(3.415) $\tau_{Z''^2} \approx C_V \bar{\rho} \Delta^2 \left( \frac{\partial \tilde{Z}}{\partial x_i} \right)^2$

(3.416) $T_{Z''^2} \approx C_V \hat{\bar{\rho}} \hat{\Delta}^2 \left[ \frac{\partial}{\partial x_i} \left( \frac{ \widehat{ \bar{\rho} \tilde{Z} } }{ \hat{\bar{\rho}} } \right) \right]^2,$

where $\hat{\Delta}$ is the characteristic test filter length scale and $C_V$ is assumed to be the same at both scales.

Notice that when the modeled forms of $\tau_{Z''^2}$ and $T_{Z''^2}$ are inserted into (3.413), $C_V$ appears inside a test filtering operation. Formally solving this system of equations for $C_V$ requires the expensive solution of an additional set of coupled integro-differential equations [60]. Alternatively, it is common practice to remove $C_V$ 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 can result in non-physical oscillations in the modeled value for $C_V$ . The square of the error involved in this approximation is $Q = (L_{Z''^2} - C_V M_{Z''^2})^2$ , where

(3.417) $L_{Z''^2} = \widehat{ \bar{\rho} \tilde{Z}^2 } - \frac{ \left( \widehat{ \bar{\rho} \tilde{Z} } \right)^2 }{ \hat{\bar{\rho}} }$

(3.418) $M_{Z''^2} = \hat{\bar{\rho}} \hat{\Delta}^2 \left[ \frac{\partial}{\partial x_i} \left( \frac{ \widehat{ \bar{\rho} \tilde{Z} } }{ \hat{\bar{\rho}} } \right) \right]^2 - \widehat{ \bar{\rho} \Delta^2 \left( \frac{\partial \tilde{Z}}{ \partial x_i} \right)^2 }.$

Minimizing this error in a least-squares fashion with respect to $C_V$ yields an expression for the modeled coefficient,

(3.419) $C_V = \frac{ L_{Z''^2} M_{Z''^2} }{ M_{Z''^2} M_{Z''^2} },$

that can be used directly in (3.410) for the filtered scalar variance.

Due to the above simplifications, the model coefficient $C_V$ can sometimes fluctuate wildly, possibly leading to numerical instabilities. A common solution to control these oscillations, and the one that is taken here, is to pass the numerator and denominator of (3.419) through a test filter, yielding

(3.420) $C_V = \frac{ \widehat{L_{Z''^2} M_{Z''^2}} }{ \widehat{M_{Z''^2} M_{Z''^2}} }.$

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