Compressible Flamelet Model

The approach employed in Loci-Stream utilizes a flamelet methodology extended to account for compressibility involving both ideal and real fluids. The ideas of this methodology were first presented by Ma et al. [MBHI2017]. The ideal-gas thermodynamics is modeled by linearizing the specific heat ratio whereas the parameters needed for the cubic Peng Robinson equation of state are pre-tabulated for the evaluation of departure functions and a quadratic expression is used to model the attraction parameter. This compressible model can account for temperature and pressure variations from the baseline flamelet table using a computationally tractable pre-tabulated combustion chemistry in a thermodynamically consistent fashion.

Governing Equations

Unlike the previous flamelet model in Loci-Stream [ThWrI2013] which was designed primarily to handle low-Mach number flow situations, the new compressible capability presented here handles all parts of the combustion domain, including the upstream fuel and oxidizer manifold sections, the combustion chamber, and the highly compressible nozzle. The governing equations used in this approach are the Favre-averaged Navier-Stokes equations, in addition to the transport equations for flamelet manifold variables Z and C, as given below.

(1)\[\frac{\partial \bar{\rho}}{\partial t} + \frac{\partial \bar{\rho} u_j}{\partial x_j} = 0\]
(2)\[\frac{\partial \bar{\rho} u_i}{\partial t} + \frac{\partial \bar{\rho} u_j u_i}{\partial x_j} = \ -\frac{\partial \bar{P}}{\partial x_i} + \frac{\partial}{\partial x_j} \left( \tilde{\tau}_{ij} - \overline{\rho u_i' u_j'} \right)\]
(3)\[\begin{split} \begin{align} \frac{\partial \bar{\rho}E}{\partial t} + \frac{\partial}{\partial x_j} (\bar{\rho} u_j E) &= \frac{\partial}{\partial x_j} \left[ \left( \frac{\lambda}{c_p} + \frac{\mu_t}{Pr_t} \right) \frac{\partial h}{\partial x_j} \right] + \frac{\partial}{\partial x_j} \left[ \sum_{k=1}^{NS} \left( \bar{\rho}D_k - \frac{\lambda}{c_p} \right) h_k \frac{\partial Y_k}{\partial x_j} \right] \\ &+ \frac{\partial}{\partial x_j} \left[ u_j \left( \tilde{\tau}_{ij} - \bar{\rho u_i' u_j'} \right) - u_j P \right] \end{align}\end{split}\]
(4)\[ \frac{\partial \bar{\rho} Z}{\partial t} + \frac{\partial \bar{\rho} u_j Z}{\partial x_j} = \frac{\partial}{\partial x_j} \left[ \left( \bar{\rho} D (EF) + \frac{\mu_t}{Sc_t} \right) \frac{\partial Z}{\partial x_j} \right]\]
(5)\[ \frac{\partial \bar{\rho} C}{\partial t} + \frac{\partial \bar{\rho} u_j C}{\partial x_j} = \frac{\partial}{\partial x_j} \left( \bar{\rho} D (EF) \frac{\partial C}{\partial x_j} \right) + \bar{\rho} \dot{\omega}_C \left( \frac{E}{F} \right)\]

where,

(6)\[ \tilde{\tau}_{ij} = \mu \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) - \frac{2}{3} \mu \frac{\partial u_k}{\partial x_k} \delta_{ij}\]
(7)\[ -\bar{\rho u_i' u_j'} = \mu_t \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) - \frac{2}{3} \mu_t \frac{\partial u_k}{\partial x_k} \delta_{ij} - \frac{2}{3} \bar{\rho} k \delta_{ij}\]

All thermochemical quantities are parameterized in terms of mixture fraction(Z) and progress variable(C), and the turbulence/chemistry interaction is modeled through a thickened flame closure. The parameterization can be represented as:

\[\tilde{Z}, \tilde{C} \xrightarrow{\text{flamelet lookup table}} \tilde{Y}_i\]

The presumed PDF closure is an approach where the turbulence closure is embedded within the flamelet tables. Other options are the Thickened Flame closure where the flamelet tables are agnostic as to the type of turbulence closure and the closure is handled within the code itself, as described in 3.3.

Thermodynamics and Transport Properties

