Alexandre Chiapolino

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

Richard Saurel

A simple and fast phase transition relaxation solver for compressible multicomponent two-phase flows

Keywords: multiphase flows, interfaces, hyperbolic systems, relaxation, phase transition

, such a solver was successfully developed to compute thermodynamic equilibrium between a liquid phase and its corresponding vapor phase. The present work extends the solver's range of application by considering a multicomponent gas phase instead of pure vapor, a necessary improvement in most practical applications. The solver proves easy to implement compared to common iterative procedures, and allows systematic CPU savings over 50%, at no cost in terms of accuracy. It is validated against solutions based on an accurate but expensive iterative solver.

Its capability to deal with cavitating, evaporating and condensing two-phase flows is highlighted on severe test problems both 1D and 2D.

Introduction

Phase transition between liquid and gas phases is a well-established phenomenon, which tends to relax thermochemical potentials toward equilibrium. Yet, finding the thermodynamic equilibrium can be tedious, as the corresponding equations are non-linear, the culprit being the non-trivial relation linking pressure (p) and temperature (T ) at saturation. In multiphase compressible flows modeling, phase transition has dominant effects in many applications, and has to be accounted for via appropriate phase transition models. Iterative equilibrium solvers are usually computationally expensive, as they require complex root-finding procedures, sometimes draining more CPU time than the flow computation itself. Their coding is not trivial, as for example, pure phase bounds require special attention.

Our previous work [START_REF] Chiapolino | A simple phase transition relaxation solver for liquidvapor flows[END_REF] presented a robust, fast and accurate alternative to the exact equilibrium solver in the case of flows consisting of a liquid phase in equilibrium with its own vapor.

Speedup on the order of 50% was achieved on computations of such flows. When only liquid and the corresponding vapor are present, p and T are linked by the saturation curve p = p sat (T ), simplifying somewhat the system to solve.

In the present paper, we aim at extending the model to liquid in equilibrium with multicomponent gases. In this case, p and T are not directly on the saturation curve, but are related through the partial pressure of the vapor component in the multicomponent gas phase. This allows the well-known existence of water vapor in air at atmospheric conditions, albeit a temperature below the boiling point.

The range of flow solvers to which the present model applies is identical to its previous version: it is designed in association with non-equilibrium hyperbolic flow models, such as Baer and Nunziato's (1986) [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] and its reduced versions. This includes the 5-equation model of [START_REF] Kapila | Two-phase modeling of deflagrationto-detonation transition in granular materials: Reduced equations[END_REF] [START_REF] Kapila | Two-phase modeling of deflagrationto-detonation transition in granular materials: Reduced equations[END_REF] and its extension for cavitating flows [START_REF] Saurel | Modelling phase transition in metastable liquids: application to cavitating and flashing flows[END_REF][START_REF] Lemartelot | Liquid and liquid-gas flows at all speeds[END_REF][START_REF] Saurel | A general formulation for cavitating, boiling and evaporating flows[END_REF], as well as formulations for thermal equilibrium two-phase mixtures such as the Homogeneous Relaxation Model (HRM) and Homogeneous Equilibrium Model (HEM) [START_REF] Downar-Zapolski | The non-equilibrium relaxation model for one-dimensional flashing liquid flow[END_REF][START_REF] Barret | Schemes to compute unsteady flashing flows[END_REF]. This paper is organized as follows. The first two sections remind the 4-equation flow model [START_REF] Saurel | A general formulation for cavitating, boiling and evaporating flows[END_REF] and the corresponding thermodynamic closure, backbone of the present study. The 4-equation two-phase model is indeed the simplest flow model for the relaxation solver. Also, a wide range of applications can be addressed with this approach [START_REF] Saurel | A general formulation for cavitating, boiling and evaporating flows[END_REF]. The third section presents the system of equations to be solved to relax thermochemical potentials, while the fourth details how to do so. The last sections are dedicated to displaying and validating the algorithm capabilities, through a series of 1D and 2D test cases involving cavitating, evaporating and condensing flows.

Flow model

Temperature-pressure equilibrium two-phase flow model

The phase transition relaxation solver may be used with models mentioned previously (with 7, 5, 4 and 3 partial differential equations) but its presentation is simplified in the context of the 4-equation model (often called HRM). When dealing with more sophisticated formulations, as for example the 7-equation model, extra ingredients have to be presented, such as velocity, pressure and temperature relaxation solvers, as done for example in [START_REF] Métayer | Dynamic relaxation processes in compressible multiphase flows[END_REF][START_REF] Lallemand | Pressure relaxation procedures for multiphase compressible flows[END_REF]. Here, there is a single step that makes the connection from the 4-equation model (modeling mixtures out of thermodynamic equilibrium) and the 3-equation model (mixtures in full equilibrium). The 4-equation model for a two-phase mixture reads:

                     ∂ρ ∂t + div(ρu) = 0, ∂ρu ∂t + div ρu ⊗ u + pI = 0, ∂ρE ∂t + div ((ρE + p)u) = 0, ∂ρ l α l ∂t + div(ρ l α l u) = 0, (1.1)
where the liquid mass conservation equation can be written alternatively as,

∂ρY l ∂t + div(ρY l u) = 0. (1.2)
Here Y l,g , α l,g , ρ l,g denote respectively the mass fraction, the volume fraction and the material density of the liquid (l subscript) and gas (g subscript) phases. ρ represents the mixture density, u represents the mixture center of mass velocity, p denotes the mixture pressure and E the mixture total energy (E = e + u 2 /2). The mixture specific internal energy is defined as e = Y l e l + Y g e g .

Mass transfer has been omitted in System (1.1) as it is addressed later.

Extension to a multicomponent gas phase

Considering multicomponent effects within the gas phase needs additional mass balance equations that can be written as, ∂ρ g α g y k ∂t + div(ρ g α g y k u) = 0, (1.3) where y k denotes the k-th gas component mass fraction within the gas phase. In the thermodynamic closure considered in this study, and detailed in the next section, the gas phase is assumed to obey the ideal gas equation of state.

As the molar volume is independent on the gas constituent, it is equivalent to consider that the constituents are ideally mixed within the gas phase (each occupying the whole gas phase volume), or assume that each component within the gas phase occupies its own separate volume.

The equivalence between these two approaches is clarified in the next section. Under the latter assumption, all mass balance equations can be written as:

∂ρY k ∂t + div(ρY k u) = 0,
where the subscript k refers to the various physical and chemical components. In the following, let us introduce the following convention:

-k = 1 for the liquid, -k = 2 for the gas component corresponding to vapor of species 1, -k = 3, . . . , N for the remaining gas components, considered non-condensable in this work.

This notation for the species conservation equation is more convenient, as Y k is now the mass fraction for the k-th species in the entire mixture (containing both liquid and gas), so that the conservation equation for each species within the gas phase is of the same form as the liquid mass conservation equation.

The extended two-phase flow model consequently reads,

                     ∂ρ ∂t + div(ρu) = 0, ∂ρu ∂t + div ρu ⊗ u + pI = 0, ∂ρE ∂t + div ((ρE + p)u) = 0, ∂ρY k ∂t + div(ρY k u) = 0, (1.4)
with

E = e + 1 2 u 2 , e = N k=1
Y k e k .

System (1.4) is clearly reminiscent of the reactive (or multicomponent) Euler equations widely used in chemically reacting flows. However, the thermodynamic closure differs significantly from the one used in gas mixtures as examined hereafter.

Thermodynamic closure

This section presents the equations of state (EOS) used for each pure component as well as the mixture equation of state. Preliminarily, let us demonstrate the equivalence between the ideal gas mixture model that follows Dalton's law, and a gas phase within each constituent is assumed to occupy its own volume, in the specific context of temperature and pressure equilibrium.

The gas phase is ideally mixed: Dalton's law

Associated mixture rules are: The ideal gas law for a given species reads,

                   T = T k : mixture in temperature equilibrium V = V k : each chemical
p partial, k V = n k RT,
and the gas mixture pressure reads,

pV = RT N k=2 n k ,
where R denotes the universal gas constant in molar units and n k is the number of moles. In this relation the temperature is deduced from the caloric equation of state for the mixture,

e = RT N k=2 Y k /W k γ k -1 ,
where W k denotes the molar mass and γ k the polytropic coefficient of chemical species k. Consequently, the pressure reads,

p = e V N k=2 n k N k=2 Y k /W k γ k -1 = e N k=2 n k W k V 1 W k N k=2 Y k /W k γ k -1 , finally, p = ρe N k=2 Y k W k N k=2 Y k /W k γ k -1 .
(2.1)

Each gas constituent occupies its own volume in pressure and temperature equilibrium

We now consider another configuration with separate chemical species. Corresponding mixture rules are, The ideal gas law for a given species now reads,

                   T = T k : mixture in temperature equilibrium v = N k=2 Y k v k :
p k V k = n k RT,
In mass units it becomes,

p k = n k W k V k R W k T = ρ k R W k T.
Consequently the specific volume of a given chemical species as a function of pressure and temperature reads,

v k = R W k T p k .
This result is inserted into the specific volume definition,

v = N k=2 Y k v k , yielding, p T = ρ N k=2 Y k R W k .
The mixture energy definition is now considered, resulting in the same relation as in Section 2.1:

e = RT N k=2 Y k /W k γ k -1 .
Eliminating the temperature by combining these two relations, the following result is obtained:

p = ρe N k=2 Y k W k N k=2 Y k /W k γ k -1 . (2.2)
It is clear that both ideal mixtures models (Dalton's law Eq. (2.1) and separate chemical species Eq. (2.2)) are equivalent when the fluids are ideal gases evolving in both temperature and pressure equilibrium. This remark simplifies dramatically the building of the mixture EOS.

Before addressing this topic, let us present the EOS for the liquid phase and its vapor.

EOS for pure constituent

In this frame, it is assumed that each gaseous constituent (k = 2, . . . N) obeys the ideal gas equation of state. The EOS coefficients for the vapor (species k = 2) must be carefully computed, in accordance with the liquid (k = 1) EOS, as to fit the phase diagram. The building of such an EOS has been addressed in Le Métayer et al. ( 2004) [START_REF] Métayer | Elaborating equations of state of a liquid and its vapor for two-phase flow models; Elaboration des lois d'état d'un liquide et de sa vapeur pour les modèles d'écoulements diphasiques[END_REF], on the basis of the Stiffened Gas (SG) EOS. In later developments, the same authors proposed the "Noble-Abel Stiffened

Gas" (NASG) [START_REF] Métayer | The Noble-Abel Stiffened-Gas equation of state[END_REF], which improves considerably the liquid specific volume accuracy by taking into account the repulsif molecular effects in addition to those already present in the SG EOS (agitation and attraction). The main formulas for the NASG EOS read for a given constituent k = 1, . . . N (liquid, vapor, gas),

                           p k (v k , e k ) = (γ k -1) (e k -q k ) v k -b k -γ k p ∞,k , T k (p k , v k ) = (v k -b k ) (p k + p ∞,k ) C v,k (γ k -1) , g k (p k , T k ) = (γ k C v,k -q ′ k )T k -C v,k T k ln T γ k k (p k + p ∞,k ) γ k -1 + b k p k + q k , c k (p k , v k ) = v k γ k v k (p k + p ∞,k ) v k -b k . (2.3)
where the following parameters are needed for each phase: γ k , p ∞,k , C v,k , q k , q ′ k and b k . These parameters are constant coefficients characteristic of the thermodynamic properties of the fluid.

Among them the coefficient b represents the covolume of the fluid. g k denotes the phase Gibbs free energy of fluid k, g k = h k -T s k with h k and s k respectively the specific enthalpy and entropy. c k denotes the speed of sound of fluid k.

The following results can be adapted to the SG EOS framework by setting b k = 0. For all gaseous constituents (k = 2, . . . N), b k = 0 and p ∞,k = 0, and the above system reduces to the ideal gas EOS. Note that, for the k-th gaseous constituent, the pressure p k is based on the subvolume containing the constituent. It is therefore equal to the pressure of the whole gas phase, and not to the partial pressure of the constituent in the gas phase. In the following, it is explicitly stated when a partial pressure is introduced.

The same strategy as in [START_REF] Métayer | Elaborating equations of state of a liquid and its vapor for two-phase flow models; Elaboration des lois d'état d'un liquide et de sa vapeur pour les modèles d'écoulements diphasiques[END_REF] and [START_REF] Métayer | The Noble-Abel Stiffened-Gas equation of state[END_REF] is repeated in order to find a formulation connecting the saturation pressure and temperature, by equating the liquid and the vapor chemical potentials g 1 = g 2 , leading to the following equation:

ln (p sat + p ∞,2 ) = A + B + E p sat T sat + C ln (T sat ) + D ln (p sat + p ∞,1 ) , (2.4)
where,

A = C p,1 -C p,2 + q ′ 2 -q ′ 1 C p,2 -C v,2 , B = q 1 -q 2 C p,2 -C v,2 , C = C p,2 -C p,1 C p,2 -C v,2 , D = C p,1 -C v,1 C p,2 -C v,2 , E = b 1 -b 2 C p,2 -C v,2 .
(2.5)

Mixture equation of state

Each constituent is assumed to follow the NASG EOS (reduced to the ideal gas EOS for gaseous constituents). However the mixture equation of state, based on mechanical and thermal

equilibrium has yet to be built. Under the assumption of mechanical and thermal equilibrium, an inherent assumption to the 4-equation model, we have:

                         T = T k ∀k, p = p k ∀k, v = N k=1 (Y k v k ), e = N k=1 (Y k e k ).
(2.6) with v the mixture specific volume.

T (K) T (K) p sat (bar) L v (kJ/kg) h 1 (kJ/kg) h 2 (kJ/kg) v 1 (m 3 /kg) v 2 (m 3 /kg)
As shown at the beginning of this section, considering (2.6) for the mixture of gases is equivalent to Dalton's law. Therefore the computed mixture pressure for the gas mixture is in agreement with both Dalton's law in the gas mixture and with the liquid-gas interface condition of equal pressures. Consequently, System (2.6) summarizes correctly the separate phase mixture model between the liquid and the gas mixture on one hand and the ideal gas mixture on the other hand.

From the expressions given in Eqs. (2.3) for pure constituents, the specific volumes and internal energies read:

       v k (p k , T k ) = (γ k -1)C v,k T k p k + p ∞,k + b k e k (p k , T k ) = p k + γ k p ∞,k p k + p ∞,k C v,k T k + q k , (2.7)
which can be rewritten to give two expressions for the temperature (since T k = T, ∀k):

T = v -N k=1 Y k b k N k=1 Y k (γ k -1)C v,k p+p ∞,k , (2.8)
T = e -N k=1 Y k q k N k=1 Y k C v,k p+γ k p ∞,k p+p ∞,k .
(2.9)

Equating these two expressions, and taking into account that p ∞,k = 0 and b k = 0, ∀k > 1, a quadratic expression for the mixture pressure is obtained as:

p = b + √ b 2 + 4ac 2a , (2.10)
with:

               a = 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 ) , (2.11)
where mixture quantities are introduced:

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 . (2.12)
Although not trivial, Eq. (2.10) ensures the strict positivity of the mixture pressure p: if Y 1 = 1, both a and c are strictly positive leading to p > (b + |b|)/a. The mixture pressure is then strictly positive. When Y 1 = 1, a > 0 and c = 0, and the strict positivity of the pressure is less trivial to demonstrate. However, combining Eqs. (2.7) and (2.11) in this limit leads to b > 0, and consequently p > 0, where p follows the NASG EOS for the pure liquid phase. Let us add that when Y 1 = 0 (absence of liquid), the pressure equation (2.10) reduces to

p = e - q v ( Cp -Cv ) Cv ,
which after substituting e = N k=2 Y k e k , and noticing that ( Cp -Cv ) = N k=2 Y k R/W k for ideal gases, leads to:

p = N k=2 Y k ρ RT /W k ,
or, in other words, the classical law of Dalton for the gas mixture. This provides an additional verification as to the possibility of considering each gas constituent to be in its own volume, as discussed in the flow model section. Albeit the apparent simplicity of the thermodynamic closure chosen, the phase transition model presented here may be extended to other thermodynamic closures provided each phase EOS is convex.

Phase transition model

Objective

When phase transition is addressed between the liquid and its vapor, only the conservation equations of the first two constituents are modified as:

     ∂ρY 1 ∂t + div(ρY 1 u) = ρν(g 2 -g 1 ) ∂ρY 2 ∂t + div(ρY 2 u) = -ρν(g 2 -g 1 ) (3.1)
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 and ν(A I , T, p) represents a relaxation parameter that controls the rate at which thermodynamic equilibrium is reached. It is a function of the specific interfacial area A I , temperature and pressure. Its determination is possible only when the interfacial area A I is available, as with droplets and bubbly flows (see [START_REF] Furfaro | Modeling droplet phase change in the presence of a multi-component gas mixture[END_REF] for example [START_REF] Furfaro | Modeling droplet phase change in the presence of a multi-component gas mixture[END_REF]) and sometimes for stratified flows. A natural way to determine this exchange area would be to use a very high resolution mesh to capture interfaces at all spatial scales. However this strategy would require tremendous computing resources and can hardly be envisaged when considering large-sized problems. Nevertheless a realistic method in specific limit situations is to consider instantaneous thermodynamical relaxation between phases by the use of additional source terms [START_REF] Métayer | Dynamic relaxation processes in compressible multiphase flows[END_REF]. In the present work, ν is considered very large, so that relaxation to thermodynamic equilibrium is immediate.

During the phase transition process, the mixture specific volume v = 1/ρ and energy e do not vary. The mass fractions for all gas species other than the vapor (Y k≥3 ) also remain constant.

However, pressure and temperature do vary, reaching their equilibrium values (p * , T * ). The phase transition model thus reduces to computing the equilibrium state (p * , T * , Y * k ), at every time step, from the state described by (p, T, v, e, Y k ).

Since pressure and temperature are functions of (v, e, Y k ), and v and e are constant during the phase transition, the primary goal of the procedure is to compute accurately Y * 1 , since Y * 2 is linked through mass conservation:

Y * 2 = 1 -Y * 1 - k≥3 Y k .
In our approach, phase transition is decoupled of transport and wave propagation using a fractional step method, essentially reducing the model to setting Y 1 = Y * 1 and Y 2 = Y * 2 after each hyperbolic step (resolution of System (1.4) without mass transfer).

