Xi Deng
email: deng.xi98@gmail.com

Dr Pierre Boivin
email: pierre.boivin@univ-amu.fr

Diffuse interface modelling of reactive multi-phase flows applied to a sub-critical cryogenic jet

Keywords: multiphase flows, reacting flows, phase transition, diffuse interface

In order to simulate cryogenic H 2 -O 2 jets under subcritical condition, a numerical model is constructed to solve compressible reactive multi-component flows which involve complex multi-physics processes such as moving material interfaces, shock waves, phase transition and combustion. The liquid and reactive gaseous mixture are described by a homogeneous mixture model with diffusion transport for heat, momentum and species. A hybrid thermodynamic closure strategy is proposed to construct an equation of state (EOS) for the mixture. The phase transition process is modeled by a recent fast relaxation method which gradually reaches the thermo-chemical equilibrium without iterative process. A simplified transport model is also implemented to ensure the accurate behavior in the limit of pure fluids and maintain computational efficiency. Last, a 12-step chemistry model is included to account for hydrogen combustion. Then the developed numerical model is solved with the finite volume method where a low dissipation AUSM (advection upstream splitting method) Riemann solver is extended for multi-component flows. A homogeneous reconstruction strategy compatible with the homogeneous mixture model is adopted to prevent numerical oscillations across material interfaces.

Having included these elements, the model is validated on a number of canonical configurations, first for multiphase flows, and second for reactive flows. These tests allow recovery of the expected behavior in both the multiphase and reactive limits, and the model capability is further demonstrated on a 2D burning cryogenic H 2 -O 2 jet, in a configuration reminiscent of rocket engine ignition.

Introduction

This study aims at providing a numerical framework for the simulation of reactive flows with distinct liquid and gas phases, as observed in cryogenic rocket engines under subcrictical conditions. In such flows, there are so many physical processes involved that it becomes highly challenging to address all of them with a consistent degree of accuracy. It is therefore of interest to develop a numerical framework with the following characteristics: (i) the framework developed should be compatible with existing gaseous combustion models, including detailed kinetic descriptions, as well as simple turbulent combustion models; (ii) the framework shall include a multi-phase thermodynamic closure, and be fully compatible with phase transition modelling; (iii) the framework should allow the presence of complex interfaces, including their dynamic creation and destruction. There should not be a limitation on density and momentum ratios across the interface.

Although the literature teems with models encompassing two of the three characteristics listed above, very few [START_REF] Gaillard | A diffuse interface lox/hydrogen transcritical flame model[END_REF][START_REF] Gaillard | Interfaces diffuses et flammes transcritiques lox/h2[END_REF] allow the simultaneous description of all three. Thus accurate description of reactive flows including multiple phases is still a challenge to numerical modelling. For instance, most studies consider the presence of liquid(s) in reactive flows in either of the following conditions: (i) The liquid phase is already in disperse form. This is made possible by studies on spray combustion modelling [START_REF] Williams | Spray combustion and atomization[END_REF][START_REF] Reveillon | Analysis of weakly turbulent dilute-spray flames and spray combustion regimes[END_REF][START_REF] Sánchez | The role of separation of scales in the description of spray combustion[END_REF]. (ii) The liquid phase is dense, but close to or above critical pressure [START_REF] Singla | Transcritical oxygen/transcritical or supercritical methane combustion[END_REF][START_REF] Oefelein | Mixing and combustion of cryogenic oxygen-hydrogen shear-coaxial jet flames at supercritical pressure[END_REF][START_REF] Schmitt | Large-eddy simulation of oxygen/methane flames under transcritical conditions[END_REF][START_REF] Ruiz | Unsteady numerical simulations of transcritical turbulent combustion in liquid rocket engines[END_REF][START_REF] Urbano | Study of flame response to transverse acoustic modes from the les of a 42-injector rocket engine[END_REF]. In that case, a cubic equation of state [START_REF] Peng | A new two-constant equation of state[END_REF][START_REF] Soave | Equilibrium constants from a modified redlich-kwong equation of state[END_REF] usually serves as thermodynamic closure despite requiring specific treatment to avoid non-physical noise generation at interfaces [START_REF] Shyue | A fluid-mixture type algorithm for compressible multicomponent flow with van der waals equation of state[END_REF][START_REF] Pantano | An oscillation free shock-capturing method for compressible van der waals supercritical fluid flows[END_REF].

Unlike in the supercritical regime, in subcritical condition liquid and gas phase must be considered as a two-phase flow whose features, including moving material interface and vaporization effects, should be accounted for. Thus, our model will be built on the existing numerical methods for describing multi-phase flows. Moreover, the selected multiphase method is required to be capable of handling phase transition and combustion. Depending on how the material interface is treated, there exist two major categories to describe multi-phase flows briefly introduced hereafter. Sharp interface methods. Sharp interface methods assume the material interface as a sharp discontinuity. The interface is sharply represented by the front-tracking method [START_REF] Tryggvason | A front-tracking method for the computations of multiphase flow[END_REF], the level-set function [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations[END_REF], or the ghost fluid method [START_REF] Fedkiw | A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF]. These methods are able to preserve the sharply resolved interface with high accuracy and give impressive results for the simulation of jet atomization. However, there are several drawbacks which make them unsuitable for simulating reactive cryogenic jets. For example, significant difficulty arises when the interface becomes too wrinkled for the front-tracking method. For the level-set method and the ghost fluid method, they still require careful considerations to ensure robustness and conservation of mass. Moreover, it is hard for them to deal with dynamic creation and demise of interface due to phase transition [START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF].

Diffuse interface methods. The interface is assumed to be a diffuse layer between the liquid and gaseous phases. Phase field models and the second gradient theory [START_REF] Jamet | The second gradient method for the direct numerical simulation of liquid-vapor flows with phase change[END_REF], for instance, enter within that category, based on original works by Cahn and Hilliard [START_REF] Cahn | Free energy of a nonuniform system. i. interfacial free energy[END_REF]. Most studies within that framework consider incompressible fluids [START_REF] Jacqmin | Calculation of two-phase navier-stokes flows using phase-field modeling[END_REF], and treating high density ratios is not trivial [START_REF] Ding | Diffuse interface model for incompressible two-phase flows with large density ratios[END_REF]. Compressible flows may be considered within the second gradient theory, but with a high computational cost, since the capillary interface structure then has to be resolved [START_REF] Jamet | The second gradient method for the direct numerical simulation of liquid-vapor flows with phase change[END_REF]. The so-called compressible one-fluid model which uses a single set of equations to describe multi-phase flow also belongs to the diffuse interface method class. For example, the so-called five-equation model [START_REF] Allaire | A five-equation model for the simulation of interfaces between compressible fluids[END_REF][START_REF] Murrone | A five equation reduced model for compressible two phase flow problems[END_REF] assumes pressure equilibrium but allows different temperatures. The homogeneous model (four-equation model) further assumes the thermal and mechanical equilibrium. A particularly interesting feature of the diffuse interface methods is that the same equations are solved everywhere, with a unique hyperbolic solver. Also, as recently demonstrated in [START_REF] Chiapolino | A simple phase transition relaxation solver for liquid-vapor flows[END_REF][START_REF] Chiapolino | A simple and fast phase transition relaxation solver for compressible multicomponent two-phase flows[END_REF], the diffuse interface method can handle large distortion of interface as well as phase transition. However, since the material interface is treated as a diffuse zone, the thermodynamic state of mixture is required to be defined.

