Framework for real-gas compressible reacting flows with tabulated thermochemistry

J. of Supercritical Fluids 101 (2015) 1–16

Contents lists available at ScienceDirect

The Journal of Supercritical Fluids journal homepage: www.elsevier.com/locate/supflu

Framework for real-gas compressible reacting flows with tabulated thermochemistry Xavier Petit, Guillaume Ribert ∗ , Pascale Domingo CORIA – UMR6614 CNRS, Normandie Université, INSA and Université de Rouen, Av. de l’Université, 76800 Saint-Etienne-du-Rouvray, France

a r t i c l e

i n f o

Article history: Received 4 November 2014 Received in revised form 16 February 2015 Accepted 16 February 2015 Available online 25 February 2015 Keywords: Real gas Compressible flow Tabulated chemistry

a b s t r a c t Combustion in liquid rocket engines happens under severe thermodynamical conditions: pressure exceeds the critical pressure of injected propellants and temperature is cryogenic. Such a situation requires an important effort of modeling: real gas effects are incorporated through cubic equations of state along with pressure-correction terms, and transport properties follow specific rules. Modeling for turbulent combustion is also an issue and is presently considered with tabulated chemistry thus reducing the number of transported variables. In this study, a framework is provided to deal with real-gas compressible reacting flows with tabulated thermochemistry. As a consequence, a cubic equation of state for tabulated thermochemistry is derived with the adopted thermodynamic relations, and the temperature computation is adapted to incorporate real gas effects through the tabulated thermochemistry approach. Two-dimensional reactive and non-reactive theoretical test cases have been performed with success to demonstrate the capacity of the new method. Finally, the simulation of a three-dimensional non-reactive single injector derived from Mayer’s experiment, that consists of a supercritical nitrogen jet injection into a warm nitrogen atmosphere, is also performed with success. The comparison between the tabulated approach and a fully coupled reference solution leads to similar results for the density values along the jet axis and for the spreading rate. © 2015 Elsevier B.V. All rights reserved.

1. Introduction Many energetic systems operate at elevated pressures in order to improve the global system efficiency. For aero-engines, the pressure in the combustion chamber is between 20 and 30 bar, but it reaches 100 bar and above for rocket engines. In the last example, the pressure exceeds the thermodynamic critical pressure, pc , of the injected propellants (typically H2 and O2 for the current Ariane 5 launcher) meaning that the fluid is in a supercritical state, i.e. there is no distinction between gas or liquid. Pressure effects have a direct impact on the Reynolds number through the kinematic viscosity: the Reynolds number increases with pressure involving a decrease of the Kolmogorov scale with consequences on the mesh grid requirements [1]. In addition, for flows under these severe thermodynamical conditions, the flame becomes very thin. Ribert et al. [2] or Pons et al. [3] have shown that diffusion laminar flames are micrometric, meaning that direct numerical simulations (DNS) are definitively out of reach even for massively parallel computers. A way to reduce the numerical CPU

∗ Corresponding author at: CORIA – UMR6614 CNRS. Tel.: +33 232959792. E-mail address: [email protected] (G. Ribert). http://dx.doi.org/10.1016/j.supflu.2015.02.017 0896-8446/© 2015 Elsevier B.V. All rights reserved.

cost is to use tabulated chemistry [4–12]. With this technique, thermo-chemical quantities ( l ) such as species mass fractions or reaction rates are tabulated as functions of a limited set of coordinates that describes physical phenomena like the chemical progress variable, mixture fraction, enthalpy, etc. Look-up tables are generated from canonical problems using complex chemical schemes solved under the low-Mach number assumption. Only few variables are then required to parametrize the whole chemical response compared to the detailed chemistry that needs a large number of species (Ns ) and elementary chemical reactions. The set of coordinates,  l , is transported along with mass, momentum and energy equations. Introducing the tabulated chemistry into large eddy simulations (LES) has already been done with the ideal gas assumption [13–16] and the recent work of Vicquelin et al. [16] raises the issues and proposes solutions to address such coupling when using a full compressible code: temperature computation and boundary conditions treatment appear as crucial elements that must be carefully validated. The first attempt to consider real gas effects through this procedure has been recently performed by Lacaze and Oefelein [17]. Their ultimate goal is to develop a combustion model that reproduces the flame behavior in a rocket engine. Equations for mass, momentum, energy and mixture fraction are considered. A model is introduced into the energy equation

2

X. Petit et al. / J. of Supercritical Fluids 101 (2015) 1–16

to recover the correct temperature when using the classical cubic equation of state for real gases. The whole procedure is validated on two-dimensional laminar counterflow diffusion flames that mimic the work of Ribert et al. [2] performed with a one-dimensional lowMach number code. However, several issues must still be addressed when a more realistic configuration that features real gas effects is simulated: as shown by Vicquelin et al. [16] the validation of the boundary conditions treatment is important for fully compressible codes. It becomes essential when real gas effects are considered because they may involve large density gradients that are challenging to capture [1,18]. Moreover the coupling between the equation of state (EoS) and the thermochemistry tabulated values must still be provided. Indeed, for ideal gases, the sole mean molecular weight of a mixture (W) is depending on the tabulated values (tab ) of the species mass fractions Yktab . It is then convenient