Expression of the thermochemical equilibrium

When only liquid and vapor are present (Y k≥3 = 0), as in our previous paper [START_REF] Chiapolino | A simple phase transition relaxation solver for liquidvapor flows[END_REF], the thermo-

chemical equilibrium directly translates into p * = p sat (T * ), independently of Y * 1 .
Here an extra complexity is added, as the vapor partial pressure in the gas phase is equal to the saturation pressure at the current temperature.

The relation between p * and T * then depends on the composition of the multicomponent gas Y k≥2 . Within ideal gas mixtures, the vapor partial pressure is directly proportional to the vapor molar fraction:

p partial = Y 2 /W 2 k≥2 Y k /W k p, (3.2)
so that the expression satisfied at thermochemical equilibrium is

p sat (T * ) = p partial = Y * 2 /W 2 k≥2 Y * k /W k p * . (3.3)
This relation is indeed essential as it allows equilibrium of liquid and multicomponent gas in conditions below saturation (for instance, existence of water vapor at ambient temperature and pressure). Case 2 T < T sat (p): the mixture can be purely gaseous when the vapor quantity does not exceed its saturation limit, to be established hereafter. If this value is reached, then liquid is present.

Solution of a simplified problem

The attainable limit for the vapor mass fraction before condensation in the multicomponent gas is obtained by solving p partial = p sat (T ). Using Eq. (3.2), one gets

Y sat 2 = p sat (T )W 2 p -p sat (T ) k≥3 Y k /W k . (3.4)
Below the saturation temperature, we then have

Y * 2 = min Y sat 2 , 1 - k≥3 Y k . (3.5)
The equilibrium mass fraction Y * 2 cannot exceed the available room (1 -k≥3 Y k ), nor can it exceed the saturation limit, after which liquid is present. Y * 1 is then deduced from mass conservation. Figure 2 shows the solution for Y * 1 and Y * 2 of this simple problem, for a mixture of liquid water, water vapor and air (N = 3), at atmospheric pressure, and a mixture temperature of 350 K. On the right-hand side of the plot, there is so much air content that water vapor is diluted enough (below the saturation limit), and no liquid is present. On the left part, Y * This simplified problem is convenient for two reasons: it illustrates well the problem to be solved, and also provides a simple way to compute initial conditions in a simulation.

Y 3 = Y air Y * 1 = Y liq Y * 2 = Y vap Y sat 2 mass fractions v = 1/ρ (m 3 /kg)
When the problem is to be solved at constant (v, e), instead of (p, T ), an analytic approach is no longer possible, and the problem has to be solved numerically. In the following section, we present an alternative to a tedious iterative process.

Thermochemical relaxation algorithm

The idea of our method is to gradually reach the exact solution (typically in 1 to 3 time steps of the flow solver), by providing a fair approximation for Y * 1 , while iterative approaches, such as the one promoted in Le [START_REF] Métayer | Dynamic relaxation processes in compressible multiphase flows[END_REF] [START_REF] Métayer | Dynamic relaxation processes in compressible multiphase flows[END_REF], directly computes the exact solution. In that direction, we will follow a similar strategy than in our previous work [START_REF] Chiapolino | A simple phase transition relaxation solver for liquidvapor flows[END_REF], extended to the fact that the relation between pressure and saturation pressure is modified compared to (3.3), and now depends on the result (since the partial pressure is a function of Y 2 ). First, let us bound the equilibrium mass fractions Y * 1 and Y * 2 :

       Y min = ǫ, Y max = 1 -Y min - N k=3 Y k . (4.1)
With ǫ → 0 for numerical purposes, typically on the order of 10 -8 . As in our previous work [START_REF] Chiapolino | A simple phase transition relaxation solver for liquidvapor flows[END_REF], we first assume Y * 1 = Y min , to check if the system has a solution without liquid. Under this assumption, pressure and temperature are computed by use of the mixture equation of state Note that, unlike in our last contribution [START_REF] Chiapolino | A simple phase transition relaxation solver for liquidvapor flows[END_REF], there is no need to check the existence of the pure liquid phase: no matter how small ǫ is chosen, there will always be a (very small) solution for Y * 2 satisfying that its partial pressure is equal to the saturation pressure. Section 5 (about the algorithm's stability) provides additional comments and clarify this last feature. If Y * 1 = Y min , 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 * ), (4.2)
leaving us with three equations for three unknowns (p * , T * , Y * 1 ). The 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 . ( 4.3)
We will now evaluate three approximate expressions for Y * 1 . Rewriting the specific volume definition from System (4.2) leads to the following approximate expression: Equations for Y m 2 (p, T ) and Y e 2 (p, T ) (expressed for the vapor phase instead of the liquid phase) can be obtained as,

Y m 1 (p, T ) = v -v g (p, T ) v 1 (p, T ) -v 2 (p, T ) , with v g (p, T ) = 1 - N k=3 Y k v 2 (p, T ) + N k=3 Y k v k (p, T ).
Y 2 = 1 -Y 1 - N k=3 Y k .
In our previous work [START_REF] Chiapolino | A simple phase transition relaxation solver for liquidvapor flows[END_REF], these two functions of (p, T ) could be reduced to functions of p in finding the equilibrium as p = p sat (T ) for a liquid-vapor mixture. As shown in Eq.

Y sat 2 (p, T ) = p sat (T )W 2 p -p sat (T ) k≥3 Y k /W k . (4.6)
The exact solution to the problem is found in pressure and temperature conditions satisfying all three Eqs. of (4.2), or, equivalently:

Y m 2 (p * , T * ) = Y e 2 (p * , T * ) = Y sat 2 (p * , T * ).
This is illustrated in Fig. 3, for a set of conditions initially out of equilibrium (30K below saturation).

In this graph, -Y m = Y m 2 (p, T ) is evaluated for the initial values of (p, x v ) (from the hyperbolic step), and

-
T = T sat (x v .p), -Y e = Y e
2 (p, T ) is evaluated for the initial values of (p, x v ) (from the hyperbolic step), and

T = T sat (x v .p), -Y sat = Y sat 2 (p, T ) is evaluated at the initial (p, T ).
As in [START_REF] Chiapolino | A simple phase transition relaxation solver for liquidvapor flows[END_REF], the idea of this method is then to pick the value with the smallest variation, under the condition that all three associated mass transfers are of the same sign. This idea is reminiscent of the well known Minmod slope limiter widely used in high order CFD solvers for gradients limitation. In other words, we introduce: where Y 2 is the initial mass fraction (from the hyperbolic step), and Y m , Y e , Y sat are three estimates for the equilibrium Y * 2 .

     r 1 = (Y m -Y 2 ) (Y e -Y 2 ) , r 2 = (Y m -Y 2 ) (Y sat -Y 2 ) , (4.7)
-If r 1 < 0, or r 2 < 0, no mass transfer happens:

Y * 2 = Y 2 -Else, the minimum mass transfer is used. For instance, if Y m is the closest value to Y 2 (out of Y m , Y e , Y sat ), then Y * 2 = Y m .
Under the latter condition, it can be implemented as:

Y * 2 = Y 2 + sgn [Y m -Y 2 ] × Min |Y m -Y 2 |, |Y e -Y 2 |, |Y sat -Y 2 | (4.8)
Summary of the procedure:

1) Search for a solution without liquid -The limit case Y * l → ǫ is considered with computation of the pressure and temperature according to the mixture equation of state (2.8) and (2.10). Additionally the corresponding saturation pressure and partial pressure for vapor are computed according to Eqs. (2.4) and (3.2).

-If the partial pressure is below the saturation pressure, no liquid is to be considered and the solution is Y * 1 = Y min and Y 2 = Y max according to Eq. (4.1).

2) If liquid is present, thermodynamic equilibrium is to be computed -Computation of the liquid mass fraction at equilibrium Y * 1 is done according to the Minmod-like procedure Eq. (4.8) using Eqs. Y m (4.4), Y e (4.5), Y sat (4.6), (4.7), and the EOS (2.7).

Algorithm's stability

In a previous contribution [START_REF] Chiapolino | A simple phase transition relaxation solver for liquidvapor flows[END_REF] in which only liquid and vapor were considered, a major challenge was to deal with the pure liquid limit. Let us consider pure liquid water, initially at 300 K and 1 atm. If a small quantity ǫ of vapor appears, the system has to go to thermochemical equilibrium, no matter how far off the initial conditions are. Because of the high heat capacity and low compressibility of liquid water, the process is mostly isothermal, and important pressure fluctuations could be observed. This motivated the introduction of a correction on the pressure used to estimate Y m and Y e , in order to improve the method capabilities to deal with the transition from pure liquid to two-phase mixture in the limit Y k≥2 → 0.

Here, there is no need for such a correction. The situation detailed above, in which a mixture can be found very far off equilibrium after the hyperbolic step, can not be encountered. In the present work's context, we can find a value for the vapor mass fraction Y 2 = Y vap , no matter the conditions given, as traces of incondensable gas are present Y 3 = Y air . At 300 K and 1 atm for a mixture of liquid water, vapor water and air, assuming Y 3 = ǫ (which can be as small as numerically allowed) will lead to a very small but non zero Y 2 . With these considerations, the lack of continuity in the transition from pure liquid to two-phase mixture, present in the liquid-vapor model, vanishes in the present context.