Thus in this study we will construct the numerical model for cryogenic reactive jets based on this diffusive interface method. The four-equation model will be extended to take into account combustion in the gas phase through inclusion of (i) molecular diffusion, (ii) detailed finite rate chemistry [START_REF] Boivin | An explicit reduced mechanism for H 2 -air combustion[END_REF], and (iii) a thermodynamic closure compatible with high temperatures [START_REF] Boivin | A thermodynamic closure for the simulation of multiphase reactive flows[END_REF].

The paper is organized as follows. First, we present the governing equations, the associated thermodynamic closure, followed by the phase transition, transport and kinetic models. Second, the numerical method used in the solution are detailed. Third, a number of 1D and 2D test cases are presented, to validate and illustrate the model capabilities. Conclusions and perspectives are drawn in the last section.

Strategy & assumptions

This section presents the governing equations and associated closure models. It follows a top-to-bottom approach: first the multi-component reactive set of conservation equations is given. Then, the thermodynamic closure is carefully explained, before details are given for the source terms for phase transition and combustion in the gaseous phase.

Governing equations

In order to simulate compressible reactive multi-component flow, a homogeneous mixture model which is also know as the four-equation model is employed and extended [START_REF] Kapila | Two-phase modeling of deflagration-to-detonation transition in granular materials: Reduced equations[END_REF][START_REF] Saurel | A general formulation for cavitating, boiling and evaporating flows[END_REF]. The homogeneous mixture model has several characteristics facilitating the simulation involving material interface and reaction [START_REF] Murrone | A five equation reduced model for compressible two phase flow problems[END_REF][START_REF] Kapila | Two-phase modeling of deflagration-to-detonation transition in granular materials: Reduced equations[END_REF][START_REF] Saurel | Modelling phase transition in metastable liquids: application to cavitating and flashing flows[END_REF][START_REF] Lund | A hierarchy of relaxation models for two-phase flow[END_REF][START_REF] Le Martelot | Towards the direct numerical simulation of nucleate boiling flows[END_REF].

Firstly, the homogeneous model assumes the mechanical and thermal equilibrium between different components, meaning the different components share the same pressure, temperature and velocity in one computational cell. Thus only one set of conservation laws for mass, momentum and energy is required, and extension to multiple components is straightforward. Discussions regarding the underlying assumptions may be found in [START_REF] Shyue | A fluid-mixture type algorithm for compressible multicomponent flow with van der waals equation of state[END_REF][START_REF] Kapila | Two-phase modeling of deflagration-to-detonation transition in granular materials: Reduced equations[END_REF][START_REF] Saurel | A general formulation for cavitating, boiling and evaporating flows[END_REF][START_REF] Shyue | A wave-propagation based volume tracking method for compressible multicomponent flow in two space dimensions[END_REF][START_REF] Saurel | A multiphase godunov method for compressible multifluid and multiphase flows[END_REF][START_REF] Abgrall | Computations of compressible multifluids[END_REF].

Secondly, the homogeneous mixture model, which belongs to the diffuse interface methods, treats the material interface as a diffuse zone. Thus the same numerical method can be applied regardless of the fluid thermodynamic state (liquid/gaseous). Most importantly, creation and destruction of material interface due to phase transition can be naturally handled [START_REF] Saurel | A multiphase model for compressible flows with interfaces, shocks, detonation waves and cavitation[END_REF].

In our work the homogeneous mixture model is extended and given as

                                 ∂ρ ∂t + div(ρu) = 0, ∂ρu ∂t + div ρu ⊗ u + pI = div(µτ), ∂ρE ∂t + div ((ρE + p)u) = div(µτ • u) -div(q), ∂ρY k ∂t + div(ρ(u + v k )Y k ) = ωc,k + ωp,k . (1)
Here, ρ, u and p are the density, the velocity vector and the pressure of the mixture respectively. The mass fraction for each component k is denoted as Y k . The source term for component k due to chemical reaction is denoted as ωc,k and due to phase transition is ωp,k . The total energy of the mixture E is defined as E = e + 1 2 u 2 with the mixture internal energy as e = N k=1 Y k e k . The viscous stress tensor τ is introduced to account for the viscosity effect, expressed as

τ = - 2 3 div(u)I + 2Π ( 2
)
where Π is the deformation rate tensor calculated as Π = 1 2 (grad(u) + (grad(u)) T ). To calculate the viscosity term, the mixture dynamic viscosity µ should also be defined. The diffusion effect for each component k is accounted for by introducing the diffusion velocity vector v k . The formulation to calculate µ and v k will be given in the following subsection. Finally, the energy flux q is calculated as

q = -λ • ∇(T ) + ρ N k=1 h k Y k v k ( 3
)
where T is the temperature of the mixture, λ is the heat conduction coefficient and h k is the enthalpy of component k.

The first term of Eq. ( 3) is the heat diffusion term expressed by Fourier's law while the second term is associated with the diffusion of species with different enthalpies.

1.2. Hybrid thermodynamic closure

Elemental volume description

Defining a thermodynamic closure in the above context comes down to expressing the pressure p = f (ρ, e, Y k ) as a function of the mixture internal energy e, the volume mass ρ, and composition Y k as required to fully close the Eulerian part of the system (1). In this paper, the elemental volume may be represented as in Fig. 1. In this approach, each component, with mass fraction Y k with respect to the elemental volume, occupies its own volume. That component, represented by Y k can be:

• A gaseous constituent -e.g. nitrogen. As pointed out in [START_REF] Chiapolino | A simple and fast phase transition relaxation solver for compressible multicomponent two-phase flows[END_REF][START_REF] Boivin | A thermodynamic closure for the simulation of multiphase reactive flows[END_REF], considering that each gas component occupies its own volume is not limiting as long as the gaseous components are assumed to follow the ideal or gas equation of state, the molar volume being a constant. It does, however, lead to significantly less error-prone implementation, as no fraction within phases is required (only fractions of the whole multi-phase mixture).

Y 1 Y 2 Y 3 Y 4 Y 5 Y 6 Y 7 Y 8
ρ, p, T • A liquid constituent. In this framework, we consider separate equation of state for the liquid and gas phases of the same constituent (e.g. liquid and gaseous N 2 ).

• A solid constituent. Although not encompassed in this work, it is possible to consider solid particles within this framework.

The main advantage of this approach is that one can choose which constituent(s) can be present in liquid phase.

For instance, for a burning fuel jet, it is quite intuitive that only the fuel can be encountered in liquid form: it is highly unlikely that any combustion products or minor species (typically present in the flame structure) become liquid since they typically appear at high temperatures. Associated CPU gains are important, as the approach reduces the number of conservation equations, and simplifies considerably the design of phase transition solvers.

In the present work, we consider that each component k follows the Noble-Abel Stiffened Gas (NASG), a recent improvement of the classical SG EOS [START_REF] Métayer | Elaborating equations of state of a liquid and its vapor for two-phase flow models[END_REF] which takes into account repulsive molecular effects. The NASG EOS allows in particular a better prediction for the liquid volume as well as the associated sound speed [START_REF] Boivin | A thermodynamic closure for the simulation of multiphase reactive flows[END_REF] compared to the SG EOS.

The NASG EOS