to directly tabulate Wtab . The equation of state is then simply cast as p = RT/W tab with  for density, T for temperature and p for pressure. R is the ideal gas constant. Defining such new EoS for real gases based on a very small set of tabulated data has not been proposed yet and, as it will be shown in the next section, is more challenging. In this study, a general framework is proposed to deal with the entire range of fluid thermodynamic states when a compressible code is coupled to a  l -dimensions look-up table for thermochemistry. In the next section, the formulation proposed by Vicquelin et al. is extended to a real gas formulation, and a new EoS is proposed to ensure the consistency between the stored information and the values of transported variables. Finally, in Section 3, reacting and non-reacting theoretical test cases are performed as well as the simulation of a non-reacting nitrogen injector under supercritical conditions to assess the whole methodology. 2. Framework for general fluid flows using a tabulated thermochemistry As stated in Section 1, coupling a compressible solver with a tabulated thermochemistry that features real gas effects requires specific ingredients: a precise temperature computation, an equation of state that can manage together real gas features and data coming from the look-up table, and a robust treatment of boundary conditions. A complete derivation of the last point can be found in [19]: the mathematical development of needed matrixes follows the usual procedure [20], but the partial derivatives must be linked to the considered EoS as shown hereafter. 2.1. Mathematical formulation

 tab

ij = 

∂ui ∂uj + ∂xj ∂xi

(1)

where U = [, ui , et ,  l ]T is the vector of conservative variables (unless stated otherwise, the tensor notation for repeated indices is adopted). ui are the velocity components, et is the total energy decomposed as the sum of the internal (e = h − p) and kinetic (ec = u2i /2) energies. h is the enthalpy of the Ns species mixture.  l is the lth coordinate used to describe the thermochemical j database. Fc represents the conservative variables flux vectors and j j T is given by: Fc = [uj , uj ui + ıji p, uj et + uj p, uj l ] . Fd are the T

diffusive fluxes, Fd = [0, −ij , −ui ij + qj , −D(∂l /∂xj )] . Assuming a unitary Lewis number, the diffusion coefficient, D, is defined based on the tabulated conductivity, tab , and the heat capacity

−

2 tab ∂uk ıij ,  3 ∂xk

(2)

term vector S = (0, 0, 0, 0, 0, ω˙ l )T with ω˙ l the chemical rate of the variable  l . For general fluid flows, cubic equations of state (EoS) [21] replace the classical ideal gas assumption. These EoS have three parameters (Eq. (3)) which quantify the attractive forces among particles (a˛), the co-volume of the particles (b), and the deviation from spherical symmetry in a molecule (namely the acentric factor (ω)). p=

RT 2 a˛(ω, T ) . − 2 W − b W 2 + ubW + wb 2

(3)

In Eq. (3), R is the ideal gas constant and W, the molar mass of the mixture. (u, w) depend on the retained model for EoS: (u, w) = (0, 0) for Van der Waals, (u, w) = (1, 0) for Redlich–Kwong or Soave–Redlich–Kwong (SRK) and (u, w) = (2, −1) for Peng–Robinson EoS, respectively [21]. The compressibility factor, Z, is written as Z = p/(RT ). Transport properties are replaced by accurate high-pressure relations such as the one proposed by Chung et al. [22] for viscosity, , and thermal conductivity, . 2.2. Temperature computation Solving Eqs. (1)–(3) together requires to compute the temperature from the internal energy, e, and still be in agreement with the thermochemical database. However, the equations governing such database are usually isobaric and written under a low-Mach number assumption. They may not be exactly representative of compressible flows because acoustic waves locally generate energy fluctuations. As a consequence, the computed / etab , energy is then slightly different from the tabulated one, e = and a special treatment must be done [15,16]. For conciseness, any variable Qtab will stand for Q(l ). For fluid flows featuring real-gas effects, the sensible energy depends on T, Yk and p (or ). Assuming that the mixture composition is fixed and coming from the database (Yk = Yktab ), i.e. the k transport equations are replaced by the transport of  l variables, the energy departure between transported and tabulated energies is given by e = e(T, ) − etab (Ttab , tab ). The differential form of e(T, ) is written



∂e  ∂e  de = dT + d ∂T ,Y ∂ T,Y k

j ∂U ∂Fcj ∂Fd + = S, + ∂t ∂xj ∂xj



where tab is the tabulated dynamic viscosity. The heat diffusion flux, qi , is then expressed as qi = −tab /Cptab (∂h)/(∂xi ). The source



The governing equations for reactive flows can be written as:

j

at constant pressure, Cptab , as D = tab /(Cptab ). The viscous stress tensor  ij is written:

(4)

k

or



∂e  with Cv = ∂T ,Y

de = Cv dT + T d,



k

Using the integral form of Eq. (5),



T

e = e − etab =

e = e − etab



T tab

−



Cv (T ,  T0



Cv (T  , ) dT +

T0



tab

∂e  and T = . (5) ∂ T,Y

T

k

becomes

T (T,  ) d

T0 



T tab

) dT − T0

T (T tab ,  ) d .