Additionally, let us note that presence of air in water is physical (dissolved), and can be estimated through Henry's law. For instance, in ambient conditions, the air mass fraction dissolved is on the order of Y 3 = 10 -5 , leading to a value of Y 2 of the same order, satisfying System (4.2). Indeed, such low values do not have significant effect on mean density and energy, which remain close to the pure phase values.

In practical computations presented in this paper, unless stated otherwise, a small value is set for Y 3 in the "pure" liquid, and Y 2 is computed accordingly, following Eq. (3.5) for the initial prescribed pressure and temperature.

Another key property responsible for the algorithm's stability is that ∀(p, T ) (a priori away from the equilibrium point), the exact solution for Y * 2 always lies between the minimum and maximum values of Y m 2 (p, T ), Y e 2 (p, T ), Y sat 2 (p, T ) This is clearly visible in Fig. 3. Although we have not proved this result because of the non-linear dependence of p sat (T ) (as in our previous contribution), we have not found a case in which this does not apply. As a consequence:

-If Y 2 > Y *
2 initially, mass transfer will be activated only if Y m , Y e and Y sat are all smaller than Y 2 . Given the above property, one of them at least is superior to Y * 2 . The algorithm then cannot "overestimate" the mass transfer.

-The same thing can be said when Y 2 < Y * 2 initially.

-Last but not least, since 0 < Y 2 < 1 -k≥3 Y k (given the hyperbolic step is properly implemented), and the exact solution also satisfies 0 < Y * 2 < 1 -k≥3 Y k , the approximate solution from the algorithm can never go out of bounds.

Numerical results

Robustness and accuracy of the present thermochemical relaxation algorithm have been illustrated in a previous contribution [START_REF] Chiapolino | A simple phase transition relaxation solver for liquidvapor flows[END_REF] for two-phase flows containing a liquid and its own vapor only. In the current context, non-condensable gases are present. The proposed relaxation solver is to be compared to an iterative algorithm such as the one promoted in Le Métayer et al.

(2013) [START_REF] Métayer | Dynamic relaxation processes in compressible multiphase flows[END_REF]. Unlike this last algorithm, the notion of vapor molar fraction and partial pressures in the gas phase must be taken into account, bringing additional difficulties. The corresponding iterative solver is consequently recalled and modified in Appendix.

To illustrate robustness and accuracy of the relaxation algorithm, a shock tube containing liquid water, its own vapor and some non-condensable air is considered, with variable initial conditions. Shock tube tests appear as excellent benchmarks as the flow contains shock waves, contact discontinuities and rarefaction fans that create some arduous conditions. In this section, the tube is 1 meter long and the initial discontinuity is located at 0.5 meter. Boundary conditions are non-reflecting ones. The computations that follow are addressed with the MUSCL Hancock method using van Leer's slope limiter and the HLLC Riemann solver (see [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: A practical introduction[END_REF] [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: A practical introduction[END_REF], [START_REF] Saurel | A general formulation for cavitating, boiling and evaporating flows[END_REF] [START_REF] Saurel | A general formulation for cavitating, boiling and evaporating flows[END_REF] for details). This section compares the results obtained with the present relaxation solver and the classic root-finding procedure method given in Appendix. Unless stated otherwise, the one-dimensional simulations are carried out on meshes with 100 cells.

Shock tube test with a mixture far from the phase bounds

A two-phase mixture with initial mass fractions set to Y 1 = 0.1 (liquid), Y 2 = 0.2 (vapor) and Y 3 = 0.7 (air) is considered throughout the entire tube with an initial pressure ratio of 2, resulting in the presence of initial density and temperature discontinuities.

In order to find an initial thermodynamic equilibrium, the temperatures of both chambers are deduced as

T = T sat (x v p) with x v = Y 2 /W 2 Y 2 /W 2 + N k=3 Y k /W k . (6.1)
The mixture density is then computed as,

ρ = 1 v with v = N k=1 Y k v k (p, T ), (6.2)
where v k (p, T ) is given by Eq. (2.7). The results are shown at time t ≈ 1 ms in Fig. 4. Perfect agreement is obtained between the present relaxation solver and the non-linear iterative method. The shock compression yields evaporation and the expansion wave results in condensation in this example.

x (m) x (m) ρ (kg.m -3 ) p (MPa) u (m.s -1 ) T (K) Y 1 Y 2

Shock tube test with a mixture with air in major proportions

Mass fractions of the previous test were given as inputs and the temperatures were deduced in order to respect initial thermodynamic equilibrium. Let us now consider a different point of view. As it is easier to measure or impose temperature in a given flow, it is now used as initial input. Pressures, velocities and the mass fraction of the non-condensable gas (air) remain inputs as well. Proportions of liquid and vapor are then deduced according to Eq. (3.5).

The following test considers a shock tube with an initial pressure ratio of Perfect agreement is obtained between the two solvers. The shock compression results in total evaporation of the liquid water and the rarefaction results in condensation. Disappearance of liquid water is computed without oscillations by the two methods.

Shock tube test with a mixture mainly made of water

Mass fractions from the previous test case are reversed by considering Y 3 = 10 -5 and T = 293 K initially throughout the tube, corresponding to a subcooled liquid. The following test considers a shock tube with an initial pressure ratio of 2. The liquid and vapor mass fractions are then deduced as Y 1 ≃ 0.9999899 and Y 2 ≃ 9.426 • 10 -8 in the left chamber and Y 1 ≃ 0.99989 and Y 2 ≃ 1.919 • 10 -7 in the right chamber. The results are shown at time t ≈ 1.5 ms in Fig. 6.

Perfect agreement is again obtained between the two solvers. The liquid phase having the highest internal heat capacity, the resulting flow is quasi-isothermal despite the slight evaporation and condensation processes appearing through the shock and rarefaction.

x (m) x (m) ρ (kg.m -3 ) p (MPa) u (m.s -1 ) T (K) Y 1 Y 2
u l = u r = 0 m.s -1 , Y lef t 1 ≃ 1.073 • 10 -2 , Y right 1 ≃ 1.186 • 10 -3 , Y lef t 2 ≃ 9.263 • 10 -3 , Y right 2 ≃ 1.881 • 10 -2 and Y lef t 3 = Y right 3 = 0.98.
Final time: t ≈ 1 ms. Mesh: 100 cells. For the sake of clarity, only 50 symbols out of 100 are plotted for the iterative method. Full liquid evaporation is correctly computed by both methods. x (m)

x (m) ρ (kg.m -3 ) p (MPa) u (m.s -1 ) T (K) Y 1 Y 2
u l = u r = 0 m.s -1 , Y lef t 1 ≃ 0.9999899, Y right 1 ≃ 0.99989, Y lef t 2 ≃ 9.45 • 10 -8 , Y right 2 ≃ 1.919 • 10 -7 and Y lef t 3 = Y right 3 = 10 -5
. Final time: t ≈ 1.5 ms. Mesh: 100 cells. For the sake of clarity, only 50 symbols out of 100 are plotted for the iterative method.

Double expansion test with a two-phase mixture mainly made of liquid water

The following test mimics cavitation in an initial subcooled liquid. The initial pressure, temperature and air mass fractions are respectively 1 bar, T = 293 K and Y 3 = 10 -5 throughout the entire tube. The liquid and vapor mass fractions are then deduced as Y 1 ≃ 0.99989 and Y 2 ≃ 1.919 • 10 -7 in the left and right chambers. The initial velocity is set to -1 m.s -1 at left and +1 m.s -1 at right. The results are shown at time t ≈ 3.5 ms in Fig. 7. x (m) x (m) Perfect agreement is again obtained between the two solvers. Appearance of vapor in the liquid is computed without oscillations by the two methods. model developed is the only possible term for vapor production. The Stiffened Gas EOS is used, the parameters for oxygen and hydrogen are given in Table 2 The 2D computations are carried out on a simplified geometry given in Fig. 9. The inlet boundary conditions correspond to two subsonic inflows:

ρ (kg.m -3 ) p (MPa) u (m.s -1 ) T (K) Y 1 Y 2
-central flow made of nearly pure liquid oxygen, at 100 K and 30 m.s -1 , pressure 3 MPa along segment AB of Fig. 9, -peripheral flow made of nearly pure gaseous hydrogen, at 150 K and 200 m.s -1 , with the pressure of 3 MPa along segment DE of Fig. 9, and one subsonic outlet, along segments GH and HI, where non-reflecting subsonic boundary condition is imposed at 3 MPa. The remaining walls are treated as symmetric boundary conditions. The mesh consists of 360, 000 triangles, with a spatial discretization varying from 0.1 mm in the jet to 0.6 mm in the far field.