For a given pure component k, the general NASG EOS with non-constant heat capacity [START_REF] Boivin | A thermodynamic closure for the simulation of multiphase reactive flows[END_REF] reads

                                                 v k (p k , T k ) = (C p,k (T k )-C v,k (T k ))T k p k +p ∞,k + b k h k (p k , T k ) = ´Cp,k (T k )dT k + b k p k + q k e k (p k , T k ) = ´Cv,k (T k )dT k + (C p,k (T k )-C v,k (T k ))T k p ∞,k p k +p ∞,k + q k s k (p k , T k ) = ´Cp,k (T k )/T k dT k -(C p,k (T k ) -C v,k (T k )) ln (p k + p ∞,k ) + q k g k (p k , T k ) = h k (p k , T k ) -T k s k (p k , T k ) c 2 k (p k , T k ) = γ k (T k )(C p,k (T k ) -C v,k (T k ))T k 1 + b k (p k +p ∞,k ) (C p,k (T k )-C v,k (T k ))T k 2 (4)
where c k is the speed of sound of fluid k, v k is the specific volume, h k is the mass enthalpy, s k is the entropy, g k is the chemical potential and γ k , p ∞,k , C v,k , q k , q k and b k are constant coefficients determined by the thermodynamic properties of the fluid.

Under thermal and mechanical equilibrium, the elemental volume depicted in Fig. 1 follow the mixing rules

                               T = T k , ∀k p = p k , ∀k v = N k=1 Y k v k , e = N k=1 Y k e k . (5)
Solving the full system above for (p, T ) as function of (v, e), as required in solving the governing equation ( 1) is tedious, but becomes rather simple under the following assumptions:

1. Only one component may be found under liquid form, and will be given index k = 1. The corresponding vapor will have the index k = 2. The other gas components are assigned indices k ≥ 3.

2. The heat capacity is assumed to be a constant at low temperatures ( of the order of the considered liquid/vapor critical temperature).

3. The liquid is assumed to be locally absent at high temperatures, where the gaseous heat capacity is no longer constant.

For fuel jets, most constituents will only appear at high temperatures, and assumption ( 1) is reasonable. In some specific applications (for instance a subcritical cryogenic liquid oxygen-methane jet), two liquid phases may be present, and it will be shown in a future work how to deal with this problem.

As shown in [START_REF] Boivin | A thermodynamic closure for the simulation of multiphase reactive flows[END_REF], the assumption (2) leads to a satisfactory agreement on a wide range of pressure and temperature, as long as the conditions remain subcritical. At high temperatures, however, the heat capacity variation must be taken into account.

Assumption (3) is also reasonable, especially in the four equation context where thermal and mechanical equilibrium is imposed everywhere: it is unlikely that liquid mixes with burnt gases at the local level. Note however that pockets of liquid may be encountered in the burnt gases: the above considerations are at the elemental volume level only.

Under the above approximations, the EOS (4) for the liquid (k = 1) reads

                               v 1 (p, T ) = (C p,1 (T )-C v,1 (T ))T p 1 +p ∞,1 + b 1 h 1 (p, T ) = C p,1 T + b 1 p + q 1 e 1 (p, T ) = p+γ 1 p ∞,1 p+p ∞,1 C v,1 T + q 1 c 2 1 (p, T ) = γ 1 (C p,1 -C v,1 )T 1 + b 1 (p+p ∞,1 ) (C p,1 -C v,1 (T ))T 2 , (6)
whereas for all gaseous constituents k ≥ 2, it reads

                               v k (p, T ) = RT k W k p k , h k (p, T ) = ´Tk T 0 C p,k (T k )dT k + h f,k e k (p, T ) = h k (p k , T k ) -p k v k (p k , T k ) = ´Tk T 0 C p,k (T k )dT k -RT k W k + h f,k c 2 k (p, T ) = C p,k (T ) C v,k (T ) .R.T. ( 7
)
where R is the universal gas constant. T 0 is the reference temperature and h f,k is the formation enthalpies for species.

W k is molecular weight for component k. Heat capacity at constant pressure C p,k (T k ) is estimated from the classical NASA polynomials [START_REF] Esch | Thermodynamic properties in polynomial form for carbon, hydrogen, nitrogen, and oxygen systems from[END_REF].

Hybrid thermodynamic closure for the mixture

From the above equations, obtaining the conservative variables as a function of the primitive variables (p, T, Y k ) is straightforward using the mixing rules

                     v = N k=1 Y k v k , e = N k=1 Y k e k , h = N k=1 Y k h k . (8)
Computing (p, T ) from the conservative variables is however non-trivial, and requires special attention, as the computation needs to be done everywhere at every time step.

As a first attempt, an iterative procedure was implemented to obtain (p, T ) from (ρ, e, Y k ). It was however found simpler to follow the following procedure, depending on the value of liquid mass fraction Y 1 :

             (p, T ) = (p liq , T liq ), if Y 1 < Y c , (p, T ) = (p gas , T gas ), if Y 1 ≥ Y c , (9)
where Y c is a cutoff value -set to an arbitrarily small value 10 -7 hereafter -and (p liq , T liq ) are obtained explicitly, as in [START_REF] Chiapolino | A simple and fast phase transition relaxation solver for compressible multicomponent two-phase flows[END_REF], assuming constant heat capacities for all constituents k. In this approach, the temperature reads

T liq = e -N k=1 Y k q k N k=1 Y k C v,k p+γ k p ∞,k p+p ∞,k , (10)
and must be computed after the pressure, obtained analytically from (ρ, e, Y k ) as

p liq = b + √ b2 + 4ãc 2ã , (11)
with

                     ã = Cv , b = e - q v -b Cp -Cv -p ∞,1 Cv -p ∞,1 Y 1 C p,1 -C v,1 , c = e - q v -b p ∞,1 Cp -Cv -Y 1 C p,1 -C v,1 , (12)
and

Cv = N k=1 Y k C v,k , Cp = N k=1 Y k C p,k , q = N k=1 Y k q k , b = N k=1 Y k b k . (13)
Where Y 1 ≥ Y c , p gas is obtained as is classical for gaseous combustion:

p gas = ρe N k=1 Y k W k N k=1 Y k /W k γ k -1 , (14)
and T gas is obtained from a Newton solver, using the last time-step temperature as initial guess.

Thermodynamic coefficients

Thermodynamic coefficients were obtained following [START_REF] Boivin | A thermodynamic closure for the simulation of multiphase reactive flows[END_REF]. In the original NASG EOS presentation [START_REF] Métayer | The noble-abel stiffened-gas equation of state[END_REF], the thermodynamic coefficients are determined by fitting the experimental saturation curves. It was however shown later [START_REF] Boivin | A thermodynamic closure for the simulation of multiphase reactive flows[END_REF] that significant errors may appear when the flow departs significantly from the saturation properties. These errors will become unacceptable if accurate assessment of the energy necessary to heat a product over wide ranges of temperatures is required, as in the examples provided in the last section. That limitation was lifted in [START_REF] Boivin | A thermodynamic closure for the simulation of multiphase reactive flows[END_REF] by providing an alternative method to compute the thermodynamic coefficients, valid over a wide pressure and temperature range.

The saturation property p sat (T ) is obtained through a fitted Antoine equation

