Numerical framework for transcritical real-fluid reacting flow simulations using the flamelet progress variable approach

TL;DR

This paper introduces a numerical framework for simulating transcritical real-fluid reacting flows using an extended flamelet progress variable approach, addressing pressure oscillations and large density ratios in high-pressure combustion scenarios.

Contribution

It develops a novel extension of the FPV model with a double-flux scheme and entropy-stable flux correction for robust transcritical flow simulations.

Findings

Successfully simulates cryogenic LOX/GH2 mixing and reacting flows.

Demonstrates robustness of the numerical schemes in multidimensional simulations.

Applicable for LES of real transcritical combustion applications.

Abstract

An extension to the classical FPV model is developed for transcritical real-fluid combustion simulations in the context of finite volume, fully compressible, explicit solvers. A double-flux model is developed for transcritical flows to eliminate the spurious pressure oscillations. A hybrid scheme with entropy-stable flux correction is formulated to robustly represent large density ratios. The thermodynamics for ideal-gas values is modeled by a linearized specific heat ratio model. Parameters needed for the cubic EoS are pre-tabulated for the evaluation of departure functions and a quadratic expression is used to recover the attraction parameter. The novelty of the proposed approach lies in the ability to account for pressure and temperature variations from the baseline table. Cryogenic LOX/GH2 mixing and reacting cases are performed to demonstrate the capability of the proposed approach in multidimensional simulations. The proposed combustion model and numerical schemes are directly applicable for LES simulations of real applications under transcritical conditions.

Equations

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

$$a$$

$$b$$

$$a_{\alpha \beta }$$

$$b_{\alpha }$$

$$c_{\alpha \beta }$$

$${\left(\frac{\partial p}{\partial T}\right)}_{v,X_i}$$

$${\left(\frac{\partial p}{\partial v}\right)}_{T,X_i}$$

$${\left(\frac{\partial a}{\partial T}\right)}_{X_i}$$

$${\left(\frac{{\partial }^{2}a}{\partial T^2}\right)}_{X_i}$$

$$G_{\alpha \beta }$$

$$K_1$$