The small scale destabilization of the liquid jet requires to extend the numerical solver presented in [START_REF] Saurel | A general formulation for cavitating, boiling and evaporating flows[END_REF] [START_REF] Saurel | A general formulation for cavitating, boiling and evaporating flows[END_REF] to second order. This was achieved following the MUSCL-Hancock scheme with the Superbee limiter, as detailed for example in Toro [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: A practical introduction[END_REF].

Two computations have been carried out on the same mesh, with the aim of studying the effect of the thermochemical relaxation on jet destabilization. The contours of mass fraction of liquid oxygen as obtained with and without evaporation effects are compared in Figs. [START_REF] Lallemand | Pressure relaxation procedures for multiphase compressible flows[END_REF].a and 10.b, showing serious changes due to phase transition.

As expected, the filaments separating the main liquid core and the gas gradually vanish as a consequence of evaporation, resulting in much steeper contours of liquid mass fraction in the vicinity of the jet, and very few pockets of liquid with significant life time. Fig. 10.d shows the contour of vapor oxygen (which indeed remains at zero in the case without phase transfer). Fig. 10.c plots the total oxygen mass fraction contour (Y 1 + Y 2 ). It can also be compared to Fig. 10.a

Although the results in terms of total mass fraction for oxygen qualitatively seem close whether or not the mass transfer is activated, it is clear that the vapor mass fraction is of utmost importance for our future work, which shall include the gaseous combustion between vapor oxygen and hydrogen.

Capillary effects are not included in these simulations, as the intense velocity gradients make them negligible in the present jet destabilization conditions. However, capillary effects may be of importance when the filaments start separating from the jet, and our future work will include them, following the approach presented in Le Martelot et al (2014) [START_REF] Le Martelot | Towards the direct numerical simulation of nucleate boiling flows[END_REF].

Explosive liquid water dispersal into air

The present relaxation solver has been presented in the context of the 4-equation formulation.

This model considers mechanical and thermal equilibrium but thermochemical disequilibrium.

However the relaxation solver is not restricted to such a model and can be extended to nonequilibrium hyperbolic flow models, such as Baer and Nunziato's (1986) [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] and its reduced versions. This last model is a 7-equation formulation assuming mechanical, thermal and chemical disequilibrium. Thereby, finding local thermodynamic equilibrium requires extra ingredients such as velocity, pressure and temperature relaxation solvers such as the ones given in [START_REF] Métayer | Dynamic relaxation processes in compressible multiphase flows[END_REF][START_REF] Lallemand | Pressure relaxation procedures for multiphase compressible flows[END_REF].

We propose here to repeat the explosive liquid water dispersal test presented in Furfaro and Saurel (2015) [START_REF] Furfaro | A simple HLLC-type Riemann solver for compressible non-equilibrium two-phase flows[END_REF] with the 4 and 7-equation models. The test consists in a cylindrical gas-liquid explosion schematically represented in Fig 11 . The results are given in Fig 12.

Nearly pure gas (1 kg/m 3 ) Atmospheric pressure (0.1 Mpa)

Nearly pure liquid water (1000 kg/m 3 ) Atmospheric pressure (0.1 Mpa)

Nearly pure dense gas (1500 kg/m 3 ) HP chamber (7000 Mpa)

Figure 11: Schematic representation of the cylindrical gas-liquid explosive system. A cylindrical explosive charge is surrounded by a liquid water layer. The internal cylinder is 160 mm long with a 52 mm radius and is initially filled with nearly pure dense gas (1500 kg/m 3 ) at high pressure (7000 Mpa). The external cylinder is 160 mm long with a 96 mm radius and is initially filled with nearly pure liquid water (1000 kg/m 3 ) at atmospheric pressure (0.1 MPa). Atmospheric conditions with nearly pure gas (1 kg/m 3 ) at 0.1 MPa are considered around both cylinders.

No noticeable changes are visible in

Conclusions

A simple relaxation solver has been built, able to deal with phase transition between a liquid phase, and a multicomponent gas phase. Several features are in favor to the present solver compared to iterative ones:

-CPU savings over 50%, -high simplicity of implementation and verification, AppendixA. Summary of the iterative algorithm

The algebraic system to be solved is,

                     v = Y * 1 v 1 (p * , T * ) + Y * 2 v 2 (p * , T * ) + k 3 Y k v k (p * , T * ), e = Y * 1 e 1 (p * , T * ) + Y * 2 e 2 (p * , T * ) + k 3
Y k e k (p * , T * ),

p sat (T * ) = Y * 2 W 2 Y * 2 W 2 + k 3 Y k W k p * , (A.1)
where the superscript * denotes the thermodynamic equilibrium state. Since Y 1 = 1 -Y 2k 3 Y k , the algebraic system reads, It is convenient to use a reference state defined with the other NASG or SG coefficients which have been determined via theoretical equations and experimental data (see Le Métayer et al.

                       f 1 (p * , T * , Y * 2 ) = v -1 -Y * 2 -
[11], [START_REF] Métayer | The Noble-Abel Stiffened-Gas equation of state[END_REF]). Let us define, (A.6)

       T ref = q 1 -v 2 C p,2 -C v,2 = B C , p ref = (C p,2 -C v,2 ) p ∞,1 -(C p,1 -C v,1 ) p ∞,2 (C p,1 -C v,1 ) -(C p,2 -C v,2 ) = p ∞,
The algebraic system thus reads, Y k e k (p * , T * ), The system F (X) = 0 has to be resolved. Its solution is given by,

                                   f 1 (p * , T * , Y * 2 ) = v -1 -Y * 2 - k 3 Y k v 1 (p * , T * ) -Y * 2 v 2 (p * , T * ) - k 3 Y k v k (p * , T * ),
f 3 (p * , p * sat , Y * 2 ) = p * sat - Y * 2 W 2 Y * 2 W 2 + k 3 Y k W k p * ,
X n = X n-1 -J F (X n-1 ) -1
F (X n-1 ), (A.9)

with n denoting the current iteration and J F (X n-1 )

-1

the inverse of the Jacobian matrix, In order to gain accuracy and robustness, it is convenient to write Eq. (A.9) as,    J(F (X n-1 ))△X = -F (X n-1 )), △X = X n -X n-1 .

J(f 1 , f 2 , f 3 , f 4 ) =       
(A.11) Thereby, for a given X, a linear system is found and the variation △X is determined via the Gauss elimination method. The variables are then updated as, X n = △X + X n-1 .

(A.12)

A solution is obtained when △X < ǫ or alternatively when |f This procedure is robust and accurate. It is used in all test problems of Section 6.

  species occupies the entire volume p = N k=2 p partial, k : the mixture pressure is the sum of partial pressures e = N k=2 Y k e k : mixture internal energy definition

  the total specific volume is the sum of fluids' specific subvolumes p = p k : pressure equilibrium among the chemical species e = N k=2 Y k e k : mixture internal energy definition

Figure 1 :

 1 Figure 1: Comparison between experimental and theoretical saturation curves for liquid water and steam with coefficients determined in the temperature range [273 -500 K]. The symbols represent the experimental data. The thick lines represent the NASG theoretical saturation curve. P sat denotes the saturation pressure, L v the latent heat, h the specific enthalpy and v the specific volume.

  Instead of solving thermochemical equilibrium at constant v * = v, e * = e and Y * k≥3 = Y k≥3 , and computing (p * , T * , Y * 1 ), let us assume p * = p, T * = T and Y * k≥3 = Y k≥3 are constant, and v and e are varying. The thermochemical equilibrium problem is straightforward to solve analytically in these conditions: Case 1 T > T sat (p): the mixture is necessarily purely gaseous. We then have Y * 1 = 0.

2

  varies linearly with Y 3 , as a consequence of the saturation relation p partial = p sat (T ). The figure also indicates the variation of the mixture specific volume, which increases up to the saturation limit (since the liquid fraction diminishes), and then decreases, because air is heavier than water vapor.

Figure 2 :

 2 Figure 2: Equilibrium mass fractions obtained from water liquid/vapor, as a function of the air mass fraction. Conditions: atmospheric pressure, T=350K. dashed line: Y * 1 , dashed-dotted: Y * 2 . The thick line represents the specific volume obtained for the mixture in these conditions.

p

  = p(v, e, Y ), T = T (v, e, Y ), since the mixture variables v and e are invariant through phase transition. Additionally, we compute the corresponding partial pressure for vapor (3.2), and compare it to the saturation pressure (2.4). If the partial pressure is below the saturation pressure, no liquid is present and the solution is Y * 1 = Y min and Y * 2 = Y max .

(4. 4 )

 4 Any mass fraction computed from this formula after the hyperbolic step will satisfy the mass conservation ∀(p, T ). Similarly, we can rewrite the specific mixture energy definition from System (4.2): Y e 1 (p, T ) = e -e g (p, T ) e 1 (p, T ) -e 2 (p, T ) , with e g (p, T ) = 1 -N k=3 Y k e 2 (p, T ) + N k=3 Y k e k (p, T ).(4.5)

  (3.3), the new relation p partial = p sat (T * ) includes an additional dependence on Y * 2 (or, equivalently, on Y * 1 ), leading to Eq. (3.3), reminded here:

  the light gray surface represents the mass conservation condition: Y m 2 (p, T ), -the medium gray surface represents the energy conservation condition: Y e 2 (p, T ), -the dark gray surface, represents the thermochemical equilibrium condition: Y sat 2 (p, T ).The algorithm consists in computing three approximations for Y * 2 :