p sat (T ) = 10 A-B C+T , (15)
with the A, B, C parameters reported in the NIST database [START_REF] Linstrom | Nist chemistry webbook[END_REF]. This avoids the iterative procedure required by the initial NASG formulation [START_REF] Métayer | The noble-abel stiffened-gas equation of state[END_REF], and decreases the numerical cost whilst keeping an excellent accuracy [START_REF] Boivin | A thermodynamic closure for the simulation of multiphase reactive flows[END_REF].

Phase transition modelling

Phase transition takes place only between the liquid and its vapor. Thus the source terms contributed to liquid and its vapor are considered as

           ωp,1 = ρν(g 2 -g 1 ), ωp,2 = -ρν(g 2 -g 1 ), (16)
where g k denotes the phase k Gibbs free energy g k = h k -T s k with h k and s k respectively the specific enthalpy and entropy. ν is a function of the specific interfacial area, temperature and pressure, and represents a relaxation parameter that controls the rate at which thermodynamic equilibrium is reached.

Instead of solving Eq. ( 16), in the present work a fast phase transition relaxation model in [START_REF] Chiapolino | A simple phase transition relaxation solver for liquid-vapor flows[END_REF] is applied. This model assumes ν is very large and relaxation to thermodynamic equilibrium is immediate. Note that finite rate phase transition models are also available [START_REF] Furfaro | Modeling droplet phase change in the presence of a multi-component gas mixture[END_REF].

During phase transition, the mass fraction of liquid and vapor varies. The pressure and temperature also change to their equilibrium values. However, the mixture specific volume v, mixture energy e and mass fractions for other species Y k≥3 remain constant.

Y 1 → Y * 1 or Y 2 → Y * 2 . ( 18
)
Instead of directly computing the exact solution with iterative method such as [START_REF] Métayer | Dynamic relaxation processes in compressible multiphase flows. application to evaporation phenomena[END_REF],the current relaxation method provides a fair approximation for Y 2 through gradually reaching the exact solution. Firstly, the relaxation solver will check if there is a solution without liquid i.e. Y * 1 = Y min where Y min = 10 -8 is a very small value. If the partial pressure is below the saturation pressure in this case, no liquid is present and the solution is

Y * 1 = Y min and Y * 2 = 1 -Y min -N k=3 Y k .
Otherwise, the following system has to be solved

                     p partial = x * v .p * = p sat (T * ), v = Y * 1 v 1 (T * , p * ) + Y * 2 v 2 (T * , p * ) + N k=3 Y k v k (T * , p * ), e = Y * 1 e 1 (T * , p * ) + Y * 2 e 2 (T * , p * ) + N k=3 Y k e k (T * , p * ), (19)
where vapor molar fraction x * v is defined as

x * v = Y * 2 /W 2 Y * 2 /W 2 + N k=3 Y k /W k . (20)
The first equation in system Eq. 19 expresses the thermochemical equilibrium condition where the vapor partial pressure in the gas phase is equal to the saturation pressure at the current temperature. The other two equations in Eq. 19 state mass conservation and energy conservation laws during phase transition. Then the solution Y * 2 can be approximated by satisfying these constraint conditions respectively as

                     Y sat 2 (p, T ) = p sat (T )W 2 p-p sat (T ) k≥3 Y k /W k , Y m 2 (p, T ) = v-v g (p,T ) v 2 (p,T )-v 1 (p,T ) , Y e 2 (p, T ) =
e-e g (p,T ) e 2 (p,T )-e 1 (p,T ) ,

where

v g (p, T ) =        1 - N k=3 Y k        v 1 (p, T ) + N k=3 Y k v k (p, T ), and e g (p, T ) =        1 - N k=3 Y k        e 1 (p, T ) + N k=3 Y k e k (p, T ).
The exact solution of the problem is to find pressure and temperature conditions satisfying

Y m 2 (p * , T * ) = Y e 2 (p * , T * ) = Y sat
2 (p * , T * ) which is not trivial to solve. The current relaxation solver approximates the solution in this way

• Y m = Y m 2 (p, T )
and Y e = Y e 2 (p, T ) are evaluated for the initial values of (p, x v ) which is before phase transition relaxation process, and T = T sat (x v .p),

• Y sat = Y sat 2 (p, T ) is evaluated at the initial (p, T ). 9

Then we select the one which produces the smallest variation. The selection process is similar to Minmod limiter [START_REF] Van Leer | Towards the ultimate conservative difference scheme. v. a second-order sequel to godunov's method[END_REF]. In other words, we introduce:

             r 1 = (Y m -Y 2 ) (Y e -Y 2 ) , r 2 = (Y m -Y 2 ) Y sat -Y 2 , (22)
where Y 2 is the initial mass fraction before relaxation phase transition. The Minmod-like selection process states:

• If r 1 < 0, or r 2 < 0, no mass transfer happens:

Y * 2 = Y 2 .
• Otherwise, the one which produces the minimum variation of mass transfer is selected.

Under the latter condition, Y * 2 is calculated as:

Y * 2 = Y 2 + sgn Y m -Y 2 × Min |Y m -Y 2 |, |Y e -Y 2 |, |Y sat -Y 2 | . ( 23
)
After Y * 2 is determined, the other variables Y * 2 , p * , T * can be updated.

Simplified transport model

For this study, we implemented the simplest transport model available, adapted to ensure appropriate behavior in the limits of pure liquid and pure multi-component gas mixtures. Future work may include the coupling with more advanced transport libraries, such as Eglib [START_REF] Ern | Eglib: A general-purpose fortran library for multicomponent transport property evaluation[END_REF] or Cantera [START_REF] Goodwin | Cantera: An object-oriented software toolkit for chemical kinetics, thermodynamics, and transport processes[END_REF].

The mixture dynamic viscosity µ and thermal diffusion coefficient are calculated by the mixture rule with µ = α l µ l + α g µ g , and λ = α l λ l + α g λ g [START_REF] Murrone | A five equation reduced model for compressible two phase flow problems[END_REF] where µ l and λ l are the dynamic viscosity and the thermal diffusion coefficient of liquid phase respectively while µ g and λ g are of gas phase. The α l and α g are volume fraction of liquid phase and gas phase. In our case, α l = α 1 and α g = N k=2 α k . For liquid phase, the µ l and λ l can be assumed to be constant and be independent of temperature. Their value can be obtained from NIST website [START_REF] Linstrom | Nist chemistry webbook[END_REF]. For high speed gas flow, the µ g is temperature-dependent and is modelled with power-law as

µ g = µ 0 ( T T 0 ) β , (25)
where µ 0 and T 0 are reference values and β is constant. Then the thermal diffusion coefficient of gas phase λ g is defined as

λ g = µ g Pr N k=2 Y k C p,k (T ) N k=2 Y k , ( 26
)
where Pr is the Prandt number. The diffusion effect for each component k is accounted for by introducing the diffusion velocity vector v k . The expression diffusion velocity vector is

v k = -D k grad(X k ) X k + v c ( 27
)
where D k is an equivalent diffusion coefficient of species k into the rest of the mixture and X k is the mole fraction which is calculated as

X k = W W k Y k
where W is mean molecular weight. To guarantee the global mass conservation, the correction velocity vector v c is introduced [START_REF] Poinsot | Theoretical and numerical combustion[END_REF]. By enforcing the global mass conservation law, the v c is calculated as