$$e(T,\rho ,X_i)=e^{\text{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},$$

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

$$h=h^{\text{ig}}-RT+K_1\left[a-T{\left(\frac{\partial a}{\partial T}\right)}_{X_i}\right]+pv.$$

$$c_v={\left(\frac{\partial e}{\partial T}\right)}_{v,X_i}=c_v^{\text{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_p^{\text{ig}}-R-K_{1}T{\left(\frac{{\partial }^{2}a}{\partial T^2}\right)}_{X_i}-T\frac{(\partial p/\partial T{)}_{v,X_i}^2}{(\partial p/\partial v{)}_{T,X_i}}.$$

$$c^2={\left(\frac{\partial p}{\partial \rho }\right)}_{s,X_i}=\frac{\gamma }{\rho {\kappa }_T},$$

$${\kappa }_T=-\frac{1}{v}{\left(\frac{\partial v}{\partial p}\right)}_{T,X_i}.$$

$$-\frac{\rho \chi }{2}\frac{{\partial }^2{Y}_k}{\partial Z^2}$$

$$-\frac{\rho \chi }{2}\frac{{\partial }^{2}T}{\partial Z^2}$$

$$\frac{\partial \bar{\rho }}{\partial t}+\frac{\partial \bar{\rho }{\tilde{u}}_j}{\partial x_j}$$

$$\frac{\partial \bar{\rho }{\tilde{u}}_i}{\partial t}+\frac{\bar{\rho }{\tilde{u}}_i{\tilde{u}}_j}{\partial x_j}$$

$$\frac{\partial \bar{\rho }\tilde{E}}{\partial t}+\frac{\bar{\rho }{\tilde{u}}_j\tilde{E}}{\partial x_j}$$

$$\frac{\partial \bar{\rho }\tilde{Z}}{\partial t}+\frac{\bar{\rho }{\tilde{u}}_j\tilde{Z}}{\partial x_j}$$

$$\frac{\partial \bar{\rho }\widetilde{Z^{{}^{\prime \prime }2}}}{\partial t}+\frac{\bar{\rho }{\tilde{u}}_j\widetilde{Z^{{}^{\prime \prime }2}}}{\partial x_j}$$

$$\frac{\partial \bar{\rho }\tilde{C}}{\partial t}+\frac{\bar{\rho }{\tilde{u}}_j\tilde{C}}{\partial x_j}$$

$$a=C_1{\tilde{T}}^2+C_2\tilde{T}+C_3,$$

$$C_1$$

$$C_2$$

$$C_3$$

$$\tilde{e}={\tilde{e}}^{\text{ig}}+{\tilde{e}}^{\text{dep}},$$

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

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

$$\frac{\tilde{\mu }}{{\tilde{\mu }}_0}$$

$$\frac{\tilde{\lambda }}{{\tilde{\lambda }}_0}$$

$$\frac{\partial U}{\partial t}{V}_{cv}+\sum _f{F}^e{A}_f=\sum _f{F}^v{A}_f+S{V}_{cv},$$

Full text

1 Introduction

Liquid rocket engines (LRE) are one of the many practical applications which operate near or above the critical point of the working fluid. The large expansion ratio needed to generate the required thrust means that the combustion chamber operates at extremely high pressures, typically between 30 and 200 bar. In liquid rocket engines, the high mass flow rate is achieved by injecting high energy density cryogenic fuel and oxidizer into the combustion chamber through an extensive array of coaxial or impinging jets, see Cheroudi [1] for a review of high-pressure injection strategies. In the most common coaxial setup, the oxidizer generally exits the lip of the injectors below the critical temperature while being above the critical pressure of the fluid. In order to initiate chemical reactions for combustion, the cryogenic fluid must undergo a transcritical phase change to a supercritical state. In the transcritical regime, the thermo-physical properties of the fluid undergoes drastic changes for minute perturbations of the baseline thermodynamic state. This highly non-linear flow behavior requires the use of a generalized state equation to account for the complex thermodynamic properties.

A physics-based understanding of real fluid effects is needed to fully address the modeling challenges for trans- and supercritical combustion flow. Early experiments sought to address the physics of a sub- to supercritical jet phase change. The transcritical effects are most clearly observed in the visualization of transcritical pure mixing [2, 3, 4, 5]. At subcritical pressures, a liquid jet undergoes a convection driven process of atomization and breakup. At supercritical pressure, the same liquid jet undergoes a diffusion driven mixing —primarily due to the lack of surface tension, increased diffusivity and reduction in evaporation enthalpy. In the diffuse region between the liquid core and the gaseous outer flow, large thermo-physical gradients are observed and the dense core length is significantly reduced under supercritical pressures [6]. Under these conditions, a dense gas/gas mixing of the propellants and oxidizers is typically observed. The delineation between atomization and diffusion driven mixing of multi-species mixtures has been characterized by Dahms and Oefelein [7].

In high pressure rocket combustion chambers, the turbulence time scale has a magnitude of 1$\mu s$ in the reactive shear layer while the length scale is 1$\mu m$ near the exit lip [8]. The flame thickness, for non-premixed combustion, decreases with the inverse of the square root of pressure and strain rate [9]. In a follow-up study, Lacaze and Oefelein [10] suggested that the flame thickness is proportional to the inverse of the square root of pressure. The basis of the flamelet formulation [11] rests on a scale separation between the turbulence and the chemical scales in both time and space. Ivancic and Mayer [8] inferred that the turbulent and chemical time scales could be of similar magnitude —a result that could invalidate the laminar flamelet assumption. To address this issue, Zong et al. [12] showed the appropriateness of the flamelet assumptions in coaxial injectors in the supercritical regime using scaling arguments. The validity of the flamelet approach has lead to a number of flamelet-based numerically studied [13, 9, 14, 10, 15, 16].

Pre-tabulated flamelet approaches have been successfully applied to a variety of combustion cases [17, 11]. The original look-up tables based on mixture fraction, $\widetilde{Z}$, and its variance, $\widetilde{Z^{{}^{\prime\prime}2}}$, tend to be inadequate for describing topologically complex flames. This is particularly true for lifted flames where the mixture fraction alone is insufficient to model the full physical complexity [18, 19]. Careful experiments in high-pressure combustors have revealed a multiplicity of flame anchoring and stabilization mechanisms [20, 21], sometimes with a lifted flame configuration. In order to capture the intrinsic physics of these flames using a computationally tractable combustion modeling approach, a progress variable needs to be transported in addition to the mixture fraction. The flamelet progress variable (FPV) approach [22, 23] has been shown to better capture the complex physics in detached flames. The FPV approach has been used for supercritical simulations by Cutrone et al. [13] and more recently by Giorgi et al. [24] in the context of Reynolds averaged simulations. The extension of the FPV model for trans- and supercritical flows in the fully compressible context with large-eddy simulation (LES) remains, for the most part, unreported, specifically with regards to the pressure and temperature coupling.

The use of LES for trans- and supercritical flows has been initially explored using one-step chemistry [25] and more recently extended for hydrogen combustion by Schmitt et al. [26] using a thickened flame approach. Other groups have proposed an extension to the linear eddy model to account for combustion in high-pressure systems [27]. Further considerations have been proposed to account for subgrid-scale modeling under supercritical conditions [28, 29], sensitivity of the state equation [30], numerical stability issues [31, 32, 33, 34, 35, 36, 37], non-reflecting boundary conditions [38, 39], and compressibility effects [40]. While direct numerical simulations (DNS) have been successfully applied to trans- and supercritical flows in idealized configurations [41, 42, 43, 44], the computational limitations restrict the applications to academic problems.

In this work, we propose an extension to the classical FPV approach [22, 23] for trans- and supercritical combustion simulations in the context of finite volume, fully compressible, explicit solvers. The novelty of the present work lies in the ability to account for pressure and temperature variations from the baseline tabulated values using a computationally tractable pre-tabulated combustion chemistry in a thermodynamically consistent fashion. In addition, we show that the solution of the laminar flamelets in mixture fraction space and the chemistry tabulation requires special considerations in order to fully model the non-linear effects in transcritical flows.

2 Thermodynamics

2.1 Equation of state

Cubic equations of state offer an acceptable compromise between the conflicting requirements of accuracy and computational tractability. [45]. The Peng-Robinson (PR) cubic EoS [46] is used in this study for the evaluation of thermodynamic quantities, which can be written as

[TABLE]

where $p$ is the pressure, $R$ is the gas constant, $T$ is the temperature, $v$ is the specific volume, and 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 [47]

[TABLE]

where $X_{\alpha}$ is the mole fraction of species $\alpha$. Extended corresponding states principle and single-fluid assumption for mixtures are adopted. [48, 49] The parameters $a_{\alpha\beta}$ and $b_{\alpha}$ are evaluated using the recommended mixing rules by Harstad et al. [50]:

[TABLE]

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

2.2 Thermodynamic properties

Thermodynamic quantities in this study are evaluated consistently with respect to (w.r.t.) the EoS used and no linearization is introduced.

2.2.1 Partial derivatives

Partial derivatives and thermodynamic quantities based on the PR EoS that are useful for evaluating other thermodynamic variables are given as

[TABLE]

2.2.2 Internal energy and enthalpy

For general 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 can be evaluated from the commonly used NASA polynomials which have a reference temperature of 298 K. The simple mixture-averaged mixing rule is used for ideal-gas mixtures.

The specific internal energy can be written as,

[TABLE]

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

[TABLE]

in which $K_{1}$ is computed through Eq. 4f. The specific enthalpy can be evaluated from the thermodynamic relation $h=e+pv$, and we have

[TABLE]

Partial enthalpies for each species, $h_{k}$, can be evaluated using the partial derivatives w.r.t. mole fractions. The procedures for evaluating partial enthalpy for real fluids are similar to those in Meng et al. [51] and therefore omitted here.

2.2.3 Specific heat capacity

The specific heat capacity at constant volume is evaluated as

[TABLE]

and the specific heat capacity at constant pressure is evaluated as

[TABLE]

2.2.4 Speed of sound

The speed of sound for general real fluids can be evaluated as

[TABLE]

where $\gamma$ is the specific heat ratio and $\kappa_{T}$ is the isothermal compressibility, which is defined as

[TABLE]

2.3 Transport properties

The dynamic viscosity and thermal conductivity are evaluated using Chung’s high-pressure method [52, 53]. This method is known to produce oscillations in viscosity for multi-species mixtures that consist of species with both positive and negative acentric factors [54, 55]. To solve this problem, a mass-fraction averaged or mole-fraction averaged viscosity evaluated based on viscosity of each individual species can be used. In this study, the negative acentric factor is set to zero only when evaluating the viscosity so that the anomalies in viscosity can be removed. This approach has similar behavior to the mole-fraction averaged approach via numerical tests.

3 Transcritical Flamelets

The large Damköhler number of the supercritical combustion [12] supports the use of laminar flamelet-based combustion models. The basic assumption of the flamelet model rests on the fact that the reaction zone remains laminar and the diffusive transport is only important in the direction normal to the flame. By recasting the governing equations as a one-dimensional similarity solution in mixture fraction space, the steady flamelet equations for the species and temperature under unity Lewis number assumption can be written as:

[TABLE]

where $Y_{k}$ is the mass fraction of species $k$, $Z$ is the mixture fraction, $\dot{\omega}_{k}$ is the reaction rate of species $k$, $\dot{\omega}_{T}$ is the heat release term, and $\chi=2D|\nabla Z|^{2}$ is the scalar dissipation rate where $D$ is the diffusion coefficient. No further assumptions are needed to the flamelet equations for a generalized equation of state (for a constant pressure combustion) apart from the modified thermodynamic relationship between temperature and density. For a given profile of the scalar dissipation (which accounts for the convective and diffusive terms normal to the flame front), the solution of the above equations represents the composition and temperature profiles within the counterflow-diffusion flame. In the low-Mach and low-pressure formulation, Peters [56] derived a scalar dissipation rate profile for a shear layer under the assumption of a unitary Chapman-Rubesin parameter that relates the local ratio of density and viscosity to its far-field values. Although formally insufficient for many combustion cases, the modeling assumptions of the scalar dissipation rate profile are robust. For transcritical conditions, we retain the same scalar dissipation rate profile.

The flamelet equations of a counterflow-diffusion flame, in mixture fraction space, are solved using the FlameMaster solver [57]. The original code was extended to incorporate the PR EoS, along with the thermodynamically consistent departure functions and thermo-physical fluid properties. The results of the one-dimensional problem are compared with the two-dimensional direct numerical simulations performed by Lacaze and Oefelein [10]. In order to avoid inconsistencies in determining the scalar dissipation profile, only the near-equilibrium flamelet solutions are compared. The setup consists of a counterflow diffusion flame with pure H2 and O2 streams, respectively, at 295 K and 120 K and a pressure of 7.0 MPa. The strain rate, defined as the velocity difference between both injectors, is set to $10^{5}$ s*-1*, which corresponds to a scalar dissipation rate of $10^{3}$ s*-1*, following the classical relationship between the strain and the scalar dissipation rate in laminar flamelets derived by Peters [56]. The high-pressure chemical mechanism by Burke et al. [58] is used, which accounts for 8 species and 27 reactions. The comparative results for temperature, composition, density and specific heat capacity are shown in Fig. 1. Specificial attention needs to be paid to the transcritical regions away from the flame front, especifically in the oxidizer stream. In this region, large variations of the thermophysical properties are observed which require an adaptive mesh refinement based on the thermophysical properties in addition to the usual gradient-based mesh adaptation techniques. Without these strong non-linear features, the full complexity of the transcritical behavior cannot be captured. The grid adaptation is especially important to capture the pseudo-boiling point [3, 59] (PBP) which contains a finite peak in the specific heat capacity, and this region acts as a proxy to phase change in subcritical thermodynamics. Note that the location of the PBP region is typically at mixture fraction around $10^{-3}$ independent of the value of scalar dissipation rate. This requires special attention during the tabulation process since a typical resolution in mixture fraction will completely miss the PBP region, and this will be discussed in details later. As can be seen in Fig. 1, when the thermodynamic features are well resolved, the low-dimensional flamelet equations in mixture fraction space can capture the essential flame structure.

4 Numerical Methods

4.1 Governing equations

The governing equations are the Favre-averaged conservation equations of mass, momentum, total energy, mixture fraction, mixture fraction variance, and progress variable, written as follows:

[TABLE]

where $u_{i}$ is the $i$th component of the velocity vector, $E$ is the total energy including the chemical energy, $C$ is the progress variable, $\mu$ and $\mu_{t}$ are the laminar and turbulent viscosity, $\lambda$ is the thermal conductivity, $D$ is the diffusion coefficient for the scalars, $\dot{\omega}_{C}$ is the source term for the progress variable, $\tau_{ij}$ and $\tau_{ij}^{R}$ are the viscous and subgrid-scale stresses which are assumed to take the form as the second term on the right-hand side of Eq. 13b, Prt is the turbulent Prandtl number, and Sct is the turbulent Schmidt number. An appropriate subgrid-scale model is needed for the computation of the turbulent viscosity $\mu_{t}$. Under the unity Lewis number assumption, the summation on the right-hand side of Eq. 13c vanishes so that the species mass fractions are not explicitly required for the energy equation. The system is closed with the PR EoS introduced in Section 2 and the flamelet-based combustion model that will be discussed in the next subsection. Moreover, the sub-grid terms associated with the EoS are neglected in this study.

4.2 Combustion model

A Flamelet/Progress Variable (FPV) approach [22, 23] is adopted in this study, in which the chemistry is pre-computed and tabulated as a series of laminar flamelet solutions. Flamelets are first computed for different values of the scalar dissipation rate at a constant background pressure and specified constant fuel and air temperatures, and then the flamelets are parametrized by the mixture fraction and the reaction progress variable. The resulting flamelet table is used for the determination of the local temperature, species, density, source term of the progress variable and other thermal-transport quantities needed by the solver. Presumed PDFs are introduced to account for the turbulence/chemistry interaction. Typically, a $\beta$-PDF is used for the mixture fraction and a $\delta$-PDF for the progress variable, which was shown to be a reasonable approximation in many studies. Since the reacting region is typically in the ideal gas regime even under transcritical combustion conditions, the PDF closures are expected to perform similarly as in previous studies for ideal gas reacting flows.

In the low-Mach number flamelet implementation, the temperature, species, and density are assumed to depend only on the transported scalars. However, when compressibility effects are taken into account, an overdetermined thermodynamic state arises from the use of the flamelet table with tabulated thermodynamics. On the one hand, the full thermodynamic state, at constant pressure, is defined within the flamelet table. On the other hand, the transport equations contain two thermodynamic variables, namely density and internal energy, which are also sufficient to fully characterize the thermodynamic equilibrium state of the fluid if the compositions are given. In order to mend the over-determined thermodynamic states, a strategy was developed by Saghafian et al. [60] in the context of ideal gas flows to account for the pressure and temperature variations arising in supersonic combustion using the FPV approach. The specific heat ratio is linearized around temperature to eliminate the costly iterative procedure to determine temperature, and also to obtain other thermodynamic quantities which are functions of temperature.

The proposed strategy for applying compressibility effects in the FPV approach has been modified to work with a generalized equation of state under transcritical conditions in the present work. The underlying strategy rests on correcting the tabulated values with the transported quantities based on the EoS used. Specifically, since PR EoS is used in this study, along with thermodynamic quantities needed for evaluation of the ideal gas thermodynamic quantities, parameters $a$, $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 Eq. 4 that are needed for the evaluation of the departure functions. The parameter $b$ is a function of species composition only and is independent of pressure and temperature. Therefore, it can be pre-tabulated within the flamelet table without any discrepancy in pressure and temperature. 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 proposed 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,

[TABLE]

where $C_{1}$, $C_{2}$, and $C_{3}$ can be determined from tabulated quantities,

[TABLE]

where subscript 0 indicates the stored baseline quantities in the table. The first and second derivatives of parameter $a$ w.r.t. temperature can be determined accordingly by taking derivatives of Eq. 14, and the proposed model corresponds to a linear and a constant approximation to the first and second derivatives, respectively. Once the parameter $a$ and its derivatives are obtained, along with the parameter $b$ and the gas constant $R$, all the partial derivatives which are needed for computing other thermodynamic quantities can be evaluated for a given mixture, and therefore, the thermodynamic state is determined. Note that similar to the quadratic model as in Eq. 14, a linear model or a constant model for the parameter $a$ can also be constructed based on the stored values in the table. The performance of the quadratic, linear, and constant model will be examined later.

As an illustration, Fig. 2 shows the results of parameter $a$ and its first derivative w.r.t. temperature, along with the approximations evaluated based on the quadratic, linear, and constant model. Pure oxygen at 100 bar is considered for this example. The reference point is assumed to be at a temperature of 300 K, as indicated by the black dot in Fig. 2. As can be seen from Fig. 2, the quadratic assumption works well in the region within 200 K and 100 K from the reference temperature for the parameter $a$ and its derivative, respectively. The linear model yields a linear profile for $a$ and hence a constant profile for its derivative. The constant model for the parameter $a$ gives zero value for its derivative. Similarly, the behavior of the second derivative of $a$, which is also important for the evaluation of real-fluid thermodynamics, can be expected for the three models considered. The quadratic model for $a$ shows superior performance in predicting $a$ and its derivatives when temperature is away from the reference value, and its performance in predicting other thermodynamic quantities will be examined in details to confirm the validity of the proposed approach.

As an example, to calculate the internal energy including the chemical energy from the temperature based on the proposed approach for PR EoS for a given mixture, i.e. fixed $\widetilde{Z}$, $\widetilde{Z^{{}^{\prime\prime}2}}$, and $\widetilde{C}$, the ideal gas part and the departure function are calculated separately,

[TABLE]

where $\widetilde{e}^{\,\text{ig}}$ and $\widetilde{e}^{\,\text{dep}}$ are the ideal-gas and departure function values of the specific internal energy. The ideal-gas value including the chemical energy of the mixture is calculated with linearized specific heat ratio [60]

[TABLE]

where $\widetilde{e}_{0}^{\,\text{ig}}$, $\widetilde{R}$, $T_{0}$, $\widetilde{\gamma}_{0}^{\,\text{ig}}$, and $a_{\gamma}^{\text{ig}}$ are all stored variables in the table for given $\widetilde{Z}$, $\widetilde{Z^{{}^{\prime\prime}2}}$, and $\widetilde{C}$, in which $a_{\gamma}^{\text{ig}}$ is the slope of the ideal-gas specific heat ratio w.r.t. the temperature. The departure function is determined as

[TABLE]

where Eqs. 14 and 15 are used for the evaluation of parameters needed for the PR EoS. Other thermodynamic quantities, such as the specific heat and speed of sound, can be evaluated similarly.

To determine primitive variables from conservative variables, a secant method is used to obtain temperature given the transported density and internal energy. With this, the pressure along with other thermodynamic quantities is evaluated as an explicit function of the density and temperature. Since the parameters needed for the real fluid EoS are pre-tabulated and approximated from the table, the computational overhead of the iteration process is acceptable.

Transport quantities are evaluated based on the method due to Chung et al. [52, 53] Since the variation of the viscosity and conductivity on pressure is small under transcritical conditions, a power-law is used to approximate the temperature dependency as follows,

[TABLE]

where $\widetilde{\mu}_{0}$ and $\widetilde{\lambda}_{0}$ are stored in the table along with their corresponding slopes, $a_{\mu}$ and $a_{\lambda}$.

Note that the proposed approach is not limited to the PR EoS. All cubic EoS’s have similar structure and a similar approach can be used for other types of cubic EoS, such as the Soave-Redlich-Kwong (SRK) EoS [61].

The developed FPV approach focuses on the thermodynamics and no special modification is taken for the chemistry part. The transcritical combustion dynamics were shown to have similar structures as for ideal gases by several studies [9, 15, 16]. Indeed, in typical engines, combustion takes place at high temperatures away from the real-fluid region which is characterized by cryogenic temperatures. Moreover, for the applications considered in this study, combustion can be considered to be at low-Mach conditions where compressibility effects can be neglected. However, for the inert and equilibrium parts of the flows without chemical source terms, compressibility may play a critical role for predicting the behaviors of the system of interest, for example flows in the injectors and through the nozzles for rocket engines. Therefore, a strategy is proposed here for practical simulations, that the flamelet table can be generated at conditions of the combustion chamber where combustion takes place, and discrepancies in temperature and pressure in other parts of the flows can be corrected by the FPV model described above. If supersonic combustion needs to be considered, methodologies developed by Saghafian et al. [60] can be used.

4.3 Numerical schemes

The massively paralleled, finite-volume solver, CharLES${}^{\,x}$, developed at the Center for Turbulence Research, is used in this study. A control-volume based finite volume approach is utilized for the discretization of the system of equations, Eq. 13:

[TABLE]

where $U$ is the vector of conserved variables, $F^{e}$ is the face-normal Euler flux vector, $F^{v}$ is the face-normal viscous flux vector which corresponds to the r.h.s of Eq. 13, $S$ is the source term vector, $V_{cv}$ is the volume of the control volume, and $A_{f}$ is the face area. A strong stability preserving 3rd-order Runge-Kutta (SSP-RK3) scheme [62] is used for time advancement.

The convective flux is discretized using a sensor-based hybrid scheme in which a high-order, non-dissipative scheme is combined with a low-order, dissipative scheme to minimize the numerical dissipation introduced. A central scheme which is fourth-order on uniform meshes is used along with a second-order ENO scheme for the hybrid scheme. A density sensor [54, 37] is adopted in this study. Due to the large density gradients across the PBP region under transcritical conditions, an entropy-stable flux correction technique, developed by Ma et al. [37], is used to ensure the physical realizability of the numerical solutions and to dampen the non-linear instabilities in the numerical schemes.

Due to the strong non-linearlity inherited in the real fluid EoS, spurious pressure oscillations will be generated when a fully conservative scheme is used [63, 37], and severe oscillations in the pressure field could make the solver diverge which cannot be solved by adding artificial dissipation. A double-flux method [35, 36, 37] is extended to the transcritical regime to eliminate the spurious pressure oscillations. The effective specific heat ratio based on the speed of sound is frozen both spatially and temporally for a given cell when the fluxes of its faces are evaluated, which renders a local system as an equivalently ideal-gas system. The two fluxes at a face evaluated in the double-flux method are not the same, yielding a quasi-conservative scheme and the conservative error in total energy was shown to converge to zero with increasing resolution. [37] A Strang-splitting scheme [64] is applied in this study to separate the convection operator from the remaining operators of the system.

5 Results and Discussions

5.1 Tabulation approach analysis

The validity and performance of the FPV model with the tabulation approach developed in the previous sections will be assessed in the following through an a priori analysis. To this end, a flamelet table is firstly generated at certain reference conditions and then a flamelet solution at different reference conditions is used as a probe for the assessment. In a compressible solver, initial and boundary conditions are typically specified through primitive variables (temperature, pressure, velocity, and species). With the FPV approach, species are determined from mixture fraction and progress variable. The solver then converts primitive variables to conservative variables and takes time advancement. To be consistent with the compressible solver, the primitive variables, namely temperature, pressure, mixture fraction and progress variable, from the probing flamelet solution are given, and thermodynamic parameters, such as gas constant, specific heat ratio, parameter $a$ and $b$ for the PR EoS, are interpolated from the table by fixing mixture fraction and progress variable. Then the species and thermodynamic quantities of interest are compared to the exact values from the probing flamelet to assess the performance of the current developed tabulation approach.

An extreme case is considered here to show the capabilities of the current model to recover real-fluid thermodynamics. The flamelet table is generated at ideal-gas conditions with $T_{\text{ox}}$ = 300 K, $T_{\text{f}}$ = 300 K, $p$ = 60 bar. A transcritical flamelet in equilibrium at relatively low scalar dissipation rate at conditions $T_{\text{ox}}$ = 100 K, $T_{\text{f}}$ = 150 K, $p$ = 100 bar is considered as the probing flamelet. Note that the flamelet table contains only ideal-gas information and directly reading variables from the table is expected to give significant errors. This situation corresponds to the case when an ideal-gas table is used for transcritical simulations, or when the table resolution in mixture fraction is insufficient, so that the PBP region is not resolved.

Figure 3 shows the results of the a priori analysis. Species and thermodynamic quantities of interest evaluated from the current FPV model and directly read from the table are compared to the exact values from the probing flamelet. The current FPV model assumes that the composition is not changed with perturbations on the reference conditions, or in other words, species are only functions of the mixture fraction and the progress variable. Thus species from the FPV model are expected to have the same values as in the flamelet table. Quadratic, linear, and constant models for the parameter $a$ are compared in terms of real-fluid thermodynamic quantities. In the following, we will go through different aspects contained in this analysis.

The species results, including major species (H2, O2, and H2O) and representative minor species (OH, O, H, and HO2), are shown in the first three subfigures in Fig. 3. It can be seen that the major species read from the table exhibit negligible difference from the exact values even considering the dramatic difference in the pressure. For minor species, small but noticeable discrepancies can be observed, but the influence on the mixture properties is expected to be insignificant due to the low magnitude of their mass fractions. These results are consistent with the findings by Saghafian et al. [60] where flamelets at different operating conditions are extensively studied for ideal gases at pressures close to ambient conditions. Similar results are found for the entire S-shaped curve with probing flamelets at different scalar dissipation rates which are not shown here.

Density, specific heat at constant pressure, and speed of sound are used as examples for comparison of the thermodynamic quantities in Fig. 3. Since the flamelet table is generated at different reference conditions from the probing flamelet, the quantities directly read from the table show significant discrepancies from the flamelet. The table is at ideal-gas conditions, and therefore, the large density value of LOX and the sharp change of density profile in the PBP region are completely missed. Similarly, the peak in specific heat in the PBP region cannot be read from the table. The behavior in speed of sound is not correct as compared to the probing flamelet in the real-fluid region as indicated by the shaded area. In the ideal-gas region, which is indicated by the white area in the last subfigures of Fig. 3 (defined by 5% deviation of compressibility factor from unity), although the gas constant is well recovered by the species information, due to the drastic difference in temperature and pressure, the table gives significant errors especially in density. This is better depicted in Fig. 4, where relative errors in density and specific heat from different models are shown. As can be seen in Fig. 4, the density directly interpolated from the table has errors exceeding 80% and 40% in the real-fluid and ideal-gas regions, respectively.

In contrast, the current FPV model with a quadratic model for the parameter $a$ shows superior performance in recovering these quantities. In the ideal-gas region, a linearized specific heat ratio is used to compensate the temperature discrepancy from the table, and this has been shown to yield good predictions with a temperature discrepancy up to 500 K [65]. Similar behavior can be seen in the current study. As shown in Fig. 4 that the relative error from the FPV model in the ideal-gas region is less than 5%. In the real-fluid region, despite the fact that extrapolation from the table is needed for recovering the thermodynamic quantities (200 K difference of temperature for the oxidizer stream), the current FPV model shows significant improvements. The quadratic, linear and constant model yield increasingly more accurate results successively towards the flamelet solution as expected. The relative error of the quadratic model is below 10% and 20% for the density and specific heat as shown in Fig. 4. Obviously when the reference table contains real-fluid information by lowering the oxidizer temperature, and with grid points in mixture fraction clustered in the oxidizer side, the performance of the current FPV model can be further improved. The operating conditions are deliberately chosen to challenge the current FPV model.

In summary, the assumption that the composition is a weak function of the reference conditions is valid under transcritical conditions. The currently developed FPV approach with quadratic expression for the attraction parameter can accurately recover the real-fluid and ideal-gas thermodynamic quantities despite variations in temperature and pressure.

5.2 Cryogenic LOX/GH2 mixing case

The proposed FPV approach is tested with a benchmark case proposed by Ruiz et al. [55] for high-Reynolds number turbulent flows with large density ratios. A two-dimensional mixing layer of liquid-oxygen (LOX) and gaseous-hydrogen (GH2) streams is simulated. Details of the configuration and the computation domain for the simulation can be found in Ruiz et al. [55] The configuration is representative of a coaxial rocket combustor, in which dense LOX is injected in the center to mix with the surrounding high-speed GH2 stream. The two streams are separated by the injector lip which is also included in the computational domain. The computational domain has a dimension of $10h\times 10h$, where $h=0.5$ mm is the height of the injector lip. A sponge layer of length 5$h$ is put at the end of the domain. A fully structured Cartesian mesh is utilized in this case with 100 grid points across the injector lip. A uniform mesh is used in axial direction. For the region within 3$h$ around the injector lip in transverse direction, a uniform mesh is adopted and stretching is applied with a ratio of 1.02 outside this region. The mesh is stretched in axial direction in the sponge layer. This results in a mesh with a total of $5.2\times 10^{5}$ cells.

Adiabatic no-slip wall conditions are applied at the injector lip and adiabatic slip wall conditions are applied for the top and bottom boundaries of the domain. A 1/7th power law for velocity is used for both the LOX and GH2 streams. Pressure outlet boundary conditions are applied after the sponge layer where acoustic waves are suppressed. The LOX stream is injected at a temperature of 100 K, and GH2 is injected at a temperature of 150 K. The pressure is set to 10 MPa which is representative of rocket combustor conditions. Note that the density ratio between LOX and GH2 is about 80. The hybrid scheme using the double-flux model is used with the RS sensor[54, 37] set to a value of 0.2. A first-order upwinding scheme is used in the sponge layer. The CFL number is set to 1.0 and no subgrid scale model is used to facilitate comparisons with the reference solutions.

The flamelet table is generated at same conditions, namely $T_{\text{ox}}$ = 100 K, $T_{\text{f}}$ = 150 K, and $p$ = 100 bar. Transcritical flamelets are calculated for the entire S-shaped curve at reference conditions. Water is used as the progress variable. Resolutions in mixture fraction and progress variable are both 100 grid points. In mixture fraction dimension, the mesh is uniform from zero to the stoichiometric value using one third of the grid points and then stretches to one with the rest of the grid points. No special treatment is conducted for the resolution in the PBP region for the reason that the current developed FPV model is insensitive to the table resolution at this region as demonstrated in the previous subsection. For the pure mixing case considered in this subsection, the progress variable, $\widetilde{C}$, is set to zero.

Figure 5 shows results for instantaneous axial velocity, density, and oxygen mass fraction. Due to the implementation of the proposed FPV model, for pure mixing case, the FPV model is expected to perform almost the same as a multi-component model in which species mass fractions are explicitly solved. Indeed, the mixture fraction for the mixing case considered here acts as the mass fraction of hydrogen. The only difference comes from the evaluation of the thermodynamic quantities which is expected to be small as demonstrated in the previous subsection. The results from the current FPV model is found to give almost the same results both qualitatively and quantitatively as the multi-component model in Ma et al. [37] As can be seen in Fig. 5, the flow field is dominated by large vortical structures in the mixing layer and three large vortical structures are separated by waves with a wavelength of approximately $5h$. The predicted structures of vortices are in good agreement with those reported in Ruiz et al. [55]. From the density field, “comb-like” or “finger-like” structures [2] can clearly be seen, which was also observed through experiments for transcritical mixing under typical rocket engine operating conditions [3, 4].

The simulation results are averaged in time to facilitate quantitative comparisons and 15 flow-through-times are used for the averaging process after steady state is reached. One flow through time corresponds to 0.125 ms [55]. Figure 6 shows the results of mean and root mean square (RMS) values for axial velocity, mass fraction of oxygen, and temperature. Statistics at different axial locations ($x/h$ = 1, 3, 5, 7) are plotted as a function of normalized transverse distance. Results from the current solver, CharLES${}^{\,x}$, are compared to those obtained from two other solvers, namely $AVBP$ and $RAPTOR$ [55]. The mean axial velocity is in good agreement between the three different solvers, while there are some discrepancies in the RMS values. Results from CharLES${}^{\,x}$ show slightly lower RMS values on the GH2 side, especially for the axial location of $x/h=3$. This is probably due to the different implementation of the sponge layer and outlet boundary conditions adopted by the different solvers. Results for the oxygen mass fraction are almost identical for the three solvers except for the small difference seen at the GH2 side which could be related to the difference seen in the velocity results. Appreciable differences are observed in the temperature statistics. The results from CharLES${}^{\,x}$ and AVBP are similar and show a narrower thermal mixing layer as compared to that of RAPTOR. Overall, the results obtained from the three solvers are in good agreement, demonstrating the capability of the proposed FPV model and the robustness of the numerical schemes for transcritical mixing cases.

5.3 Cryogenic LOX/GH2 reacting case

The cryogenic LOX/GH2 mixing case described in the previous subsection is then ignited to further assess the performance of the proposed FPV model for reacting cases. The computational domain, mesh resolution, boundary conditions, and numerical schemes are kept the same as in the mixing case. The same flamelet table is adopted for the reacting case. As expected, large acoustic waves are generated right after the ignition takes places, and first-order upwinding scheme is used in the sponge layer to ensure that the pressure waves exit the domain and eventually a stable flame is obtained.

Figure 7 shows instantaneous results for temperature, mass fractions of H2, O2, and OH. The black curve in the temperature subfigure in Fig. 7 indicates the stoichiometric value of mixture fraction. Pink curves in Fig. 7 correspond to the PBP location which is characterized by the peak of the specific heat capacity. As can be seen from Fig. 7, the proposed FPV model is able to predict the transcritical flame robustly due to the double-flux model and entropy flux correction technique adopted in the numerical solver. No spurious oscillations in pressure or velocity were observed. Due to the adiabatic boundary conditions applied at the injector lip, the flame is attached. The flame exhibits laminar behavior close to the injector tip and is more wrinkled and turbulent downstream. The structure of a transcritical flame is found to be similar to that under ideal-gas conditions through flamelet studies. The PBP region is spatially separated from the reaction zone and the real-fluid pseudo-boiling process is found to take place in the region characterized by almost pure oxygen [66, 16]. This can be seen from the current simulation results that the PBP location as indicated by the pink curve has no interaction with the flame. The transcritical process happens to almost pure oxygen as can be seen from the oxygen mass fraction results in Fig. 7.

Figure 8 compares instantaneous density fields from both the cryogenic mixing and reacting cases. The legends are in logarithmic scale to emphasize the small structures at relatively low density values. As can be clearly seen from Fig. 8, the density fields show considerable differences between the two cases. Lower density can be seen in the reacting case due to the flame in between the two streams. The LOX stream in the reacting case shows suppressed vortical structures. Whereas in the mixing case, large vortical structures can be seen and the dense LOX stream penetrates into the GH2 stream by a height of more than $h$ above the injector lip. The heat release from the flame not only separates the two streams which results in an almost pure species pesudo-boiling process, but also suppresses the development of the vortical structures in the mixing layer. This is in consistent with the simulation results by Ruiz [67].

The proposed FPV model has been demonstrated to have the capability for transcritical reacting flow simulations. The current numerical scheme is directly applicable to LES of real applications under transcritical conditions.

6 Conclusions

A FPV approach is developed for trans- and supercritical combustion simulations in the context of finite volume, fully compressible, explicit solvers. The PR cubic EoS is used for the consistent evaluation of the thermodynamic quantities under transcritical conditions. Capabilities to calculate transcritical flamelet solutions have been implemented in the FlameMaster [57] solver and validated against DNS results. The double-flux model developed for transcritical flows [37] is used to eliminate spurious pressure oscillations caused by the nonlinearity inherent in the real-fluid EoS. A hybrid scheme with entropy-stable flux correction technique [37] is used to deal with the large density ratio in transcritical combustion cases. For the FPV approach, parameters $a$ and $b$ in the cubic EoS are pre-tabulated for the evaluation of the departure functions and a quadratic model is used to recover the attraction parameter $a$. The ideal-gas values are calculated from a linearized specific heat ratio model. [60] The novelty of the proposed approach lies in the ability to account for pressure and temperature variations from the reference tabulated values using a computationally tractable pre-tabulated combustion chemistry in a thermodynamically consistent fashion. An a priori analysis was conducted to assess the performance of the proposed transcritical FPV approach and it is shown that the assumption that the composition is a weak function of the reference conditions is valid under transcritical conditions, and the currently developed FPV approach can accurately recover the real-fluid and ideal-gas thermodynamic quantities despite perturbations in temperature and pressure. The solution of the laminar flamelets in mixture fraction space and the chemistry tabulation requires special considerations in order to account for the full non-linear effects in transcritical flows. The transcritical FPV approach works robustly and accurately even with an insufficient resolution in mixture fraction in the table. Cryogenic LOX/GH2 mixing and reacting cases are performed to demonstrate the capability of the proposed approach in multidimensional simulations. The current combustion model and numerical schemes are directly applicable to LES of real applications under transcritical conditions.

Acknowledgments

Financial support through NASA with award numbers NNX14CM43P and NNM13AA11G are gratefully acknowledged. The authors would like to thank Dr. Guilhem Lacaze for sharing data for comparison.