Augmenting the above governing PDEs is a flamelet tabulation, based on the Peng-Robinson equation of state (PR-EoS) that provides both ideal-gas reference state and real-fluid thermodynamic information about the mixture (for a given Z and C) so that the density and temperature of the fluid may be recovered from the solved state variables (p, E, Z and C) using an efficient bracketed secant method iterative process. The Peng-Robinson equation of state [PeRo1976] is employed for the evaluation of thermodynamic quantities. It can be written as:

\[p = \frac{RT}{v-b} - \frac{a}{v^2 + 2bv - b^2}\]

where p is the pressure, R is the gas constant, T is the temperature, v is the specific volume. The attraction parameter(a) and effective molecular volume(b) are dependent on temperature and composition to account for effects of intermolecular forces. For mixtures, the parameters a and b are evaluated as:

\[a = \sum_{\alpha=1}^{N_s} \sum_{\beta=1}^{N_s} a_{\alpha\beta} X_{\alpha} X_{\beta}\]
\[b = \sum_{\alpha=1}^{N_s} X_{\alpha} b_{\alpha}\]

where \(X_\alpha\) is the mole fraction of species \(\alpha\). Extended corresponding states principle and single-fluid assumption for mixtures are adopted [5, 6]. The parameters a and b are evaluated using the recommended mixing rules by Harstad et al. [HaMiB2004].

\[a_{\alpha\beta} = 0.457236 \left( \frac{RT_{c,\alpha\beta}}{P_{c,\alpha\beta}} \right)^2 \left(1 + c_{\alpha\beta} \left[1 - \left(\frac{T}{T_{c,\alpha\beta}}\right)^{\frac{1}{2}}\right]\right)^2\]
\[b_{\alpha} = 0.077796 \frac{RT_{c,\alpha}}{P_{c,\alpha}}\]
\[c_{\alpha\beta} = 0.37464 + 1.54226\omega_{\alpha\beta} - 0.26992\omega_{\alpha\beta}^2\]

where \(T_{c,\alpha\beta}\) and \(P_{c,\alpha\beta}\) are the critical temperature and pressure of species \(\alpha\) and \(\beta\), respectively. The critical mixture conditions for temperature \(T_{c,\alpha\beta}\), pressure \(P_{c,\alpha\beta}\) and acentric factor \(\omega_{\alpha\beta}\) are determined using the corresponding state principles [PoPrC2001].

Partial derivatives and thermodynamic quantities based on the PR-EoS that are required for evaluating other thermodynamic variables can be derived analytically, as given below:

\[\left( \frac{\partial p}{\partial T} \right)_{v,X_i} = \frac{R}{v-b} - \frac{\left( \frac{\partial a}{\partial T} \right)_{X_i}}{v^2 + 2bv - b^2}\]
\[\left( \frac{\partial p}{\partial T} \right)_{T,X_i} = -\frac{RT}{(v-b)^2} \left\{ 1 - 2a \left[ RT(v+b) \left( \frac{v^2 +2bv -b^2}{v^2-b^2} \right)^2 \right] ^{-1} \right\}\]
\[\left( \frac{\partial a}{\partial T} \right)_{X_i} = -\frac{1}{T} \sum_{\alpha=1}^{N_s} \sum_{\beta=1}^{N_s} X_\alpha X_\beta a_{\alpha\beta} G_{\alpha\beta}\]
\[\left( \frac{\partial^2 a}{\partial T^2} \right)_{X_i} = 0.457236 \frac{R^2}{2T} \sum_{\alpha=1}^{N_s} \sum_{\beta=1}^{N_s} X_\alpha X_\beta c_{\alpha\beta} \left(1 + c_{\alpha\beta}\right) \frac{T_{c,\alpha\beta}}{P_{c,\alpha\beta}} \left( \frac{T}{T_{c,\alpha\beta}} \right)^{\frac{1}{2}}\]
\[G_{\alpha\beta} = \frac{c_{\alpha\beta} \left( \frac{T}{T_{c,\alpha\beta}} \right)^{\frac{1}{2}}} {1 + c_{\alpha\beta} \left(1-\frac{T}{T_{c,\alpha\beta}}\right)^{\frac{1}{2}}}\]
\[K_1 = \int_{+\infty}^{v} \frac{d\nu}{\nu^2 + 2b\nu - b^2} = \frac{1}{\sqrt{8b}} \ln \left( \frac{\nu + (1 - \sqrt{2})/b}{\nu + (1 + \sqrt{2})/b} \right)\]