v c = N k=1 D k W k W grad(X k ). ( 28
)
The diffusion coefficients for species k are approximated with constant Schmidt number assumption as

D k = µ ρSc k , (29)
a fair approximation to begin with [START_REF] Poinsot | Theoretical and numerical combustion[END_REF] in reacting flows. Note that, in the following, D 1 = 0 for the liquid constituent.

Combustion modelling

The chemical source term ωc,k is generally evaluated through considering a chemical system of N species involving

M reactions as k=N k=1 υ k, j Z k k=N k=1 υ k, j Z k for j = 1...M (30)
where Z k is a chemical symbol for species k, υ k, j and υ k, j are the molar stoichiometric coefficients. Then the resulting source terms are given by

ωc,k = W k M j=1 (υ k, j -υ k, j )Q j (31)
where Q j is the progress rate of reaction j and can be evaluated through

Q j = K f, j N k=1 [X K ] υ k, j -K r, j N k=1 [X K ] υ k, j . (32)
with K f, j and K r, j defined as the forward and reverse rates of reaction j. [X k ] is the molar concentrations of species k in the gas phase. The central part of combustion modelling is to calculate reaction rates of K f, j and K r, j which can be usually modeled with the empirical Arrhenius law.

In the present work, a finite rate reduced model for the H 2 -O 2 reaction is introduced to capture the flame structure in transient flow. The 12-step skeletal mechanism for H 2 -O 2 combustion [START_REF] Boivin | An explicit reduced mechanism for H 2 -air combustion[END_REF] derived from the detailed San Diego mechanism [START_REF] Williams | Chemical-kinetic mechanisms for combustion applications[END_REF] is employed in the current numerical solver. The mechanism, summarized in Tab.1, involves eight reacting species. The associated rate were obtained from the San Diego mechanism website [START_REF] Williams | Chemical-kinetic mechanisms for combustion applications[END_REF], and the transport properties listed in Tab. 2 from Cerfacs database [49].

Numerical methods

Finite volume method

The partial differential equation system Eq. 1 is solved with the finite volume method on unstructured grids. The integral form of Eq. 1 with divergence theorem reads where Ω is the unstructured cell element and Γ is the cell boundary. U represents the vector of conservative variables in Eq. , and F inv (U) and F visc (U) are inviscid and viscous numerical fluxes across cell boundaries Γ respectively. W

∂ ∂t ˆΩ UdΩ + ˆΓ F inv (U) dΓ = ˆΓ F visc (U) dΓ + ˆΩ WdΩ, (33) 1 H
is the source term. The viscous numerical fluxes are evaluated with the second order central difference scheme. For high speed compressible flow, the accuracy and stability of numerical solution is mainly influenced by evaluation of inviscid flux F inv (U). The inviscid flux should be calculated without introducing excessive numerical dissipation or numerical oscillations.

Low dissipation Riemann solver with homogeneous reconstruction strategy

The numerical dissipation errors of evaluation of inviscid flux come from the reconstruction scheme and the Riemann solver [START_REF] Deng | A finite volume multi-moment method with boundary variation diminishing principle for euler equation on threedimensional hybrid unstructured grids[END_REF][START_REF] Xie | A hybrid pressure-density-based mach uniform algorithm for 2d euler equations on unstructured grids by using multi-moment finite volume method[END_REF]. A reconstruction scheme higher than second order is not trivial on unstructured grids [START_REF] Xie | A multi-moment finite volume method for incompressible navier-stokes equations on unstructured grids: volume-average/point-value formulation[END_REF][START_REF] Xie | A multi-moment constrained finite volume method on arbitrary unstructured grids for incompressible flows[END_REF]. Thus we will introduce low dissipation Riemann solver to reduce the numerical dissipation errors.

Advection upstream splitting method (AUSM) proposed in [START_REF] Liou | A new flux splitting scheme[END_REF] is a robust and accurate method to treat both linear and nonlinear waves in complex flow [START_REF] Niu | Computations of two-fluid models based on a simple and robust hybrid primitive variable riemann solver with ausmd[END_REF][START_REF] Niu | Development of a less-dissipative hybrid ausmd scheme for multi-component flow simulations[END_REF]. The AUSM scheme splits the inviscid numerical flux into convective and acoustic part and adaptively reconstructs the flux according to the local Mach number. One sequel of AUSMtype Riemann is SLAU (Simple Low-dissipation AUSM) [START_REF] Shima | Parameter-free simple low-dissipation ausm-family scheme for all speeds[END_REF][START_REF] Kitamura | Towards shock-stable and accurate hypersonic heating computations: A new pressure flux for ausm-family schemes[END_REF] scheme which reduces the numerical dissipation in low Mach region and is free from reference parameters. As shown in [START_REF] Deng | Multimoment finite volume solver for euler equations on unstructured grids[END_REF], the SLAU is less dissipative than classic Riemann solver such as HLLC (Harten, Lax, and van Leer contact). Thus we extend SLAU scheme in our solver to reduce the numerical dissipation. According to the idea of SLAU, the inviscid numerical flux F inv (U) can be split into

F inv (U) = m + |m| 2 Ψ L + m -|m| 2 Ψ R + PN, (34)
with vectors defined as

Ψ = [1, u, H, Y 1 , ...Y k ] T , N = [1, n, 0, 0, ...0] T ,
where subscripts L/R denote the left and right states of physical fields at the cell-interface and H stands for total enthalpy H = (ρE + p)/ ρ. n is the normal vector to the cell boundary. The mass flux m is

m = 1 2 ρ L V L + V + + ρ R V R -V -- χ c (p R -p L ), (35)
where V = u • n stands for the velocity normal to cell boundary. Other quantities in Eq.( 35) are computed as follows,

V + = (1 -φ) V + φ |V L | , V -= (1 -φ) V + φ |V R | , V = ρ L |V L | + ρ R |V R | ρ L + ρ R , χ = 1 -M 2 , c = c L + c R 2 M = min        1.0, 1 c |u L | 2 + |u R | 2 2        , Ma L/R = V L/R c, φ = -max [min (Ma L , 0) , -1] • min [max (Ma R , 0) , -1] .
Similarly, the pressure flux P is computed by

P = p L + p R 2 + f + p -f - p 2 (p L -p R ) + |u L | 2 + |u R | 2 2 f + p + f - p -1 ρc, (36)
with

f ± p =          1 2 1 ± sign(Ma) , if |Ma| 1 1 4 (Ma ± 1) 2 (2 ∓ Ma), otherwise , ρ = ρ L + ρ R 2 .
It can be seen that numerical flux is adaptive to local Mach number Ma. The numerical dissipation is reduced as Ma becomes small.

The left and right states of physical variables (ρ L/R , u L/R , p L/R , Y k,L/R ) are reconstructed with the second order TVD-like scheme proposed in [START_REF] Deng | Limiter-free discontinuity-capturing scheme for compressible gas dynamics with reactive fronts[END_REF] . However, as shown in [START_REF] Deng | High resolution multi-moment finite volume method for supersonic combustion on unstructured grids[END_REF] this will lead to non-physical numerical oscillations across the material interface. These oscillations are caused by the non-linearity of TVD schemes [START_REF] Deng | High resolution multi-moment finite volume method for supersonic combustion on unstructured grids[END_REF][START_REF] Johnsen | Preventing numerical errors generated by interface-capturing schemes in compressible multi-material flows[END_REF]. To overcome this issue, we use a homogeneous reconstruction strategy which is compatible to the homogeneous fluid model.