Figure 3 : 2 (

 32 Figure 3: Illustration of the algorithm for the multicomponent phase transition solver. From light gray to dark gray: Y m (p * , T * ), Y e (p * , T * ), Y sat (p * , T * ). The exact solution of the problem is found at the intersection of the three surfaces. The value from hyperbolic step corresponds to the initial conditions of the algorithm: Y 2 = 2.Y * 2 (far from equilibrium) at p = 1 atm, T = T sat (p) -30 K = 343 K, and Y 3 = 0.1. Here, computations lead to Y * 2 = 0.1, Y m = 0.11, Y e = 0.9, Y sat = 0.03, so that the algorithm retains Y m , within 10% of the exact solution.

Figure 4 : 1 = Y right 1 = 0 . 1 , Y lef t 2 = Y right 2 = 0 . 2 and Y lef t 3 = Y right 3 = 0 . 7 .

 4110122023307 Figure 4: Comparison of the present relaxation algorithm (thick lines) versus the non-linear iterative solver (symbols). Shock tube test with a two-phase mixture made of liquid water, vapor water and air. The dotted lines represent the solutions without phase transition. The dashed lines represent the initial conditions: p l = 2•10 5 Pa, p r = 10 5 Pa, u l = u r = 0 m.s -1 , Y lef t 1

  2 and initial temperature of T = 293 K throughout the tube. The air mass fraction is initially set to Y 3 = 0.98 in the whole tube. The liquid and vapor mass fractions are then deduced as Y 1 ≃ 1.073 • 10 -2 and Y 2 ≃ 9.263 • 10 -3 in the left chamber and Y 1 ≃ 1.186 • 10 -3 and Y 2 ≃ 1.881 • 10 -2 in the right chamber. The results are shown at time t ≈ 1 ms in Fig. 5.

Figure 5 :

 5 Figure 5: Comparison of the present relaxation algorithm (thick lines) versus the non-linear iterative solver (symbols). Shock tube test with a two-phase mixture with air in major proportions. The dotted lines represent the solutions without phase transition. The dashed lines represent the initial conditions: p l = 2 • 10 5 Pa, p r = 10 5 Pa, u l = u r = 0 m.s -1 , Y lef t

Figure 6 :

 6 Figure 6: Comparison of the present relaxation algorithm (thick lines) versus the non-linear iterative solver (symbols). Shock tube test with a two-phase mixture mainly made of liquid water. The dotted lines represent the solutions without phase transition. The dashed lines represent the initial conditions: p l = 2 • 10 5 Pa, p r = 10 5 Pa, u l = u r = 0 m.s -1 , Y lef t

Figure 7 : 1 = Y right 1 ≃ 0.99989, Y right 2 = Y right 2 ≃ 1 = 10 - 4 .

 711221104 Figure 7: Comparison of the present relaxation algorithm (thick lines) versus the non-linear iterative solver (symbols). Double expansion test with a mixture mainly made of liquid water. The dotted lines represent the solutions without phase transition. The dashed lines represent the initial conditions: p l = 2 • 10 5 Pa, p r = 10 5 Pa, u l = -1 m.s -1 , u r = +1 m.s -1 , Y lef t 1

1 Figure 8 : 1 ≃ 1 . 1 ≃ 1 .186 • 10 - 3 , Y lef t 2 ≃ 9 . 2 ≃ 1 . 3 = Y right 3 =

 181111103292133 Figure 8: Comparison of the present relaxation algorithm (thick lines) versus the non-linear iterative solver (symbols). Shock tube test with a two-phase mixture mainly made of air, computed on coarse (left figure) and fine (right figure) meshes to assess method robustness. The dashed lines represent the initial conditions: p l = 2 • 10 5 Pa, p r = 10 5 Pa, u l = u r = 0 m.s -1 , Y lef t 1

XFigure 9 :

 9 Figure 9: Geometrical data for half the computational domain of the evaporating liquid jet.

a) without phase change: Y 1 b 2 Figure 10 :

 1210 Figure 10: Fragmentation of a liquid O 2 jet by a coaxial high-speed H 2 current. The results from Figs b,c and d are from the same computation, whereas Fig a is taken at the same time, from a computation without phase transition. All results present mass fraction contours, on the same color map. Figs a) and b) compare the liquid mass fraction of oxygen without and with phase transfer. Fig. d) shows the contour of vapor oxygen resulting from the liquid jet evaporation. Fig. c) represents the sum Y 1 + Y 2 (Fig b+c), showing a total mass fraction of oxygen, to be compared to the case without mass transfer (Fig. a).

  Fig 12 illustrating the capabilities of the relaxation solver to model phase change in both 7-equation and 4-equation models . An evaporation front is clearly visible with both models.

Figure 12 :

 12 Figure 12: Vapor mass fraction created during a 2D cylindrical gas-liquid explosive dispersal. The figure compares the results obtained with the 4-equation model (mechanical and thermal equilibrium, top) and with the 7-equation model (total disequilibrium, bottom) after relaxation of the pressures, velocities, temperatures following the relaxations solvers presented in [9] and phase transition consideration through the relaxation method presented in this paper. The results computed by both models are in close agreement.

k 3 Y 3 Y- k 3 Y k e 1 3 Y 2 )

 333132 k v 1 (p * , T * ) -Y * 2 v 2 (p * , T * )k k v k (p * , T * ), f 2 (p * , T * , Y * 2 ) = e -1 -Y * 2 (p * , T * ) -Y * 2 e 2 (p * , T * )k k e k (p * , T * ), f 3 (p * , p * sat , Y * 2 ) = p * sat -However the relation between p sat (T * ) and T * ,ln (p sat + p ∞,v ) = A + B + E p sat T sat + C ln (T sat ) + D ln (p sat + p ∞,l ) , (A.3)is non-linear and no analytical relation can be found. An additional equation is consequently added to System (A.2).

f 2 - k 3 Y

 23 (p * , T * , Y * 2 ) = e -1 -Y * 2 k e 1 (p * , T * ) -Y * 2 e 2 (p * , T * )k 3

f 4 7 ) 1

 471 (p * sat , T * ) = p * sat + p ∞,2 p ref + p ∞,2 -Four variables are unknown, p * , T * , Y * v , p * sat and System (A.7) is non-linear calling for an iterative method such as Newton-Raphson's. Let us then denote, (p* , T * , Y * 2 ) f 2 (p * , T * , Y * 2 ) f 3 (p * , p * sat , Y * 2 ) f 4 (p * sat , T * )

1 |

 1 < ǫ, |f 2 | < ǫ, |f 3 | < ǫ, |f 4 | < ǫ, with ǫ → 0.

Table 1 :

 1 For liquid water and steam, the NASG fluid parameters, determined in the [273-500K] temperature range are given in Table1. In the same table, the coefficients for air 2 are given (in the frame of the ideal gas assumption). Noble Abel Stiffened Gas (NASG) coefficients for water and air determined in the temperature range [273 -500 K] These parameters are used in the computational examples of the present paper. A comparison between the NASG EOS and experimental data is displayed in Fig 1, showing good agreement.

	Coefficients	Liquid water Water vapor air
	C p (J/kg/K)	4285	1401	1007
	C v (J/kg/K)	3610	955	719
	γ	1.19	1.47	1.4
	P ∞ (Pa)	7028 10 5	0	0
	q (J/kg) q ′ (J/kg/K)	-1177788 0	2077616 14317	0 0
	b (m 3 /kg)	6.61 10 -4	0	0
	W (g/mol)	18	18	29

Table 2 :

 2 . Stiffened Gas coefficients for the simulation of the evaporating liquid jet.

	Coefficients	Liquid phase (O 2 ) Vapor phase (O 2 ) Hydrogen (H 2 )
	C p (J/kg/K)	1702	780	14256
	C v (J/kg/K)	695	531	10183
	γ	2.45	1.47	1.4
	P ∞ (bar)	1062	0	0
	q (J/kg) q ′ (J/kg/K)	-258000 0	6900 -9280	-1.2.10 6 0

  1 -Dp ∞,2 D -1 .Eq. (A.3) is now written with the built reference statep ref , T ref , ln (p ref + p ∞,2 ) = A + B + E p ref T ref + C ln (T ref ) + D ln (p ref + p ∞,1 ) . (A.5)38Subtracting Eq. (A.5) to Eq. (A.3), it is straightforward to find another formulation connecting the saturation pressure and temperature,

	p sat + p ∞,2 p ref + p ∞,2	= exp C	T ref T sat	+ ln	T sat T ref	-1	p sat + p ∞,1 p ref + p ∞,1	D	exp E	p sat T sat	-	p ref T ref	.
													(A.4)

   ∂f 1 (p * ,T * ,Y * 2 ) ∂p * ∂f 1 (p * ,T * ,Y * 2 ) ∂T * ∂f 1 (p * ,T * ,Y *

			2 )	∂f 1 (p * ,T * ,Y * 2 )	
	∂f 2 (p * ,T * ,Y * 2 )	∂f 2 (p * ,T * ,Y * 2 )	∂Y * 2 ∂f 2 (p * ,T * ,Y * 2 )	∂p * sat ∂f 2 (p * ,T * ,Y * 2 )	 
	∂p * ∂f 3 (p * ,p * sat ,Y * 2 )	∂T * ∂f 3 (p * ,p * sat ,Y * 2 )	∂Y * 2 ∂f 3 (p * ,p * sat ,Y * 2 )	∂p * sat ∂f 3 (p * ,p * sat ,Y * 2 )	  	(A.10)
	∂p *	∂T *	∂Y * 2	∂p * sat	 
	∂f 4 (p * sat ,T * )	∂f 4 (p * sat ,T * )	∂f 4 (p * sat ,T * )	∂f 4 (p * sat ,T * )
	∂p *	∂T *	∂Y * v	∂p * sat