(6)

X. Petit et al. / J. of Supercritical Fluids 101 (2015) 1–16

3

Table 1 [VORT] Initial thermodynamical conditions, p0 = 5.6 (MPa). Name

z (–)

c (–)

T (K)

Ux (m/s)

Uy (m/s)

M (–)

Re (–)

M1 M2 M3

0.1 0.1 0.1

0 0 0

146 146 146

60 300 500

0 0 0

0.09 0.45 0.75

4.5 × 106 22.5 × 106 37.5 × 106

In Eq. (5), T can be re-arranged as:

 ∂e  

T = Yk ∂k  Ns

k=1

(7)

. T,j = / k

Assuming moderate1 fluctuations of pressure, T and Cv depend on tabulated density (tab ) and temperature (Ttab ) as: Cv (T  , ) = Cv (T  , tab ) and

T (T,  ) = T (T tab ,  ),

(8)

with tab and Ttab evolving with the table parameters defined by the transport of  l . In addition, Cv (resp. T ) is considered having constant values on a small temperature (density) range which is a reasonable hypothesis away from the critical point. e is now written:



T

Cv (T tab , tab ) dT

+ T tab

T (T tab , tab ) d = Cvtab [T − T tab ] + Ttab [ − tab ].

T = T tab +

[e − etab ] − T tab [ − tab ] Cvtab

.

(10)

RT 2 a˛ − , W1 W2

tab2

tab

with W1 = Wtab − btab and W2 = W tab + ub W tab  + wb 2 , to ensure that the data coming from the look-up table is also consider in the pressure computation. Moreover, a˛ = a˛(T, Ytab ) and a first order truncated Taylor expansion around T = Ttab is done, leading to:



∂a˛ ∂T

tab · (T − T tab ).

(12)

As a consequence, Eq. (11) becomes p=



l

tab2

(uW tab btab + 2wb

RT 2 a˛tab (∂a˛/∂T ) − − W1 W2 W2

tab

T

+

(∂a˛/∂T ) W2

tab tab T

,

(13)

with consequences on partial derivatives, i.e. with an impact on the speed of sound evaluation and the treatment of boundary conditions.

1 For substantial pressure fluctuations, a Taylor expansion of Cv [16] and T with temperature may be required.

tab

)a˛tab 2

tab

2 (ub

tab2

W tab + 2wb

2[(da˛)/(dT )] (T − T tab ) W2

−

W22

)[(da˛)/(dT )]

tab

(T − T tab )

W22

,

(15)



tab R 2 [(da˛)/(dT )] ∂p  = − ,  W1 W2 ∂T ,

(16)

l

and,





+

⎣

2



W2

 +

2 RT W12

 −

=

l= / k

⎡

(11) 2

(14)

dk .

l= / k

k

l

RT RTbtab 2a˛tab ∂p  = + − 2  W1 W2 ∂ T, W1

∂p  ∂l T,,

Eq. (3) is used when dealing with a multicomponent formulation. However, in the context of tabulated thermochemistry, variables such as the molecular weight or co-volume of molecules are coming from the look-up table. As a consequence, the EoS is changed as:



It becomes then straightforward to recover the desired partial derivative for pressure:

(9)

2.3. Equation of state for tabulated thermochemistry

a˛  a˛tab +



l

+

Finally, the temperature can be evaluated from tabulated values and the transport equations of e and :

p=



 ∂p  ∂p  ∂p   dp = d + dT + ∂ T, ∂T , ∂k T,,



T tab T

2

tab

(ub W tab + 2wbtab )d and the differential form of pressure with the general cubic EoS (Eq. (11)):

+

e =



Such partial derivatives are now expressed using the differential forms of dW1 = dWtab − dbtab − btab d and tab dW2 = (2W tab + btab u)dW tab + (W tab u + 2wb 2 )dbtab +

2 W2



2 a˛ W22

(2W tab + btab u) −

⎤  tab  tab   ∂a˛  ⎦ ∂T  ∂T  ∂k  Yk

+

2 a˛ W22

 ∂a˛tab  ∂k 

(W

tab 2

u + 2wb

 − l = / k

 )



∂W tab  ∂k 

l= / k



∂btab  ∂k 

 

2 (T − T tab ) W2



W12

l = / k

 tab

RT

l= / k

tab    (17). 

∂ (∂a˛)/(∂T ) ∂k

l = / k

These partial derivatives are required to express the Jacobean matrix M, that links conservative variables U to primitives variables V as M = ∂U/∂V, to treat the Navier-Stokes Characteristic Boundary Conditions, that formally have the same expression than with an ideal gas assumption, and to evaluate the speed of sound, a that will take the form of a2 = Cp /Cvtab × (∂p)/(∂)|T,l . More details may be found in [19]. 3. Results The above models have been implemented into the solver SiTCom-B that already features real gas effects [23–25] and NSCBC3D treatment [14,26,27]. The different types of inlet and outlet conditions (subsonic, supersonic, etc.) are managed with the same procedure than in [24]. It usually assumes that the flow is normal to the boundary plane (NSCBC-1D). However, Lodato et al. [26] recently extended this formulation to account for convection and

4