Instead of reconstructing (ρ

L/R , u L/R , p L/R , Y k,L/R ), the (T L/R , u L/R , p L/R , Y k,L/R
) are reconstructed with TVD schemes.

Then the state variable ρ L/R is calculated with the EOS. It can be seen in following numerical tests that this simple strategy is consistency with homogeneous fluid model which assumes the thermal and mechanical equilibrium across the material interface.

In obtaining the characteristic speeds of the Riemann solver [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF], we use the classical Wood approximation [START_REF] Wood | A textbook of sound[END_REF],

expressed as

1 ρc 2 = N k=1 α k ρ k c 2 k , ( 37
)
where α k is the volume fraction of component k and is obtained as α k = ρY k /ρ k , the specific density ρ k being obtained with [START_REF] Oefelein | Mixing and combustion of cryogenic oxygen-hydrogen shear-coaxial jet flames at supercritical pressure[END_REF] at the local pressure and temperature. It is noteworthy that the behavior of mixture sound speed is nonmonotonic according to the Wood approximation. Formulations ensuring the monotonicity can be found in [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF].

Advection of isolated material interface

In order evaluate the ability of the numerical solver to maintain the equilibrium of velocity and pressure fields across the interface, a simple interface-only problem without phase transition is considered in this test. The numerical solutions of mixture density and pressure at t = 0.1 s (after one period) are presented in the Fig. 3.

Although the numerical model produces a diffusive zone for the material interface, it is able to maintain the mechanical equilibrium across the interface without any numerical oscillations in pressure field. Thus the current numerical solver is consistent with the homogeneous mixture model which assumes the mechanical and thermal equilibrium between different phases.

Multi-component Sod shock tube problem

The Sod shock tube problem where the exact solution is available is employed here to validate our multi-component flow solver. The exact solution can be obtained by solving exact Riemann problem with γ = 1.4 for single-species ideal gas [START_REF] Sod | A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws[END_REF]. Therefore we solve the multi-component flow with a gas mixture consisting of 21% O 2 and 79% N 2 , which results in a similar γ. Following the initial condition suggested in [START_REF] Wang | Partial characteristic decomposition for multi-species euler equations[END_REF], we design the following shock tube as

(T 0 , p 0 ) =         
(375 k, 101325 Pa) 0 ≤ x ≤ 0.5 (300 k, 10132.5 Pa) otherwise

. ( 39
)
The computation is conducted until t = 6 × 10 -4 s which is comparable to the dimensionless case at t = 0.2. The numerical solution of normalized pressure and density are presented in Fig. 4 in which the exact solution is also included. It can be seen that the numerical solution agrees well with the exact solution in terms of the position of shock and rarefaction wave.

Shock tube test with water/air mixture

In the following tests, several shock tube benchmark problems are employed to verify the ability of the developed numerical solver to solve flow structures containing shock waves, contact discontinuities and rarefaction fans with phase transition effect. We use these tests to show the correct implementation of fast relaxation phase transition model under the current unstructured solver with hybrid thermodynamic closure and low dissipative Riemann solver. The mixtures consisting of liquid water, water vapor and air are considered in the shock tube computation. The coefficients of NASG for these components are listed in Table 4.

In this test, a two-phase mixture far from the phase bounds with initial mass fractions Y 1 = 0. 4: Thermodynamic parameters for water and air. The calculation of coefficients is based on [START_REF] Boivin | A thermodynamic closure for the simulation of multiphase reactive flows[END_REF].

p L = 0.2 MPa and p R = 0.1 MPa is set at the left and right side of domain, which results in initial density and temperature discontinuities. A thermodynamic equilibrium state is considered as the initial condition. Thus the leftside temperature and right-side temperature are calculated as T L = T sat (x v p L ) and T R = T sat (x v p R ) respectively.

Then the mixture density can be computed with the defined EOS. The computation is conducted until t = 1 ms. The simulation results with and without phase transition are presented in Fig. 5. The numerical solution without phase transition is consisting of a right-moving shock wave, a contact discontinuity and a left-moving rarefaction fan, which agrees with the typical solution of the Sod shock tube problem [START_REF] Sod | A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws[END_REF]. With phase transition effect, evaporation is yielded with the shock compression; and condensation is caused due to the rarefaction expansion wave. These results agree well with the phenomena observed in [START_REF] Chiapolino | A simple phase transition relaxation solver for liquid-vapor flows[END_REF].

Shock tube test in an air dominated mixture

Unlike from the previous problem in which initial temperature is deduced from given pressure and mass fraction, in this test the initial temperature and pressure is given while the mass fraction is deduced from thermodynamic equilibrium. An initial pressure ratio of 2 and initial temperature of T = 293 K are set throughout the tube. The mass fraction of air is initially set to Y 3 = 0.98 in the whole tube. The results at time t = 1 ms are shown in Fig. 6. Again, the typical flow structures of the shock tube problem are reproduced by the current numerical solver. Evaporation is produced along with the shock wave while condensation happens along with the rarefaction fan.

Cavitation test with double expansion waves

In this test, we simulate a situation similar to the cavitation process. A cavitation bubble will be produced velocity is set to -1 m.s -1 at left and +1 m.s -1 at right. The numerical results at time t ≈ 3.5 ms are presented at Fig. 7.

It can be seen that a caviation bubble is created due to expansion waves. The current numerical solver reproduces the cavitation phenomena, which agrees with the similar simulations in [START_REF] Chiapolino | A simple phase transition relaxation solver for liquid-vapor flows[END_REF].

Free-propagating pre-mixed flame

Resolving free-propagating pre-mixed flame structure requires proper implementation of the thermodynamic closure for high temperature gas, diffusion fluxes for heat, momentum and multi-species, and combustion modelling for gaseous mixture. Thus the benchmark test of freely propagating flame is employed here to validate the current numerical model. Considering the initial thermodynamic equilibrium condition, we initialize the premixed flame as Table 5. It is noteworthy that in order to test the robustness of the current solver, a very small amount of liquid oxygen Y O 2 (liquid) = 1.0 × 10 -8 is added in the whole domain, which should not bring any obvious influence to the simulation results. The computation is conducted until the steady state is achieved. The simulation results are presented in Fig. 8 where the reference solution from Cantera [START_REF] Goodwin | Cantera: An object-oriented software toolkit for chemical kinetics, thermodynamics, and transport processes[END_REF] is included. It can be seen that the simulation results reproduce the correct flame structure and agree well with the reference solution. Thus the combustion phenomena can be correctly captured by the current numerical solver. x(m) T (K)

Mass fractions

Unsteady simulation of liquid oxygen and gaseous hydrogen rocket engine