For real fluids, thermodynamic quantities are typically evaluated from the ideal-gas value plus a departure function that accounts for the deviation from the ideal-gas behavior. The ideal-gas enthalpy, entropy and specific heat are evaluated from the NASA polynomials at a reference temperature of 298 K. The specific internal energy can be written as

\[e(T, \rho, X_i) = e^{ig}(T, X_i) + \int_{0}^{\rho} \left[ p - T \left( \frac{\partial p}{\partial T} \right)_{\rho, X_i} \right] \frac{d\rho}{\rho^2}\]

where superscript “ig” indicates the ideal-gas value of the thermodynamic quantity, and Eq. can be integrated analytically for PR-EoS to give:

\[e = e^{ig} + K_1 \left[ a-T \left( \frac{\partial a}{\partial T} \right)_{X_i} \right]\]

where \(K_1\) is evaluated using the expression above. The specific enthalpy can be expressed as:

\[h = h^{ig} - R T + K_1 \left[ a-T \left( \frac{\partial a}{\partial T} \right)_{X_i} \right] + p v\]

The specific heat capacity at constant volume and constant pressure, respectively, are evaluated as:

\[c_v = \left( \frac{\partial e}{\partial T} \right)_{v, X_i} = c^{ig} - K_1 T \left( \frac{\partial^2 a}{\partial T^2} \right)_{X_i}\]
\[c_p = \left( \frac{\partial h}{\partial T} \right)_{p, X_i} = c^{ig} - R - K_1 T \left( \frac{\partial^2 a}{\partial T^2} \right)_{X_i} - T \frac{\left( \frac{\partial p}{\partial T}\right)_{v,X_i}^2}{\left( \frac{\partial p}{\partial v}\right)_{v,X_i}}\]

The speed of sound for a real fluid is given by:

\[c^2 = \left( \frac{\partial p}{\partial \rho} \right)_{s, X_i} = \frac{\gamma}{\rho \kappa_T}\]

where \(\gamma\) is the specific heat ratio and \(\kappa_T\) is the isothermal compressibility defined as:

\[\kappa_T = -\frac{1}{v} \left( \frac{\partial v}{\partial p} \right)_{T, X_i}\]

The specific heat ratio is linearized around temperature to eliminate the costly iterative procedure to determine temperature, and to obtain other thermodynamic quantities which are functions of temperature. The underlying strategy rests on correcting the tabulated values with the transported quantities based on the EoS used. Specifically, because PR-EoS is employed, along with thermodynamic quantities needed for evaluation of the ideal gas thermodynamic quantities, parameters a and b, and the first and second derivatives of the parameter a w.r.t. temperature are needed for the calculations of the partial derivatives in the expressions derived above, which are required for the evaluation of the departure functions. However, the parameter a, along with its derivatives, is a function of both the species composition and the temperature, and thus may not be consistent with the temperature corresponding to the transported variables. The following procedure is adopted for the evaluation of the parameter a and its derivatives: the dependence of the parameter a on temperature is assumed to be a quadratic function as follows:

\[a = C_1 \tilde{T}^2 + C_2 \tilde{T} + C_3\]

where the coefficients \(C_1\), \(C_2\) and \(C_3\) are obtained from tabulated quantities:

\[C_1 = \frac{1}{2} \left( \frac{\partial^2 a}{\partial T^2} \right)_0\]
\[C_2 = \left( \frac{\partial a}{\partial T} \right)_0 - 2 C_1 T_0\]
\[C_3 = a_0 - C_1 T_0^2 - C_2 T_0\]

where subscript 0 indicates the stored baseline quantities in the table. The real-fluid energy is then evaluated as

\[\tilde{e} = \tilde{e}^{ig} + \tilde{e}^{dep}\]

where \(\tilde{e}^{ig}\) and \(\tilde{e}^{dep}\) are the ideal-gas and departure function values of the internal energy. The ideal-gas value including the chemical energy of the mixture is calculated with linearized specific heat ratio:

\[\tilde{e}^{ig} = \tilde{e}_0^{ig} + \frac{\tilde{R}}{a_{\gamma}^{ig}} ln \left( 1 + \frac{a_{\gamma}^{ig} \left(\tilde{T} - T_0 \right) }{\tilde{\gamma}_0^{ig} -1} \right)\]

where \(T_0\) , \(\tilde{e}_0^{ig}\), \(\tilde{R}\), \(a_{\gamma}^{ig}\) and \(\tilde{\gamma}_0^{ig}\) are parameterized with \(\tilde{Z}\), \(\tilde{Z}^{'' 2}\), and \(\tilde{C}\) and stored in the flamelet table. The departure function is given by:

\[\tilde{e}^{dep} = K_1 \left[ a - \tilde{T} \left( \frac{\partial a}{\partial T} \right)_{X_i} \right]\]

where the equations above detailing the quadratic fit of the parameter a are used to compute the parameters required for PR-EoS. Temperature and density are obtained by a bracketed secant iteration method from the computed pressure and energy expressions above.

Transport quantities are evaluated based on the method due to Chung et al [[ChLeS1984] [ChLlS1988]]. A power-law is used to approximate the temperature dependency:

\[\frac{\tilde{\mu}}{\mu_0} = \left( \frac{\tilde{T}}{T_0} \right)^{a_{\mu}}\]
\[\frac{\tilde{\lambda}}{\lambda_0} = \left( \frac{\tilde{T}}{T_0} \right)^{a_{\lambda}}\]

where \(\mu_0\) and \(\lambda_0\) are the reference viscosity and thermal conductivity evaluated at the reference temperature \(T_0\) in the flamelet table. The exponents \(a_{\mu}\) and \(a_{\lambda}\) are estimated using a finite-difference approximation of the slope change in the variables over a small temperature difference.

The new compressible flamelet methodology discussed above is thermodynamically consistent in the entire flow path of a rocket engine (from oxidizer and fuel manifolds to the exit of the nozzle) and circumvents the need for ad-hoc compressibility corrections of the FPV model in Stream. This new model is called the FPVC model.

Thickened Flame Closure

For this type of closure, the flamelet tables are parameterized by Z and C only, and the turbulence closure is handled within the code. See [ThWrI2019] for a detailed description of the thickened flame closure.

The basic idea in using Thickened Flame Closure is to move all turbulence closure for the chemical source terms out of the flamelet tables, bringing it instead into the governing PDEs being solved by the CFD solver. The concepts underlying thickened flame closure evolved from the analytical investigation of the fundamental conservation equations governing laminar premixed flames as well as the DNS solution of such flames [LePoV2000]. In general, when attempting to numerically simulate the propagation of a premixed flame on a grid of finite resolution, one must be able to resolve the flame to avoid numerical difficulties. Resolving the flame implies that the flame thickness must be captured over multiple cells in the grid. Generally, however, an actual flame is much thinner than the minimum cell size in the mesh. From basic analytical work [LePoV2000] it was found that if the molecular diffusivities in the governing equations are multiplied by a flame thickening factor, F, and at the same time the reaction rates in the governing equations are divided by the same factor, an equivalent thickened flame will result. This flame will be F times thicker than the original flame (and able to be resolved on the mesh) but will have the same flame speed as the original flame. Subsequent research [CoDuV2000], however, found that the resulting thickened flame did not respond in the same manner as the original flame to eddies in the flow field. Typically, a flame exposed to swirling eddies in the flow will stretch and fold (or wrinkle). As a result of this, the surface area of the flame front separating unburned fuel and oxidizer is increased, leading generally to an increased net reaction rate and thus an increased heat release per unit volume around the flame front. To counter the reduction in stretching and wrinkling of the thickened flame compared to the original flame, a so-called efficiency function, E, was introduced into the formulation which effectively enhances the reaction rate term in the governing equations. The resulting governing equations for mixture fraction and progress variable in the flamelet formulation, when using thickened flame closure, appear as follows, where multiplicative factors due to the closure model appear in the diffusion terms:

\[\frac{\partial \bar{\rho} Z}{\partial t} + \frac{\partial \bar{\rho} u_j Z}{\partial x_j} = \frac{\partial}{\partial x_j} \left[ \left( \bar{\rho} D (\text{EF}) + \frac{\mu_t}{Sc_t} \right) \frac{\partial Z}{\partial x_j} \right]\]
\[\frac{\partial \bar{\rho} C}{\partial t} + \frac{\partial \bar{\rho} \tilde{u}_j C}{\partial x_j} = \frac{\partial}{\partial x_j} \left( \bar{\rho} D(\text{EF}) \frac{\partial C}{\partial x_j} \right) + \bar{\rho} \dot{\omega}_C \left( \frac{E}{F} \right)\]

In the original thickened flame closure model [CoDuV2000], the flame thickening factor was specified as a constant value in the range \(5<F<30\). The efficiency function, E, on the other hand was determined on a point-by-point basis in the flow using a rather complex procedure designed to quantify the local flame wrinkling due to resolved and unresolved (sub-grid) flow motions and its effective increase in the local heat release. Vreman et al [VAOG2008] designed a new variation of the thickened flame closure approach which is substantially easier to implement and is more amenable to the use of a variety of turbulence models (RANS, DES, LES). Noting that the multiplicative factor EF which multiplies the molecular diffusivity in the C equation above could be tied to the underlying turbulence model (giving a spatially-varying thickening of the flame), this product is given by the following relation:

\[F^{1+\alpha} = \text{EF} = \frac{D + \mu_t / (\bar{\rho} Sc_t) + D_n}{D}\]

One can see that the product EF is essentially the ratio of the total diffusion to the molecular diffusion. This total diffusion strictly includes the original laminar diffusion (D), the turbulent diffusion arising from the use of a turbulence model (\(\mu_t\) is the eddy viscosity resulting from the turbulence model and \(Sc_t\) is the turbulent Schmidt number) as well as the diffusion (\(D_n\)) arising from the numerical treatment of the convection term in the governing equations (this last term is commonly neglected). As a result of this definition of EF no modification of the existing diffusion terms in Stream need be made as the existing turbulent diffusion coefficient is exactly the product D*EF (referring to the equations above). One needs simply to be able to compute the factor (\(\frac{E}{F}\)) which then multiplies the chemical source terms. To this end, Vreman et al [VAOG2008] proposed a simple model which relates the efficiency function E to the flame thickening factor as

\[E = F^{\alpha}, \quad 0 < \alpha < \frac{2}{3}\]

where the theoretical limits of \(\alpha\) were originally established by Colin [CoDuV2000]. Using this relation, the final expression for the factor (E/F) which multiplies the chemical source term in the governing equations of Loci-STREAM assumes the following form:

\[\frac{E}{F} = (F)^{\alpha - 1}\]

To visualize the effective smoothing that is being applied to the chemical source term, the table below shows values of (\(\frac{E}{F}\)) for various viscosity ratios (total viscosity/laminar viscosity) for given values of \(\alpha\). The general trend is evident whereby increasing viscosity ratios results in more damping of the chemical source term, whereas increasing flame wrinkling efficiency results in less damping of the chemical source term. For RANS simulations, in which the flame may be many cells thick, local damping factors can range anywhere from 0.1 < (\(\frac{E}{F}\)) < 0.5. For DES and LES simulations in which the viscosity ratios are much lower, damping factor will generally fall in the range 0.8 < (\(\frac{E}{F}\)) < 0.95. Reversion to proper DNS behavior is satisfied in the sense that as the eddy viscosity approaches zero, resulting in a viscosity ratio of unity, the damping factor (\(\frac{E}{F}\)) also approaches unity, implying no modification of the laminar reaction source term values.

Values of E for different values of viscosity ratio and \(\alpha\)

Viscosity ratio

E = \(F^0\)

E = \(F^{\frac{1}{2}}\)

E = \(F^{\frac{2}{3}}\)

0.01

0.991

0.997

0.998

0.1

0.913

0.970

0.982

1

0.513

0.801

0.875

10

0.095

0.457

0.624

20

0.050

0.369

0.549

50

0.021

0.274

0.459

100

0.010

0.218

0.401

References

[PoPrC2001]

Poling, Bruce E., John M. Prausnitz, and John P. O’Connell. 2001. Properties of Gases and Liquids. 5th ed. New York: McGraw-Hill Education.

[LePoV2000] (1,2)

Legier, J & Poinsot, Thierry & Veynante, D. (2000). Dynamically thickened flame LES model for premixed and non-premixed turbulent combustion. Proceedings of the Summer Program.