X. Petit et al. / J. of Supercritical Fluids 101 (2015) 1–16

pressure gradients in boundary planes, resulting in a 3D-NSCBC approach. The introduction of these additional transverse terms requires a specific treatment that is presently addressed in the case of flows featuring real gas effects. The parameter ˇ will be

introduced as a transverse damping parameter following [26]. In the rest of the study, only the SRK EoS is considered meaning that the coefficients u and w required in the general form of the cubic EoS are fixed to u = 1 and w = 0. Four test cases are considered to assess

Fig. 1. [VORT] Time evolution of a vortex convection (case M1) described by isolines that are plotted for |u1 − Ux | = 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4 and 4.5 m/s, with Ux = 60 m/s. The field of p* is also displayed. Left column: NSCBC-1D; right column: NSCBC-3D.

X. Petit et al. / J. of Supercritical Fluids 101 (2015) 1–16

the methodology described above. They are identified by acronyms given in square brackets and consider a mixture of liquid oxygen (LOx) at 120 K and methane at 288 K for a pressure of 5.6 MPa, with a high Reynolds number to be representative of typical supercritical test-benches:

5

• [VORT] is the 2D convection of a non-reacting vortex exiting the domain of calculation at a fixed composition of a LOx/CH4 mixture. • [MIXT] is the 2D convection of a non-reacting mixture fraction pocket and exiting the domain of calculation.

Fig. 2. [VORT] Time evolution of the relative error (Eq. (20)) at three locations (case M1): red for NSCBC-1D, black for NSCBC-3D. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

6

X. Petit et al. / J. of Supercritical Fluids 101 (2015) 1–16

• [REAC] is the 2D convection of a reacting fresh gas LOx/CH4 mixture into the corresponding burned gases and exiting the domain of calculation. • [JET] is the 3D non-reacting simulation of a cold nitrogen jet. The objectives of these test cases are to evaluate the amplitude of the generated noise when a perturbation crosses a boundary. Indeed, such noise must be as small as possible to avoid any contamination of the whole domain of simulation. In addition, when the perturbation leaves the domain of simulation, the background pressure must be recovered. The transported variables required to parametrize the tabulated chemistry are the mixture fraction z, progress variable c defined as eq eq c = Yc /Yc with Yc = YCO2 + YCO + YH2 O + YH2 and Yc the equilibrium response of the progress of reaction. The tabulated chemistry is built based on a library of laminar premixed flames [7,8] generated with the REGATH numerical code [2,3,28]. For the JET case, only the mixing of nitrogen between two thermodynamical states is tabulated. The following terms are then stored in a look-up table then used in the four test cases: W, T, Cv , e, ω˙ c , , , T , b, a˛ and d(a˛)/dT. The solver SiTCom-B that is used to compute the following test cases features a 4th-order skew-symmetric like numerical scheme [29] for the spatial integration, and a 4th order Runge-Kutta method for the time integration. A numerical procedure adds artificial viscosity [30,31] to zones where a sensor [24] detects either strong density or strong pressure gradients, thus damping spurious numerical oscillations. 3.1. Two-dimensional vortex convection [VORT] A two-dimensional domain (Lx × Ly = 2 ×1 cm2 ) is simulated with 800 × 400 homogeneous mesh cells. Inlet and outlet are treated with the real gas NSCBC (x-axis). Periodic conditions are used for faces perpendicular to the y-axis. Three simulations, called M1, M2 and M3, have been performed with the parameters given in Table 1. Mach (M) and Reynolds (Re) numbers are based on the convective velocity and the vortex length. The pressure is set to p = 5.6 MPa. z = 0.1 and c = 0 correspond to a mixture of oxygen and methane (YCH4 = 0.1 and YO2 = 0.9) with a temperature of 146 K. The compressibility factor is around Z = 0.18 meaning that real gas effects cannot be neglected. A velocity potential (˚) is used to initialize the velocity field: Ux (x, y) =

∂˚ ∂y

and Uy (x, y) = −

with ˚=

√

 2 Uvort

exp

∂˚ , ∂x

(x − x0 )2 + (y − y0 )2 − 2r 2

(18)

 + Ux y + Uy x.

(19)