With above benchmark tests which respectively show the current numerical solver is able to solve moving interface, phase transition and combustion phenomena, this model is now applied to simulate the liquid oxygen and gaseous hydrogen rocket engine under subcritical injection. The flow consists of a coaxial liquid oxygen jet surrounded by a high-speed hydrogen flow. Under subcrtical condition, the rocket engine is characterized by the liquid-gas interface, which is in contrast to transcritical or supercritical condition. Due to shear stresses in the two-phase flow between the two phases, the liquid phase will be destabilized into filaments and droplets. Then the combustion takes place in the gas phase after evaporation of liquid jet. Here, we will show the current numerical solver is able to reproduce this complex process which is challenging to numerical models. Firstly, the instantaneous simulation results before ignition are shown in Fig. 10. The results of mass fraction show the liquid jet is destabilized into small filaments due to shear stresses, which is the typical flow structures of evaporating jet and is reproduced by the current numerical solver. The vapor oxygen is only produced across the material interface along with the evaporating process. The vorticities caused by shear flow, which may be diffused by numerical diffusion errors, are resolved by implementation of the high resolution scheme and the low dissipative Riemann solver in the current numerical solver. Moreover, the pressure field shows that there is no obvious numerical oscillation across the material interface, which is important since as shown in [START_REF] Deng | High resolution multi-moment finite volume method for supersonic combustion on unstructured grids[END_REF] numerical oscillation across material interface will lead to non-physical destabilization.

Then the cryogenic jet is ignited to investigate the combustion process. The instantaneous field of gaseous oxygen and temperature field are presented in Fig. 11. The reaction takes place after the liquid oxygen evaporates and is mixed with gaseous hydrogen. Then the typical non-premixed flame forms in the gaseous oxygen and hydrogen mixing layer.

The non-premixed flame is attached around the cryogenic jet, which is regular near the expansion inlet and becomes irregular as jet is destabilized. The flame structure is further investigated in Fig. 12 where the distribution of mass fraction of different species and temperature is plotted. The plot is made along the vertical line in the half upper domain before the jet is destabilized. These results illustrate the physical process during evaporation and combustion.

A mixture zone consisting of liquid oxygen and vapor oxygen is observed across the liquid interface. The evaporated oxygen then reacts with hydrogen, producing water vapor and a diffusion flame with a maximum temperature around 3300 K. The distributions of radicals are consistent with the structure described in [START_REF] Gaillard | A diffuse interface lox/hydrogen transcritical flame model[END_REF], showing that the delicate flame structure is resolved with current solver. In a summary, the current numerical solver reproduces the processes taking place in the cryogenic jet, which includes destabilization of moving interfaces, evaporation and combustion.

Concluding Remarks

In this work, we make contributions to construction of a numerical model for the simulation of multiphase reactive flows in the sub-critical cryogenic. This model is built on the diffusive interface method and is an extension of the four-equation model presented in [START_REF] Saurel | A general formulation for cavitating, boiling and evaporating flows[END_REF]. The complex multi-physics processes have been accounted for by integrating phase transition [START_REF] Chiapolino | A simple and fast phase transition relaxation solver for compressible multicomponent two-phase flows[END_REF] model, a hybrid thermodynamic closure strategy, a simplified transport model and a detailed chemistry model for H 2 -O 2 combustion. Then the model is solved by non-oscillatory finite volume method where a low dissipation Riemann solver is extended for multi-component flow.

The obtained model reproduces the expected results in both the multiphase and reactive limits. This has been extensively tested on shock tubes including phase transition, and on the propagation of planar premixed flames, proving that the hybrid thermodynamic closure implementation presented here is in fact producing the expected results.

Lastly, the stability of the model is shown on an unsteady simulation of a cryogenic liquid jet in conditions far below critical conditions known to be highly challenging numerically [START_REF] Gaillard | A diffuse interface lox/hydrogen transcritical flame model[END_REF].

Given the limitations of the current model, several points will be addressed in our future work to construct a more general model for reactive multi-phase flow. Firstly, the capillary effects or atomization models will be included when further considering distinct material interface characteristics in sub-critical condition. Secondly, the NASG EOS and phase transition model will be extended to account for multi-liquid mixture. Thirdly, temperature non-equilibrium model will be considered to further reduce the pressure oscillations across cold liquid and hot gas. x(m) T (K)

Mass fractions

Figure 1 :

 1 Figure 1: Thermodynamic closure: each component occupies its own volume, regardless of its state. (p, ρ, T ) are mixture averaged properties.

  A twodimensional tube [0, 1] × [0, 0.5] consisting of uniform triangular grids is created as the simulation domain. The mesh size is ∆h = 1 100 in our calculation. The mesh for the shock tube problem is illustrated in Fig. 2. The problem consists of a square liquid column in vapor advected with a uniform velocity u = u 0 = 10 m/s under equilibrium pressure p = p 0 = 10 6 Pa and uniform temperature T = T 0 = 80 K. For initial condition, liquid is set in the region of x ∈ [0.4, 0.6] m and the vapor is filled elsewhere. The mass fraction of liquid is set as Y 1 = 0.9999998 for the liquid region and Y 2 = 0.9999998 in the vapor region. A very small mass fraction for non-condensable gas Y 3 = 10 -7 is added to the whole domain. Periodic boundary conditions are used on the left and right boundaries during the computations.

Figure 3 :

 3 Figure 3: Numerical solutions for the advection of isolated material interface problem. The mixture density field is presented in (a) and pressure field in (b).

Figure 5 :Figure 6 :

 56 Figure 5: Numerical solutions for the shock tube test with a mixture far from the phase bounds. The solution without phase transition is presented with the solid line. The solution with phase transition is in the dashed line.

Figure 7 :

 7 Figure 7: Numerical solutions for cavitation test with double expansion waves. The solution without phase transition is presented with the solid line. The solution with phase transition is in the dashed line.

Figure 8 :

 8 Figure 8: Freely propagating flame: Temperature profile (thick line), H 2 ( ), O 2 ( ), H 2 O ( * ), H ( ), HO 2 (•) mass fractions. Fresh gases are in stoichiometric proportion, at 300K and atmospheric pressure. Cantera reference (plain line), and the current numerical solver (red dashed line).

A

  two-dimensional computation is conducted to simulate liquid oxygen and gaseous hydrogen jet. The half of geometry of computation domain is shown in the Fig. 9. A central cold flow made of nearly pure liquid oxygen, at 100 K and 30 m.s -1 , pressure 3 MPa is imposed as the inlet subsonic boundary along segment AB of Fig. 9, while along segment DE of Fig. 9 a peripheral flow made of nearly pure gaseous hydrogen, at 150 K and 200 m.s -1 , with the pressure of 3 MPa is imposed. Along segments GH and HI, a non-reflecting subsonic boundary at 3 MPa is considered. The remaining boundaries are treated as symmetric. The mesh is composed of around 3.6 × 10 5 triangular elements. The NASG EOS is used for liquid oxygen and low temperature oxygen vapor of which parameters are listed in

Figure 10 :

 10 Figure 10: Numerical results for evaporating liquid oxygen jet before ignition. Results for mass fraction of liquid oxygen and vapor oxygen, vorticity field and pressure field are presented respectively.

Figure 11 :

 11 Figure 11: Numerical results for liquid oxygen and gaseous hydrogen jet after ignition. The mass fraction of gaseous oxygen and temperature field are presented.

Figure 12 :

 12 Figure 12: Radial profiles in the jet at x = 10mm: Temperature (thick, right axis) LO 2 (solid); GO 2 (dashed), H 2 O (dotted), H 2 (dot-dashed) mass fractions in the top subplot; , OH (solid), O(dashed) and H(dot-dashed) mass fractions in the bottom subplot.

  Thus the fast phase transition relaxation uses a fractional step method to get the equi-and T * are functions of (v, e, Y * k ), and v and e are constant during the phase transition. Also the Y * 1 and Y * 2 are linked with the formulation Y * 2 = 1 -Y * 1 -k≥3 Y k . Thus the equilibrium state (Y * 1 , Y * 2 , p * , T * ) can be fully determined if either Y * 1 or Y * 2 is known. Then the main relaxation process of Eq. (17) can now be

	librium state (Y * 1 , Y * 2 , p

* , T * ) from the state (Y 1 , Y 2 , p, T ). The fast phase transition relaxation process reads