Note that we only consider mass transfer between liquid and vapor, so the reference energies of the other gas components have no importance. If other mass transfers are to be considered (between gas species), reference energies have to be set appropriately.

Acknowledgements

Part of this work has been carried out in the framework of the Labex MEC (ANR-10-LABX-0092) and of the A*MIDEX project (ANR-11-IDEX-0001-02), funded by the Investissements d'Avenir French Government program managed by the French National Research Agency (ANR). We also acknowledge funding from ANR through project ANR-14-CE22-0014. Damien Furfaro and François Fraysse are gratefully acknowledged for carrying out the simulations presented in Fig. 12.

-increased robustness with high-order methods. This is reached without any noticeable loss of accuracy. This research can be continued in many directions. Among them, consideration of supercritical fluids in the same theoretical frame appears important. Another perspective deals with the consideration of combustion effects within the gas phase.

Computational time, efficiency and simplicity

All computational examples confsidered in this paper led to the same observation: the present relaxation solver is much faster than the root-finding iterative algorithm given in Appendix. In our previous work [START_REF] Chiapolino | A simple phase transition relaxation solver for liquidvapor flows[END_REF] CPU savings from 5 to 50% were obtained, depending on the conditions.

Due to the increased complexity of the iterative algorithm for the multi-component version, the CPU saving is now 50% in all reported cases. The main argument for this time gain is related to its simplicity: the relaxation solver is direct whereas the iterative method requires solving a non-linear-algebraic system that may cause difficulties as a result of non-linearities. Besides, the root-finding iterative method requires the calculation of the saturation pressure via Eq. ( 2.4) at each iterative step, which itself requires an iterative method.

In practical applications, fine discretization is usually used on zones of interest. However, as the method is not iterative, and converges after some time steps, it is interesting to address its robustness and accuracy in situation with both coarse space and time resolutions. To this end, Fig 8 repeats the test case already presented in Fig 5, with a 10-cell mesh and first order Godunov numerical scheme. These results are very reasonable, illustrating the robustness of the method and its correct behavior even on coarse meshes. Results for the same test case are also given in Fig 8 for a higher-order solver (MUSCL-Hancock scheme with the Superbee limiter), on a 1000-cell mesh. Again, perfect agreement is found between the present relaxation solver and the non-linear iterative method. Robustness of the present method thus seems quite independent of mesh resolution and order of accuracy of the numerical scheme.

Multi-dimensional illustrations

Evaporating liquid jet

In this section, the capabilities of the flow model are illustrated on an evaporating liquid jet configuration in conditions typical of cryotechnic rocket engines during the ignition phase (for which the engine has not yet reached supercritical conditions). The flow consists of a coaxial liquid oxygen jet surrounded by a high-speed hydrogen flow, injected in conditions above the saturation point of the inner oxygen core, which then evaporates whilst being destabilized. Such a case is very challenging, because there is initially no vapor oxygen, and the phase transfer


Philippe Tigréat

Carlos Rosar

Kos Lassance

Xiaoran Jiang

Vincent Gripon

Philippe Tigreat

Xiaoran Liang

Claude Berrou

Assembly output codes for learning neural networks

Keywords: Assembly coding, Clustered Clique Networks, ECOC, Deep Learning, Coding theory, Classification

come

I. INTRODUCTION

Automatic learning systems are different from storing systems in that they aim at generalizing to unknown inputs. This happens through the extraction of the core features of learned data, and works as long as the unknown data to extrapolate to follows a similar distribution. The system thus learns a dictionary of features that is targeted to be well suited for the task at hand, e.g. classitication. These features are meant to correspond to the most relevant building blocks underlying the training inputs, that unseen data samples would likely also be made of.

In supervised Deep Learning networks, an output is calculated through a pipeline of vector-matrix products between input data and connection weights, intertwined with non-linear mathematical operators. Each level of the network is associated with a set of features that is more and more abstract as one moves towards the upper layers. During learning, an error is calculated from the difference between the resulting output and an objective vector specific to the class of the input. A gradient is then calculated from this error for the whole set of connection weights, and the hierarchy of features thus gets optimized through gradient descent for the task of classifying the data-set examples.

Combining classifiers has been an extensive area of research for a few decades [START_REF] Tulyakov | Review of classitier combination methods[END_REF] and several algorithms have been shown to bring improved decision by leveraging the diversity brought by an assembly of systems. Among these methods, boosting is a way to combine opmlOns of experts by weighting them based on their respective estimated accuracy. These base classifiers can be ditferentiated using various strategies, like training them on various subsets of the data or by providing them with ditferent sub-parts of an ensemble of learned feature detectors.

Another way to combine classifiers is to split the problem into a set of binary problems. A base classitier will then focus on classifying inputs between two of the initial classes, as in One-Vs-One (OVO), or between one class and the rest as in One-Vs-All (OVA). A way to implement these strategies is to provide a classitier with objective vectors that are not always specific of a single class but can be associated with a set of classes. Much attention is paid here to the Error-Correcting Output Coding (ECOC) method, which splits a multi-class problem into several two-way classification problems between meta-classes. It is shown in [START_REF] Bagheri | Error correcting output codes for multiclass classification: application to two image vision problems[END_REF] that ECOC can reach a better classification performance on an image data-set as compared to other multi-class methods.

The approaches presented here are derived from ECOC for multi class problems. These methods allocate assemblies of output neurons to the different input classes, with potential overlap between the codes of two classes. A clustering of the output layer is also applied, with only one active neuron per cluster for each target vector, and a local soft-max [3] process applied in each cluster at test time. The soft-max operator has the effect of normalizing to 1 the sum of energies of output neurons in each cluster. Experimental results suggest that the number of classes sharing an output node impacts performance, and so does the minimal distance between class codes. This finding is maintained when output codes are repeated so as to ensure that the ditferent tested networks have very similar numbers of parameters.

The outline of the paper is as follows. Section II provides theoretical considerations on the advantages of different output codes. Section III explains the methodology used in training the networks. Section IV presents experimental results.

II. CODING THEORY

Prior to experimenting with assembly codes as output for neural networks, a theoretical analysis can provide insights about which assemblies should perform better.

Consider a classifier with P classes to identify. The simplest way to make the classifier express its decision is to assign a single output node (the so-called grandmother cell) to each class. We propose to replace these P nodes with n= (�) nodes representing all the combinations (or assemblies) of m classes among P. Let us detine the coding rate R of IOg2( p) and n: log2( p)

the corresponding code as the ratio between R = (�) .

(1)

To calculate the minimum Hamming distance, let us consider any two classes among P. The number of assemblies that contain neither one nor the other is ( p � 2 ) and the number of assemblies that contain both is