Ux and Uy are the two components of the convective velocity along x and y axis. Ux is set to 60, 300 and 600 m/s for the simulations M1, M2 and M3, respectively. The vortex magnitude (Uvort ) is fixed to 5 m/s for the three simulations. r = 0.1 cm and represents a vortex radius. x0 and y0 are the coordinates of the vortex center. The initial pressure and temperature fields were not modified to match with the initial flow modifications. However, it did not generate any strong perturbations at the beginning of the simulation. The three simulations have been performed using either the NSCBC-1D (ˇ = 1.0) or NSCBC-3D (ˇ ∈ [0, 1[) formalisms. For the latter, the best value has been found for ˇ = M where the Mach number M is coming from the target configuration (see Table 1). The coefficient for the pressure relaxation term is set to p = 2 .105 s−1 . In Fig. 1 the vortex, which is described with iso-lines of streamwise velocity, is leaving the domain of calculation on the right.

Fig. 3. [MIXT] Initial flow field of mixture fraction.

The pressure field is described with the relative error, p* , that is defined as p∗ (x, y, t) =

p(x, y, t) − p0 , p0

(20)

where p0 = 5.6 MPa is the target pressure, i.e. once the vortex left the domain of calculation p* → 0. Using the NSCBC-1D formalism leads to a deformation of the vortex iso-lines, an increase of the vortex length and |p* | reaches 5%. Using the NSCBC-3D formalism keeps the vortex shape and magnitude; |p* |< 0.3 %. In Fig. 2 the streamwise normalized profiles of pressure are plotted at different locations (y = −0.1, y = 0.0 and y = +0.1 cm) and a comparison between the two NSCBC formulations is performed with the configuration M1. Very good results have been obtained when running the cases M2 and M3 as shown in [19] with a value of |p* |< 0.3 % found with the NSCBC-3D treatment meaning that the exit of the vortex does not modify the incoming flow. 3.2. Two-dimensional convection of a mixture fraction pocket [MIXT] A two-dimensional domain (Lx × Ly = 2 ×1 cm2 ) is simulated with 800 × 400 homogeneous mesh cells. Inlet and outlet are treated with the real gas NSCBC (x-axis). Periodic conditions are used for faces in the y-axis. The initial field of mixture fraction z is computed with p0 = 5.6 MPa, Ux = 60 m/s, Uy = 0 m/s, c = 0 and z = zm +

z0 − zm [1 + tanh(Cr (r − R0 ))], 2



(21)

with r = (x − x0 )2 + (y − y0 )2 . r represents the distance from the pocket center whose coordinates are (x0 , y0 ) = (0, 0). R0 is the pocket radius and the constant Cr = 2.5 103 m−1 controls the profile stiffness of the hyperbolic tangent. The reference mixture fraction is z0 = 0.1. As the pressure is set to p = 5.6 MPa, the present simulation deals with a mixture of oxygen and methane (YCH4 = 0.1 and YO2 = 0.9) at 146 K, i.e. Z = 0.18. zm is fixed to 0.3 and represents the maximum mixture fraction found at the center of the pocket of mixture. zm = 0.3 corresponds to a mixture of oxygen and methane (YCH4 = 0.3 and YO2 = 0.7) at 167 K, i.e. with a compressibility factor Z = 0.3. The initial fields are shown in Figs. 3 for mixture fraction. In this simulation, the Reynolds number is equal to Re = 4.5 106 and the Mach number is M = 0.09. In Fig. 4 the density and reduced pressure fields (p* ) are plotted for different times of the simulation with NSCBC-1D. No deformation is observed on the density isolines (750 kg/m3 , 700 kg/m3 , 650 kg/m3 , 600 kg/m3 , 550 kg/m3 , 520 kg/m3 , 500 kg/m3 and 490 kg/m3 ) that describe the fluid pocket. Using a NSCBC-3D treatment leads to similar results. This is confirmed in Figs. 5 and 6 where the streamwise profiles of  and p* are compared with both types of boundary conditions for different times of simulation.

X. Petit et al. / J. of Supercritical Fluids 101 (2015) 1–16

7

Fig. 4. [MIXT] Snapshots of density (left) and pressure (p* on the right) fields with NSCBC-1D. Isolines of density describe a pocket of light fluid exiting the domain of calculation (750 kg/m3 , 700 kg/m3 , 650 kg/m3 , 600 kg/m3 , 550 kg/m3 , 520 kg/m3 , 500 kg/m3 and 490 kg/m3 ).

3.3. Two-dimensional convection of a reacting fresh gases pocket [REAC] A two-dimensional domain (Lx × Ly = 2 ×1 cm2 ) is simulated with 800 × 400 homogeneous mesh cells. Inlet and outlet are

treated with the real gas NSCBC (x-axis). Periodic conditions are used for faces in the y-axis. The initial field is computed with p0 = 5.6 MPa, Ux = 60 m/s, Uy = 0 m/s, z = 0.5 and c = cm +

c0 − cm [1 + tanh(Cr (r − R0 ))], 2

(22)

8

X. Petit et al. / J. of Supercritical Fluids 101 (2015) 1–16

Fig. 5. [MIXT] Centerline density profile for six times of calculation with NSCBC-1D and NSCBC-3D treatment.



with r = (x − x0 )2 + (y − y0 )2 . r represents the distance from the pocket center whose coordinates are (x, y) = (0, 0). R0 is the pocket radius and Cr is set to 2.5 103 m−1 . cm is fixed to 0 and represents the minimum progress variable found in the simulation, i.e. inside the fresh gases pocket. This configuration corresponds to a mixture

of oxygen and methane (YCH 4 = 0.5 and YO2 = 0.5) at 219 K, i.e. Z = 0.77. The reference progress variable encountered in the burned gases is c0 = 1 and corresponds to the combustion of the fresh gases pocket. Z = 1 and the adiabatic flame temperature is Tad = 1450 K

X. Petit et al. / J. of Supercritical Fluids 101 (2015) 1–16

9

Fig. 6. [MIXT] Centerline normalized pressure (p* ) profile for six times of calculation with NSCBC-1D and NSCBC-3D treatment.

when using Lindstedt’s mechanism [32]. In Fig. 7, the temperature of the tabulated premixed flamelet is plotted as a function of the progress variable. The temperature increases till a maximum close to Tmax = 1880 K before decreasing to Tad , that is a classical behavior when dealing with oxy-combustion. The initial fields are shown in Fig. 8 for c and source term of c. In this

simulation, the Reynolds number is equal to Re = 16 × 103 and the Mach number is M = 0.06. In Fig. 9, the fields of source term of c are plotted for different times of calculation with NSCBC-1D (ˇ = 1.0) and NSCBC-3D (ˇ = 0.625) treatments. The pocket of fresh gases exits the domain of calculation without any spurious waves for both cases.

10

X. Petit et al. / J. of Supercritical Fluids 101 (2015) 1–16

Fig. 7. [REAC] Premixed laminar flamelet.

This is confirmed in Figs. 10 and 11, where the time evolution of the progress variable c and reduced pressure (p* ), are plotted along the x-axis, respectively. The hyperbolic tangent (Eq. (22)) used for the initial profile of c evolves with time, exhibiting a local maximum, because the flame front has its own speed of propagation (flame speed [20]) that is oriented outward the pocket. As a consequence the temperature follow the trend and almost reach Tmax . The reduced pressure never exceeds 0.6% meaning that the pressure level is maintained. 3.4. Three-dimensional large-eddy simulation of a supercritical single jet [JET] The topic addressed in this test case is inspired by the work of Mayer et al.’s experiment [33]. A cold nitrogen (pc = 3.34 MPa and Tc = 126.2 K) is injected into a warm nitrogen environment (298 K) and the ambient pressure in the chamber is 3.98 MPa. The diameter of the injector is dinj = 2R = 2.2 mm and the chamber length is 1000 mm with a diameter of 122 mm. Velocity, temperature, density and speed of sound at the inlet boundary are 5.4 m/s, 137 K, 163.5 kg/m3 and 202 m/s, respectively. This configuration is modified as in [24] with an inlet velocity of 50 m/s, all other parameters being the same as in the original experiment. The compressibility factor, Z, for inlet conditions is around 0.6 but close to unity for ambient fluid. 3D simulations are performed for a domain of calculation defined by Lx × Ly × Lz = 85 × 60 × 60 mm. The mesh

arrangement corresponds to Nx × Ny × Nz = 464 × 280 × 280 ≈ 36.4 millions of points. Only a small part of the experimental chamber is simulated to save the time of calculation, but the boundary conditions are adapted to still have a representative simulation of the configuration. Real-gas NSCBC treatment is applied with fixed velocity and  l values at inlet, and prescribed pressure for outlets. Walls are considered as adiabatic walls. A turbulent pipe profile as in [24] is prescribed at inlet. The dynamic Smagorinsky model [34] is used to evaluate the turbulent viscosity term. For this nonreacting case, the only control parameter of the look-up table is enthalpy, h, as used in [35,36] for example. It describes a mixing of N2 at 137 K (hinlet or z = 1) with N2 at 298 K (hrest or z = 0). The mixing variable defined as z = (h − hrest )/(hinlet − hrest ). A simulation is also performed without resorting to a tabulation of thermodynamics to compare both approaches. In Fig. 12(a), the instantaneous flow is viewed through the density field. High density gradients are located along a potential core that starts from the injection till 20dinj . The jet enters into the cavity, quickly destabilizes and opens as a part of the light surrounding nitrogen penetrates into the main core jet. In Fig. 12(b), the maximum velocity is found in the injection plane because of the inlet turbulent profile. In Fig. 12(c) the mixing between the dense and light fluids leads to the creation of dense pockets that seem to persist in the main core jet with high speed velocity. These results are fully similar to those found in [24] for which a fully coupled

Fig. 8. [REAC] Initial flow field: (a) progress variable and (b) source term of the progress variable.

X. Petit et al. / J. of Supercritical Fluids 101 (2015) 1–16

11

Fig. 9. [REAC] Snapshots of source term of c (left) and pressure (p* on the right) fields with NSCBC-3D. Isolines of (c = 0.05, 0.1, 0.2, 0.3, or 0.95) describe a pocket of fresh fluid exiting the domain of calculation filled of burned gases.

formulation is used. This is confirmed in Fig. 13(a) for the normalized density evolution (+ = [ − rest ]/[inlet − rest ]) on the jet axis, and in Fig. 13(b) for the jet opening. This latter is based on density evaluated using the half width, half maximum (HWHM)

indicator [24,37,38]. Cross-stream profiles of average density are used to define the HWHM length, L , of the jet, taken at the median value between dense and light regions along the x-axis, meaning that L /D = 0.5 as L = D/2 at the injection plane.

12

X. Petit et al. / J. of Supercritical Fluids 101 (2015) 1–16

Fig. 10. [REAC] Time evolution of the progress variable at y = 0 with NSCBC-1D and NSCBC-3D treatments for a pocket of fresh fluid exiting the domain of calculation filled of burned gases.

X. Petit et al. / J. of Supercritical Fluids 101 (2015) 1–16

13

Fig. 11. [REAC] Time evolution of the relative pressure (p* ) at y = 0 with NSCBC-1D and NSCBC-3D treatments for a pocket of fresh fluid exiting the domain of calculation filled of burned gases.

14

X. Petit et al. / J. of Supercritical Fluids 101 (2015) 1–16

Fig. 12. Mayer et al.’s supercritical jet configuration: instantaneous flow field.

X. Petit et al. / J. of Supercritical Fluids 101 (2015) 1–16

15

(1) small impact of acoustic waves on the flow when a perturbation crosses a boundary, and (2) convergence back to the initial state when the perturbation leaves the domain of simulation. Finally, the comparison between the tabulated and fully coupled approaches on a supercritical jet highlights the ability of the proposed methodology to handle complex flows. Future studies will then be focused on reacting configurations where the development of a turbulent combustion model will be a major issue.

Acknowledgments Xavier Petit was financially supported by SAFRAN Snecma, Space Engines Division, and CNES Launchers Directorate. Computations were performed using HPC resources from GENCI under the allocation 2013-020152 and CRIHAN.

References

Fig. 13. [JET] Comparison between tabulated chemistry (—) and multi-species (symbols) formulations for (a) density values on jet axis and (b) spreading rate.

4. Conclusion Incorporating a tabulated thermochemistry into a fully compressible solver must be performed with care. Indeed, the tabulated values needs a correction because of the acoustic waves that are not accounted for in the look-up table. In this study, the emphasis is put on situations where real gas effects occur. With an established formulation, a cubic EoS usually replaces the ideal gas EoS, and transport and thermodynamic properties are corrected with pressure-dependent terms. With a tabulated approach, such a cubic EoS must be adapted to incorporate the data coming from the precomputed look-up table. This issue has been presently addressed with the addition of only four additional stored variables ( T , b, a˛ and da˛/dT) compared to an ideal gas formulation with tabulation. T is the partial derivative of the internal energy with respect to density for a fixed temperature and composition. T plays a key role in the computation of temperature when a compressible solver uses tabulated thermochemistry featuring real gas effects. Indeed,

T acts as a correction term of the tabulated temperature in order to account for the impact of acoustic waves. This correction is needed since the internal energy depends on temperature and pressure (or density) in contrast to the ideal gas approach for which only the temperature is considered. Neglecting this correction leads to an erroneous value of temperature. To validate the proposed developments and their integrations into the derivation of boundary conditions, three customary test-cases have been simulated with success: convection of a 2D vortex into a uniform gas composition, a non-reacting pocket of mixture fraction, and a reacting fresh gas pocket into burned gases. They all fill out the requirements that are

[1] V. Yang, Modeling of supercritical vaporization, mixing, and combustion processes in liquid-fueled propulsion systems, Proceedings of the Combustion Institute 28 (2000) 925–942. [2] G. Ribert, N. Zong, V. Yang, L. Pons, N. Darabiha, S. Candel, Counterflow diffusion flames of general fluids: oxygen/hydrogen mixtures, Combustion and Flame 154 (2008) 319–330. [3] L. Pons, N. Darabiha, S. Candel, G. Ribert, V. Yang, Mass transfer and combustion in transcritical non-premixed counterflows, Combustion Theory and Modeling 13 (2009) 57–81. [4] N. Peters, Laminar diffusion flamelet models in non-premixed turbulent combustion, Progress in Energy Combustion Science 16 (1984) 319–339. [5] D. Bradley, L. Kwa, A. Lau, M. Missaghi, S. Chin, Laminar flamelet modeling of recirculating premixed methane and propane-air combustion, Combustion and Flame 71 (1988) 109–122. [6] U. Maas, S. Pope, Simplifying chemical kinetics: intrinsic low-dimensional manifolds in composition space, Combustion and Flame 88 (1992) 239–264. [7] O. Gicquel, N. Darabiha, D. Thevenin, Laminar premixed hydrogen/air counterflow flame simulations using flame prolongation of ILDM with differential diffusion, Proceedings of the Combustion Institute 28 (2000) 1901–1908. [8] J. van Oijen, F. Lammers, L. de Goey, Modeling of premixed laminar flames using flamelet generated manifolds, Combustion Science and Technology 127 (2001) 2124–2134. [9] G. Ribert, O. Gicquel, N. Darabiha, D. Veynante, Tabulation of complex chemistry based on self-similar behaviour of laminar premixed flames, Combustion and Flame 146 (2006) 649–664. [10] K. Wang, G. Ribert, P. Domingo, L. Vervisch, Self-similar behavior and chemistry tabulation of burnt-gas diluted premixed flamelets including heat-loss, Combustion Theory and Modeling 14 (2010) 541–570. [11] G. Ribert, K. Wang, L. Vervisch, A multi-zone self-similar chemistry tabulation with application to auto-ignition including cool-flames effects, Fuel 91 (2012) 87–92. [12] G. Ribert, L. Vervisch, P. Domingo, Y.-S. Niu, Hybrid transported-tabulated strategy to downsize detailed chemistry for numerical simulation of flames, Flow Turbulence and Combustion 92 (2014) 175–200. [13] A. Kempf, H. Forkel, J.-Y. Chen, A. Sadiki, J. Janicka, Large-eddy simulation of a counterflow configuration with and without combustion, Proceedings of the Combustion Institute 28 (2000) 35–40. [14] P. Domingo, L. Vervisch, D. Veynante, Large-eddy simulation of a lifted methane jet flame in a vitiated coflow, Combustion and Flame 152 (2008) 415–432. [15] J. Galpin, A. Naudin, L. Vervisch, C. Angelberger, O. Colin, P. Domingo, Largeeddy simulation of a fuel-lean premixed turbulent swirl-burner, Combustion and Flame 155 (1–2) (2008) 247–266. [16] R. Vicquelin, B. Fiorina, S. Payet, N. Darabiha, O. Gicquel, Coupling tabulated chemistry with compressible CFD solvers, Proceedings of the Combustion Institute 33 (1) (2011) 1481–1488. [17] G. Lacaze, J. Oefelein, A non-premixed combustion model based on flame structure analysis at supercritical pressures, Combustion and Flame 159 (2012) 2087–2103. [18] J. Oefelein, V. Yang, Modeling high-pressure mixing and combustion processes in liquid rocket engines, J. Propulsion and Power 14 (1998) 843–857. [19] X. Petit, Turbulence–chemistry interaction in a LOx–CH4 flame (Ph.D. thesis), INSA de Rouen, 2014, pp. 150–169 (in French). [20] T. Poinsot, D. Veynante, Theoretical and Numerical Combustion, 3rd ed., 2012. [21] B. Poling, J. Prausnitz, J. O’Connell, The Properties of Gases and Liquids, 5th ed., McGraw Hill, 2001. [22] T. Chung, M. Ajlan, L. Lee, K. Starling, Generalized multiparameter corresponding state correlation for polyatomic, polar fluid transport properties, Industrial and Engineering Chemistry Research 27 (1988) 671–679. [23] G. Ribert, D. Taieb, X. Petit, G. Lartigue, P. Domingo, Simulation of supercritical flows in rocket-motor engines: application to cooling channel and injection system, EUCASS Book Series 4 (2012) 205–226.

16

X. Petit et al. / J. of Supercritical Fluids 101 (2015) 1–16

[24] X. Petit, G. Ribert, P. Domingo, G. Lartigue, Large-eddy simulation of supercritical fluid injection, J. Supercritical Fluids 84 (2013) 61–73. [25] http://www.coria-cfd.fr/index.php/SiTCom-B, 2014. [26] G. Lodato, L. Vervisch, P. Domingo, Three-dimensional boundary conditions for direct and large-eddy simulation of compressible viscous flows, J. Computational Physics 227 (2008) 5105–5143. [27] P. Domingo, L. Vervisch, S. Payet, R. Hauguel, DNS of a premixed turbulent V flame and LES of a ducted flame using a FSD-PDF subgrid scale closure with FPI-tabulated chemistry, Combustion Flame 143 (2005) 566–586. [28] S. Candel, T. Schmitt, N. Darabiha, Progress in transcritical combustion: Experimentation, modeling and simulation, in: 23rd ICDERS Conf., 2011. [29] F. Ducros, F. Laporte, T. Souleres, V. Guinot, P. Moinat, B. Caruelle, Highorder fluxes for conservative skew-symmetric-like schemes in structured meshes: application to compressible flows, J. Computational Physics 161 (2000) 114–139. [30] A. Jameson, W. Schmidt, E. Turkel, Numerical solution of the Euler equations by finite volume method using Runge-Kutta time stepping schemes, in: AIAA Paper 1259:1281, 1981. [31] S. Tatsumi, L. Martinelli, A. Jameson, Flux-limited schemes for the compressible Navier–Stokes equations, AIAA J. 33 (2) (1995) 252–261.

[32] W. Juchmann, H. Latzel, D. Shin, G. Pieter, T. Dreier, H. Volpp, R. Lindstedt, K. Leung, Absolute radical concentration measurements and modeling of lowpressure CH4 /O2 /NO flames, Proceedings of the of the Combustion Institute 27 (1998) 469–476. [33] W. Mayer, J. Telaar, R. Branam, G. Schneider, J. Hussong, Raman measurements of cryogenic injection at supercritical pressure, Heat and Mass Transfer 39 (2003) 709–719. [34] J. Smagorinsky, General circulation experiments with the primitive equations, Monthly Weather Review 91 (1963) 99–164. [35] J. van Oijen, L. de Goey, Modelling of premixed laminar flames using flamelet-generated-manifolds, Combustion Science and Technology 161 (2000) 113–137. [36] G. Ribert, M. Champion, O. Gicquel, N. Darabiha, D. Veynante, Modeling nonadiabatic turbulent premixed reactive flows including tabulated chemistry, Combustion and Flame 141 (2005) 271–280. [37] M. Oschwald, A. Schik, Supercritical nitrogen free jet investigated by spontaneous Raman scattering, Experiments in Fluids 27 (1999) 497–506. [38] R. Branam, W. Mayer, Characterization of cryogenic injection at supercritical pressure, J. Propulsion and Power 19 (2003) 342–355.