(Y 1 , Y 2 , p, T ) → (Y * 1 , Y * 2 , p * , T * )

(17)

Since p *

Table 1 :

 1 The 12-step skeletal mechanism for the combustion of H 2 -O 2[START_REF] Boivin | An explicit reduced mechanism for H 2 -air combustion[END_REF]. Up-to-date rates are available[START_REF] Williams | Chemical-kinetic mechanisms for combustion applications[END_REF].

	+ O 2	OH + O	7	HO 2 + OH → H 2 O + O 2
	2 H 2 + O	OH + H	8	H + OH + M	H 2 O + M
	3 H 2 + OH		H 2 O + H	9	2 H + M	H 2 + M
	4 H + O 2 + M → HO 2 + M 10 2 HO 2 → H 2 O 2 + O 2
	5 HO 2 + H → 2 OH	11 HO 2 + H 2 → H 2 O 2 + H
	6 HO 2 + H		H 2 + O 2	12 H 2 O 2 + M → 2 OH + M
	µ l		1.9650 × 10 -4	λ l	0.1518
	µ 0	1.8405 × 10 -5	β	0.6759
	Pr	0.7500		Sc H 2	0.2100
	Sc H	0.1400		Sc O 2	0.8000
	Sc OH	0.5300		Sc O	0.5300
	Sc H 2 O 0.6000		Sc HO 2 0.8000
	Sc H 2 O 2 0.8200		Sc N 2	1.0000

Table 2 :

 2

(λ l , µ l ) for liquid oxygen -its diffusion parameter is set to zero. Gas phase related coefficients: power-law viscosity coefficients (SI units), Prandtl number and Schmidt numbers for each species of the 12-step mechanism.

Table 3 :

 3 Numerical errors and convergence rates of the density perturbation transport test.

	grid size	L 1 error	L 1 order
	1/20	3.66 × 10 2 -
	1/40	1.02 × 10 2 1.84
	1/80	3.12 × 10 3 1.71
	1/160	9.11 × 10 4 1.78

Figure 2: Computational grids for shock tube problems.

Table

  Numerical solutions for multi-component shock tube problem with 21% O 2 and 79% N 2 . The black symbols represent the normalized quantity of pressure and density. The black solid lines represents the exact solution.

						numerical					numerical
		1				exact		1				exact
		0.8						0.8
	pressure	0.6					density	0.6
		0.4						0.4
		0.2						0.2
		0	0.2	0.4	0.6	0.8	1	0	0.2	0.4	0.6	0.8	1
				x						x
				(a) normalized pressure				(b) normalized density
	Figure 4: Coefficients	H 2 O (l)		H 2 O (g)	air
				C p (J/kg/K)	4185.3		1908.3	1007
				γ			1.0123		1.3281	1.4
				P ∞ (Pa)	1835 × 10 5		0	0
				q (J/kg)		-1143030	1957400	0
				b (m 3 /kg)	9.2003 × 10 -4	0	0
				W (g/mol)		18		18	29
				A			4.6543		-
				B			1435.264		-
				C			-64.848		-

1 (liquid water), Y 2 = 0.2 (vapor water) and Y 3 = 0.7 (air) is considered throughout the shock tube. An initial pressure jump of

Table 5 :

 5 Initial conditions: 1-D domain is initialized with fresh gases corresponding to (0 : L/2) and burnt gases (L/2 : L)

	Variables fresh gases	burnt gases
	T	300 K	2385 K
	p	1 atm	1 atm
	Y H 2 Y H	2.852 × 10 -2 1.145 × 10 -3 0 6.983 × 10 -5
	Y O 2 (gas) Y OH	2.264 × 10 -1 7.474 × 10 -3 0 5.458 × 10 -3
	Y O	0	3.838 × 10 -4
	Y H 2 O Y HO 2 Y H 2 O 2	0 0 0	2.403 × 10 -1 1.074 × 10 -6 1.444 × 10 -10
	Y N 2	7.451 × 10

-1 7.451 × 10 -1

Table 6 ,

 6 while NASA polynomials are used to calculate C p for other gaseous components.

	I

Table 6 :

 6 Thermodynamic parameters for liquid oxygen and oxygen vapor. The calculation of coefficients is based on[START_REF] Boivin | A thermodynamic closure for the simulation of multiphase reactive flows[END_REF]

Acknowledgment

The presence of Xi Deng at M2P2 was made possible by Labex MEC (ANR-10-LABX-0092) and the A*MIDEX project (ANR-11-IDEX-0001-02), funded by the Investissements d'Avenir. This work was granted access to the HPC resources of Aix-Marseille Université funded by the project Equip@Meso (ANR-10-EQPX-29-01).

Implicit large eddy simulation with dissipated non-oscillatory finite volume schemes

Instead of using explicit large eddy simulation (LES) which requires additional sub-grid scale model from filtering the equation system, we are employing a method called implicit large eddy simulation (ILES) [START_REF] Grinstein | Implicit large eddy simulation: computing turbulent fluid dynamics[END_REF] method here. We solve the equation system with dissipated non-oscillatory finite volume (NFV) numerical schemes. The embedded numerical dissipation in NFV schemes is used in replace of the explicit SGS models. As shown in [START_REF] Ritos | Implicit large eddy simulation of acoustic loading in supersonic turbulent boundary layers[END_REF], the ILES apporach is able to produce results close to those from explicit LES.

Time integration and solution procedure

In each time step ∆t, we first advance the solution with the Strang splitting scheme employed in [START_REF] Fedkiw | High accuracy numerical methods for thermally perfect gas flows with chemistry[END_REF]. Then a phase transition relaxation model is adopted to update the solution. We summarize the solution procedure in the ∆t as following: (i) update the solution with the reaction source term in a half time step ∆t/2; (ii) use the obtained solution as initial data, advance the solution with the advection and diffusion flux in a time step ∆t with a third-order Total Variation Diminishing (TVD) low-storage Runge-Kutta scheme [START_REF] Gottlieb | Total variation diminishing runge-kutta schemes[END_REF]; (iii) update the solution again with the reaction source term over a half time step ∆t/2; (iv) apply phase transition relaxation model to recalculate the variables Y k ,p and T in thermodynamic equilibrium.

Numerical experiments

Accuracy test for solving multi-component Euler equation

In order to test the convergence rate of the developed model, we solve the propagation of mass fraction disturbances of gas oxygen and hydrogen with the multi-component Euler equation. The initial temperature and pressure are specified uniformly through the whole computational domain, and a sinusoidal perturbation is given to the mass fraction, as follows

The computational domain [-1, 1] × [-1, 1] is divided into uniform triangular elements. The convergence studies are conducted by gradually refining grids with periodic boundaries. The evolution time is t = 0.2 s. The numerical errors and convergence rates with respect to the partial density of ρY O 2 for the multi-component Euler equation are summarized in Table 3. We observe that the convergence rate with the constructed model is nearly second order.