( P-2
The product of R and dmin, called the merit factor F, is deduced as:

2 m (p -m ) log2 (p ) F = R dmin (P-l)P (3)
and its maximal value is obtained for m = P/2:

Plog2(P) F max 2( P_I) . ( 4
)
The corresponding minimum distance may be expressed as: clusters such that each class appears once and only once in each cluster. Therefore, there are' = P/m nodes in each cluster, assuming that P is a multiple of m (Fig. 1). This structure has two advantages.

( dmJmax = (p � 1 ) � (
Firstly, it complies exactly with the clustered clique-based associative memory proposed in [START_REF] Gripon | Sparse neural networks with large learning diversity[END_REF] which offers the possibility to store a number of patterns proportional to F (for larger sought diversities, the sparse scheme proposed in [START_REF] Kamary Aliabadi | Storing sparse messages in networks of neural cliques[END_REF] may be contemplated). However, the 2 286 optimal value I = 2 (as deemed by the optimization of F) is too low for a clique-based implementation and a trade-off has then to be found. For instance, , = 5 (that is, m = 2) seems a good choice for P = 10. The second advantage is that a cluster containing disjoint assemblies and therefore probabilities, through the soft-max principle [START_REF] Arribas | A model selection algorithm for a posteriori probability estimation with neural networks[END_REF], may be rigorously used for both learning and testing.

Ill. METHODOLOGY

Assembly codes

Couple cells are a method we introduce for encoding the class using a distributed code. Each cell no longer reacts to a single class but is specific of a couple of classes. Hence since the experiments are performed on data-sets with 10 classes, there are 45 possible couplings a given output neuron can be associated with. The whole set of 45 couplings is used here. Moreover, it is possible to partition these 45 couplings into 9 clusters of 5 couple cells where each individual class is represented exactly once in each cluster. This way the soft-max methodology can be applied inside of every cluster, between 5 competing hypotheses that are mutually exclusive.

Quintuplets are yet another code where each neuron is associated with a combination of classes, this time 5 among 10. The same approach as for the couples is used, by using all 252 possible quintuplets and partitioning them into 126 clusters each containing 2 complementary quintuplet cells.

The first two important factors considered to choose assembly codes are usability in a clique-based architecture and the merit factor. Alongside with the grandmother cell (our baseline), couples of classes (best trade-off cliques/merit factor) and quintuplets of classes (best merit factor, but larger output network) are tested. But only comparing these parameters is not enough, because one could argue that couples and quintuplets work better because they have a larger number of parameters for the last layer of the neural network. To avoid this we also train networks using repeated codes to reach output length of equal or comparable size. This repetition allows us to make a fair comparison between networks with virtually the same number of parameters in spite of using different output codes. Couples of classes have an interesting characteristic, in that there are plenty of different possible ways to partition 45 couple cells into 9 clusters, each featuring the 10 classes. Therefore it is possible to design an output layer of 1260 couple cells parted in 252 clusters where no cluster configuration is repeated twice.

During the training phase the categorical cross-entropy is used as the loss function. For classification a majority voting is done where each active output node votes for its associated set of classes. The assembly codes are summarized in table 1. ) and the last one is an output layer. At each connection layer, between the input and the activation, a batch normalization [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF] and a Gaussian noise with mean 0 and standard deviation 0. 3 are used. To finish the connection layer, a rectifier activation is used. The output layer has the length of the output code. A batch normalization is applied between the input and the activation and as in the case of the shallow network it uses a per-cluster soft-max 3 287 activation. The number of parameters of each network is given by the formula: 1000*(input length) + 750000 + 250*(n). This means: 1534000+250*(n) for MNIST and 3822000 + 250*(n) for SVHN. This deep network has 150 epochs to learn where it optimizes the weights with the ADAM optimizer [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF], using an initial learning rate of 0.002 that has an annealing phase of 50 epochs where the learning rate decays linearly to O.

Assembly Code

IV. RESULTS

Experimenting with codes

Assigning the same number of output nodes to the different classes may not be ideal for a real-world data-set. Following this idea, it may be interesting to allocate different amounts of the output material to the different combinations of classes, as for instance the distributions of examples can be more correlated in a subset of classes than on average in the whole data-set. t. ';-1-�0;;':. 2 ;------; ; 0:':; . 3----; ;' 0. 7 4 ---;O;";. 5 :------; ;" 0.';:-6-�0,.. . 7 ;------; ; 0':: . 8--0::". 9 ;;------,-1 1. 0 proportion of Is in output codes We see here that the error rate is maximal when there are less than 10% of ones in the output codes. It decays with higher values down to a minimum reached when there are around 60% of ones. Above that proportion, the error rate raises again. The minimal distance is obtained when generating as many ones as zeros on output as shown on figure 3, whereas the classifier performs better with output codes made of 60% up to 80% of ones. The error rate is also about 1. 7 times lower for 90% of ones as compared to the case with 10% of ones, while the minimal distance is the same in both cases. This asymmetry is due to the way the class label is selected at test time, where a majority voting procedure is applied in which each output unit getting a value of I increments the score of all classes it is associated with. In this setting, the higher the number of ones in the output targets, the more connection weights end up being involved in the decision process at test time. The classifier thus makes use of a finer features-output mapping in this case. With too many ones however, e. g. 90%, the distance between codes becomes too low which affects performance.

MNlST -Shallow Network

The first test with the assembly codes defined in section 3. 3, applies the shallow network to the MNIST data-set. The goal of this test is to compare the assemblies in similar settings, rather than achieving a competitive result to the state-of-the-art. The results indicate that Quintuplets are better than Grandmother cells and also that Couples are better than Grandmother cells, but it is inconclusive in the comparison between Couples and Quintuplets (less than two misclassified images of difference on average between the best networks of the two assemblies). Another test is conducted on the MNIST data-set, now trying to emulate the results from [START_REF] Pezeshki | Deconstructing the Ladder Network Architecture[END_REF]. It is summarized in figure 5. Despite respecting the hierarchy of Quintuplets 2: Couples 2: Grandmother cells, the results are too close to take any conclusions (less than one image on average between 5Q and l26G). Finally, applying the deep network to the SVHN data-set allows us to obtain more significant results than the ones obtained over MNIST. SVHN is more difficult to classifY and has been less extensively studied. The results respect the hierarchy drawn from coding theory in section 2 (Quintuplets 2: Couples 2: Grandmother Cells), with an average distance of 0. 11 % (� 26 images) between the worst quintuplets network and the best couples network and 0. 33% (� 85 images) between the worst quintuplets network and the best grandmother cells.

iii iii IG iii iii 25G f f 126G t t SC iii ! lC iii iii 6C f f 28C iji iji 10 iii iii 50
On both MNIST and SVHN. these results show the interest of the quintuplets code which gives better performance than most other output codes for a comparable number of parameters. The only case where a code outperforms the quintuplets is the "special couples" (SC) code of non-repeated couple cells clusters, which beats 5Q on MNIST. It is also worth noting that the non-repeated couples perform better than the repeated 28C on both data-sets.

V. CONCLUSIONS AND PERSPECTNES

A way to represent categories in multi-class problems is presented, that departs itself from the usual "grandmother cell" approach. Experimental results show that an assembly of neurons representing meta-classes can do a better job as output for a neural network. On the widely studied MNIST data-set, we use as starting point a multi-layer perceptron network that is close in performance to the state-of-the-art, and show that the proposed assembly codes can improve its accuracy. Furthermore, these results are in accordance with predictions drawn from coding theory in that a higher minimal Hamming distance between code words typically results in a better training of the classitier. The comparison also holds when code repetition is used to adjust the lengths of the outputs, so that all compared networks have about the same number of parameters.

However one of the goals underlying this work is to combine deep neural networks with clustered clique-based associative memories, and the quintuplet contiguration which gives the best results here would not be a good tit for this use. Indeed, even if it has a good minimal distance, the resulting clique patterns would have too many nodes in common. In this lO-class case, couple cells would be a good trade-off with a high minimal Hamming distance and a good suitability with Clustered Cliques Networks. Similar trade-off's may also be found for higher-dimensional classitication problems.

A prospect of our ongoing work is also to be able to classify data sets with a high number of classes. This way the compression power of the clustered clique code may be exploited, which allows to support a big dictionary of code words with a limited amount of material while keeping a good word recognition ability. The assembly output deep learner presented is a first step in this direction. Sets of images of words are a very interesting application case as it can have tens of thousands of categories.

5 )Figure 1 .

 51 Figure 1. The output layer of the classifier is organized in c clusters, each one having I nodes. These I nodes represent disjoint combinations of m classes among P with I = Plm. Now we propose that the n assemblies are distributed among c

3. 1 .

 1 MNJST MNIST [7] is a data-set of grey-scale images of handwritten digits that is widely used in deep learning. The data-set is made of 60000 images targeted for training and 10000 test examples. Many published works use the first 50000 examples from the training set to actually train the network, saving the last 10000 examples to perform cross validation. Some papers however present networks trained on the whole set of 60000 training images, as is the case in [8]. Here the 50000/1 0000 split of the training set is used for the experiments of section 4.1, whereas in sections 4.2 and 4.3 the networks are trained with all 60000 examples.

  3.2. SVHNSVHN[START_REF] Netzer | Reading digits in natural images with unsupervised feature learning[END_REF] is a data-set made of color images of digits captured in real-world situations, e.g. house numberings. It contains 73257 training images and 26032 test images. We use 60000 examples from the train set to actually train our networks, and the remaining 13257 serve for cross-validation.

3. 4

 4 Neural network settings For tests, the neural networks used are multi-layer perceptrons. Two architectures of network are used, one that is shallow and the other deep. The shallow network has only 1 hidden layer, while the deep network has 5 hidden layers. The results are presented by the mean and standard deviation over 10 executions with different weight initialization, noise and image order. In sections 4. 1 and 4. 2 the neural network is shallow and its only hidden layer is composed of 2000 feature units. Training lasts for 200 epochs with a constant learning rate of 0. 1. The model "baseline + noise" from [8] is chosen as a base for the deep network used in sections 4.3 and 4.4. This is a network that gets close to the current state-of-the-art in the task of permutation-invariant MNIST. It has 6 layers, where the first 5 of them are fully connected layers (with sizes 1000-500-250-250-250

An

  experimental scheme inspired from the quintuplet configuration is tested, with 252 output neurons parted in 126 clusters of size 2 where a soft-max is applied. Instead of rigorously associating an output cell with each possible quintuplet however, 10 binary output codes are generated following a binomial distribution. This distribution is modulated to make vary the average number of ones in the output codes. Used as output for a shallow network, the average number of ones has an impact on classification performance, as shown by figure2.

3. 0

 0 ,---" "T" "" "-----r--r--. .,. .. .. ..-�-. .. .. .. .,.--�-�-_ __,

Figure 2 :

 2 Figure 2: Influence of the proportion of ones in the output target vectors on classification error rates.

  130 ,---.---_,----,_--�--_,----,_--�--_.-

Figure 3 :

 3 Figure 3: Measured mean Hamming distance between generated class codes depending on the proportion of ones.

Figure 4 :

 4 Figure 4: Results summary for a shallow network on the MNIST data-set

Figure 5 :

 5 Figure 5: Results summary for a deep network on the MNlST data-set

  Number of parameters in millions

Figure 6 :

 6 Figure 6: Results summary for a deep network on the SVHN data-set

  ) . The minimum Hamming distance of the

	m-2

, ;::: g ! :H; :� H:'��r o l�'

Table 1 :

 1 Summary of the tested assembly codes. Assembly Code

* number: Code repeated number times. ': 252 non-repeated clusters

ACKOWLEDGEMENT

The authors would like to thank NVIDIA for providing us with a free graphics card allowing to speed up computations for the experiments performed during this work. The software developed to perform the simulations presented here was based on libraries Theano [12] and Keras [13].

This work was supported by the European Research Council under Grant ERC-AdG2011 290901 NEUCOD 978-1-5090-3401-7116/00 ©2016 IEEE