A HOMOGENEOUS MODEL FOR COMPRESSIBLE THREE-PHASE FLOWS INVOLVING HEAT AND MASS TRANSFER.

. A homogeneous model is proposed in order to deal with the simulation of fast transient three-phase ﬂows involving heat and mass transfer. The model accounts for the full thermodynamical disequilibrium between the three phases in terms of pressure, temperature and Gibbs enthalpy. The heat and mass transfer between the phases is modeled in agreement with the second law of thermodynamics, which ensures a stable return to the thermodynamical equilibrium. The set of partial diﬀerential equations associated with this model is based on the Euler set of equations supplemented by a complex pressure law, and by six scalar-equations that allow to account for the thermodynamical disequilibrium. It therefore inherits a simple wave structure and possesses important mathematical properties such as: hyperbolicity, unique shock deﬁnition through Rankine-Hugoniot relations, positivity of the mixture fractions. Hence the computation of approximated solutions is possible using classical algorithms, which is illustrated by an example of simulation of a steam-explosion.


Introduction
The steam explosion phenomenon occurs in industrial plants when some heated materials (solid or molten solid) come into contact with cold water [4]. The brutal heat transfer from the heated material to the liquid leads to a sudden and brutal production of steam. This steam expands quickly and strong shock waves are produced in the liquid phase, which propagate inside the devices and may damage some of them. When the steam bubble expands in an open domain, for example in a pool with a free surface, some hot materials (solid, steam and/or liquid) are expelled at high velocity. Such a phenomenon occurs in the steel industry (foundry), causing casualties and damages. In the framework of the safety demonstration for the nuclear power plants, some specific scenarii involving steam explosion are studied. This is for instance the case for the Fuel Coolant Interactions (FCI) which occur in the Reactivity Initiated Accident (RIA) [25] or when the corium comes into three-phase mixture within the whole flow. On the one hand, the second law of thermodynamics allows to define: the thermodynamical properties of the mixture (the mean pressure and the mean temperature) and the time-evolution of the setting of the three phases within a volume of mixture. Thanks to the dissipation property associated with the second law of thermodynamics, this evolution represents a stable return to the thermodynamical equilibrium when considering a closed volume of mixture. On the other hand, the time-evolution of a mixture-volume is classically described through the first law of thermodynamics and Newton's law. We then end up with a system of equations which is based on the Euler set of equations associated with a complex pressure law and supplemented by six equations that account for the thermodynamical disequilibrium between the phases.
In Section 1 the model is built on the basis of [26,31]. The closures that are then obtained ensure some good mathematical and physical properties for the system of equations: hyperbolicity of the convective part of the system, uniqueness of the definition of the shocks, entropy dissipation, stable return to thermodynamical equilibrium. These properties are gathered in Section 2. They are a serious advantage when the goal of the model is to be used to perform numerical simulations including shock patterns. The model proposed in Section 1 remains a quite general three-phase flow model, and in Section 3 we introduce some specific features that are mandatory for the steam-explosion modeling. We then specify the fact that mass transfer only occurs between the liquid water and its vapour phase, and we introduce the heating of a solid phase through a source term. At last, in Section 4, we present the numerical simulation of the heating of a solid phase mixed with liquid which leads to steam generation and strong pressure waves. In this simulation, the water phases are described using the IAPWS 97 thermodynamical laws [36].

Modeling three-phase flows with a homogeneous model
In this section, we build a homogeneous model which describes a mixture of three-phase. The assumptions are introduced throughout the section. In order to write the model we proceed in two steps by adopting a Lagrangian point of view. We first propose to model the thermodynamic behaviour of a volume of the mixture using the second law of thermodynamic. This first step of the modeling process follows the process proposed in [26,31] for two-phase flows. Then the first law of thermodynamics and Newton's law are applied to describe the evolution of this volume within the whole flow.

Some definitions and assumptions
Let us consider a volume V (in m 3 ) of the three-phase mixture which is associated with a mass M (in kg) and an internal extensive energy E (in J). Each phase i = {1, 2, 3} occupies a volume V i , has a mass M i and an internal energy E i . We assume the following properties for the mixture.
(H 1a ) The geometric repartition of the phases inside the volume V is not taken into account.
(H 1b ) The surface tension is neglected.
(H 2 ) The three phases are not miscible.
(H 3 ) Vacuum occurrence is not considered here. With these assumptions, the volume V, the mass M and the internal energy E can be written: It allows to treat naturally the cases where only one or two of the three phases are present. The hypothesis (H 2 ) and (H 3 ) are mandatory to write the first equation of (1) on the volumes. In [18,32], the miscible case has been investigated and it leads to a different system (1). Assumption (H 3 ) implies that we consider that V, M and E are non-negative: (V, M, E) ∈ (R * + ) 3 .

The second law of thermodynamics
In this section we use the second law of the thermodynamics to define the time evolution of the quantities (V i , M i , E i ) for an isolated mixture, that is for a fixed (V, M, E). We assume that the extensive phasic entropies η i (in J/K) are defined such that the following properties hold: Remark 1. When M i = 0 (i.e. when phase i exists), assumption (H 6 ) allows to define a specific entropy (in J/K/kg) from the extensive entropy η i by setting a = 1/M i . The specific entropy s i , thus, only depends on where the second equality assumes an abuse of notation with respect to the dimension of the entropy. Hence, η i is a complete Equation Of State (EOS), from which we define the pressure P i , the temperature T i and the Gibbs enthalpy µ i (or Gibbs free enthalpy, in J/kg): and It should be noticed that the assumption (H 7 ) is equivalent to ensure that the temperature T i are non-negative. Moreover, these definitions imply the classical Gibbs relation used in the Classical Irreversible Thermodynamics (CIT) theory: Remark 2. Recalling assumption (H 6 ), we have: Hence, by deriving this relation with respect to a and by applying a = 1, we get: and therefore using the definitions (2), (3) and (4) we obtain the relation The thermodynamic behaviour of the phase i is defined by the entropy η i and the Gibbs relations (5) and (6). We now assume that the extensive entropy of the mixture η is the sum of the extensive entropy of each phase: where for the sake of simplicity, we set Thanks to the assumptions (H 4 ) − (H 6 ) on the phasic entropies, the mixture entropy η satisfies the properties: 9 , η(aW ) = aη(W ). The details of the proof of these properties can be found in Appendix 5.1. By deriving the mixture entropy defined by (7) and by using the phasic Gibbs relations (5) we find: This relation can be rewritten in terms of the mixture quantities (V, M, E) by using the chain-rule Relation (9) is the Gibbs relation for the mixture, from which the mixture temperature T , the mixture pressure P and the mixture Gibbs enthalpy µ can be defined. Indeed we have: which also means that the mixture temperature T is non-negative.
Until now, we have considered the extensive mixture entropy η, which is defined on (R + ) 9 . We propose now to introduce the intensive entropy. For this purpose, let us define H(M), the subset of (R + ) 9 such that: It can be proved thatη is strictly concave on H(M). The detail of the proof is given in Appendix 5.2.
Some thermodynamical properties of the mixture have been examined above. As in [3] we choose to assume that the time-evolutions of the quantities (V i , M i , E i ) for a fixed (V, M, E) (i.e. for an isolated system) are of the form: where the quantities (V i , M i , E i ) and the time-scale λ > 0 have to be defined. The second law of thermodynamics applied to our system states that: when it is isolated, the mixture entropy must increase. In other words, when dV = dM = dE = 0, the models chosen for dV i , dM i and dE i must lead to an increase of the mixture entropyη. The quantities (V i , M i , E i ) and the time-scale λ are chosen to comply with the second law of thermodynamics. Sinceη is strictly concave, the plane which is tangent toη at any pointW of H(M) is aboveη. This can be written: By derivingη with dV = dM = dE = 0 we get: where Then, thanks to the inequality (15), we have: When the maximum W is reached in the interior of the domain D(V, M, E), i.e. when the three phases coexist, the derivative of the entropy with respect to V i , M i and E i vanish and by the Gibbs relation (19) we get that the pressure, the temperature and Gibbs enthalpy of all the phases are equal: When the maximum W is not reached in the interior of D(V, M, E), the three phases do not coexist. The maximum is then reached on a bounadry of the domain and at least one phase is not present. In such a case, the equilibrium state may be composed of two phases, say phases i and j = i, such that their pressures, temperatures and chemical potentials are equal: with V k = 0, M k = 0, E k = 0 for k / ∈ {i, j}. If there does not exist a couple of phases ensuring (21) and if (20) has no solution, then the equilibrium state corresponds to a single-phase state containing the phase i that possesses the maximum entropy η i (V, M, E). Assumption (H 10 ) is in fact mandatory and implicitly admitted in Section 1.2 to write the Gibbs relation for the mixture entropy. Assumption (H 11 ) allows to define some specific quantities in Section 2. It leads to the equation:

The first law of thermodynamics and Newton's law
Assumption (H 12 ) is classical and reads: For the sake of simplicity, we only consider here the force due to the pressure gradient. The momentum equation arises from (H 13 ) and can be written: where the pressure P is the same that the one defined in (11). At last, the first law of thermodynamics (H 14 ) states that the variation of the energy E is due to the work of the external forces and to the heat Qdt (in J) supplied to the system by its surroundings. Since we only consider here the forces due to the pressure, the variation of the energy is:

The set of PDE in intensive form
The equations (14), (22), (23), (24) and (25) define the evolution of the quantities The mass conservation (22) allows us to write the model using specific quantities (per unit of mass). Therefore we define the specific volume of the mixture τ = V/M (in m 3 /kg), the specific energy of the mixture e = E/M (in J/kg). The specific entropy S (in J/K/kg) is defined as the entropy per unit of mass: Thanks to the properties of the entropyη and the entropies η i , we have ∀W ∈ H(M): which for y i = 0 and by using the notation of Remark 1 gives: The volume fraction α i = Vi V of phase i, the mass fraction y i = Mi M of phase i and the energy fraction z i = Ei E of phase i play an important role in the model since they describe how the phases are mixed to compose (V, M, E).
In the following, they will be respectively denoted by α i , y i and z i . The set of equations (1) implies that we have: and The Gibbs relation for the specific entropy S can be deduced from the Gibbs relation (9), and it is: with the mixture temperature and the mixture pressure: The equations (14) on can also be re-written using the fractions: where the equilibrium fractions are At last, the equations of the previous subsection (23), (24) and (25) read using the specific quantities: whereQdt = Q/Mdt is the specific heat (in J/kg) supplied to the system. In the following, we setQ = 0 and a specific emphasis on the heating source term is proposed in Section 3. The derivative dφ of a variable φ corresponds here to the derivative along a streamline of the flow, which can also be written: hence, the set of equations (27), (31) and (32) can be written in conservative form: where ρ = 1/τ is the mixture density, and E = e + |U | 2 /2 is the specific total energy of the mixture. The fraction vector Y gathers the fractions of phase 1 and 2: Y = (α 1 , y 1 , z 1 , α 2 , y 2 , z 2 ). The fractions of the third phase are deduced from Y through the relations (27). The source-term vector Γ Y is then: The temperature law and the pressure law for the mixture are given by definitions (10) and (11). By using the specific quantities they read:

Properties of the whole model
In this section, we present the main mathematical properties of the model. We focus here on the properties that are mandatory for a model to be used in a numerical simulation process [16].
Without any loss of generality, and since the system is invariant under frame rotation, we consider here for the sake of simplicity system (33) for a one-dimensional space variable x, that is: with Y = (α 1 , y 1 , z 1 , α 2 , y 2 , z 2 ) and E = e + U 2 /2. The closure relation for the pressure is given by relations (10) and (11): where . The phasic pressure and temperature laws P i and T i must be specified by the user.
The sound speed c of system (36) is defined as: Using formulas (37) for the mixture pressure P and the mixture temperature T , it can be written: where d 2 s i stands for the Hessian matrix of the phasic entropies (τ i , e i ) → s i (τ i , e i ): We recall that the phasic sound speeds are defined as: It must be emphasized that this mixture celerity c is, thus, not a barycenter of the phasic celerities c i . We already know from (13) that the mixture temperature T is non-negative. Hence, if the specific phasic-entropies (τ i , e i ) → s i (τ i , e i ) are strictly concave, the square of the mixture sound-speed c 2 is non-negative. Assumption (H 5 ) on the concavity of the phasic entropies η i , implies that the specific entropies s i (see Appendix 5.3) are also concave and thus that c 2 ≥ 0. In order to enforce the strict hyperbolicity of the model, the user has to specify strictly concave EOS (τ i , e i ) → s i (τ i , e i ) for each phase.
The convective part of system (36) is based on the Euler set of equations associated with a complex mixture EOS. Hence, it inherits from the eigenstructure of the Euler system and is composed of three waves: a contact discontinuity associated with the velocity U , and two genuinely non-linear waves U ± c. Provided that c > 0, the fractions Y are constant across the shock waves since they are associated with the contact discontinuity U . It can be proved (see Appendix 5.3) that the specific entropy (τ, e) → S(Y, τ, e) is strictly concave, which ensures that the shock waves associated with the genuinely non-linear waves U ± c are uniquely defined through the Rankine-Hugoniot relations.
Concerning the equilibrium fraction Y , it has been shown in Appendix 5.3 that the specific entropy S(Y, τ, e) is strictly concave with respect to Y on [0, 1] 3 . Therefore, there exists a unique equilibrium fraction Y (τ, e) corresponding to the maximum of S at a given (τ, e). We now recall the following classical lemma.

Lemma.
Let Ω be a subset of R. Let Ψ, a, Π and U be some sufficiently regular applications, with the following properties: and such that: Suppose that for all x b ∈ ∂Ω, the boundary of Ω, This lemma can be used to prove the following property.
Property. For regular solutions, and under the assumptions of the Lemma, the fraction α i , y i and z i remain in [0, 1].
In order to prove the property above, we proceed for the volume fractions α i , and the same demonstration holds for y i and z i . Since i α i = i α i = 1, we have from (36): with a i = −1/λ and Π i = α i /λ. The equilibrium volume fraction α i belongs to [0, 1] and the time scale λ must be chosen non-negative, so that we obviously have Π i ≥ 0. The lemma can be straightforward applied for Ψ = α i which remains positive under the assumptions of the lemma. This proof can obviously be extended to the fractions y i and z i which, thus, also remain positive.

The specific case of the steam explosion
The model proposed in the sections above deals with general three-phase flows. It can for example be used to perform simulations involving the same material in liquid, vapor and solid state (for the simulation of the sudden depressurization of CO 2 pipes for instance). Relations (20) then define the triple point when the three phases coexist in a stable manner. When only two phases coexist in a stable manner, say phase 1 and phase 2, system (20) reduces to: which defines the saturation curve between phase 1 and phase 2 (of course the saturation curves between phase 1 and phase 3, and between phase 2 and phase 3 are defined by permuting the indices in (43)). In this section we propose further assumptions that allow to deal with the steam explosion [4]. We are interested here in situations where one of the three phases represents a material, and the two others represent the same fluid in the liquid state and in the vapour state. In the following, phase 1 stands for the liquid phase, phase 2 for the vapour phase and phase 3 for the inert phase.
Since the phase 3 is inert, its partial mass M 3 /M will remain constant: Moreover, due to the mass transfer between the liquid phase and the vapour phase, each variation of the liquid mass is balanced by the variation of the vapour phase, and conversely, which can be written: These two additional rules modify the entropy equation (19) into the following equation: Hence the Gibbs enthalpy of the inert phase µ 3 does not play any role in the entropy dissipation. Nevertheless, its pressure and its temperature are still part of the relaxation process, since the equilibrium state is now defined as: where M 3 = M 3 is constant along the streamlines. Therefore, the system of equations (36) and the closure laws (37) are not modified. Nevertheless, the equilibrium fractions must be computed using the intensive form of the relations (45), that is: where y 3 (x, t) = y 3 (x, t) for every point (x, t).
In order to deal with steam explosions, another mandatory ingredient for the model is to cope with external heating source terms. In particular, some complex chemical reactions occurring in the bulk of the inert phase may lead to an increase of its internal energy. For the sake of simplicity, these terms have been omitted in the previous sections and they are introduced here. In the following, we focus on the time variation of the mixture, and without loss of generality, we set U = 0. Let us assume that the heating of the mixture is such that: (i) the specific volume of each phase is constant, d(τ i ) = 0; (ii) the partial mass of each phase is constant, d(α i ρ i ) = 0; (iii) the internal energy of each phase is such that, d(α i ρ i e i ) = α i ρ i q i dt; where q i dt is the specific heat received by phase i. Thanks to the assumptions (i) and (ii), the specific volume of the mixture is constant d(τ ) = 0, and the volume and the mass fractions of the three phases are constant: d(α i ) = 0 and d(y i ) = 0. We also get from (iii) that the variation of the specific energy of the mixture ρe = ( i α i ρ i e i ) is: whereQ has been introduced in Section 1.4. Then the time evolution of the energy fraction z i = (α i ρ i e i )/(ρe) is: Hence, when the phases are heated by an external source q i , system of equations (36) becomes: The source terms Γ αi , Γ yi and Γ zi are the source terms associated with the thermodynamical relaxation Γ Y defined in (34). In our particular case, we recall that Γ y3 = 0 since the phase 3 is inert.
Remark 3. The positivity results of Section 2 for z i still hold provided that q i /e i andQ/e remain bounded. Indeed, the equation for z i of system (48) can also be written: with a i = −1/λ + (q i /e i −Q/e) and Π i = z i /λ. The lemma of Section 2 can then still be applied here.
Accounting for the phasic heating source terms represents an important feature of the model. The latter possesses three energy equations through: the energy fraction equations and the mixture energy equation, and this allows the user to specify how the energies of the three different phases vary according to external sources. This is typically not the case for the classical homogeneous model [9], which is widely used for industrial simulations.

An example of numerical simulation: steam explosion during a Reactivity Induced Accident
As an illustration of the capability of the model (48), we propose in this section a simple heating test case which may be seen as a simplified situation of RIA. A sketch of the test case has been plotted in Figure 1. We consider that some fuel particles are released from the fuel rod which has a radius R 1 = 5 10 −3 m. Within the ring delimited by R 1 and R 2 = 6 10 −3 m, the liquid water contains fuel particles with a volume fraction α f = 0.01, thus the liquid volume-fraction is α l = 0.99 and the steam volume-fraction α v = 0. The computational domain is the ring [R 1 , R 3 = 2 10 −2 m]. At the beginning of the simulation, t = 0, we assume that there is no vapor in the domain, and that the liquid and the fuel particles are at the same pressure P = 155.0 bars and at the same temperature T = 613 K. We recall that the saturation temperature for the water at 155.0 bars is equal to 618 K. We also assume that the initial velocity is equal to zero.
A constant heating source term q f is then applied at t = 0 to the fuel-particle phase according to system (48). This source term arises from the chemical reactions that occur in the fuel particles. It has been arbitrarily chosen equal to q f = 2 10 10 W/kg. It should be mentioned that in a more realistic RIA situation, the heating source term q f is negligible with respect to the heat transfer due to the initial temperature disequilibrium between the liquid water and the fuel particles. Nevertheless, accounting for such disequilibrium requires a realistic time scale λ. This is the case when considering pressure disequilibrium for steam-liquid configurations as reported in [23], and this is also the case for temperature disequilibrium in the case of three-phase flows. For the sake of simplicity, we consider here an instantaneous relaxation time scale λ = 0 which unfortunately does not allow for strong initial temperature-disequilibrium.
The EOS for the liquid water and for the steam are based on the IAPWS 97 formulation [36], whereas the fuel particles are modeled using a Stiffened Gas EOS [28]. The specific entropy of the fuel-particle phase then reads: where C V,f , Q f , Π f , γ f , and s 0,f are parameters. The corresponding pressure and temperature laws are: Since there is no mass transfer involving the phase f , the parameter s 0,f is useless and it has thus been set here to zero: s 0,f = 0. The other coefficients have been estimated using the data for uranium dioxide [1] at a temperature of 623 K. The specific enthalpy of the uranium dioxide, h f = e f + P f τ f , is given as its difference to the reference specific enthalpy at T = 298 K. Unfortunately, it seems that no information can be found on the latter. Hence we have arbitrarily chosen the value h 0 = 5.0 10 4 J/kg. The other data, collected in [1], are: ρ f (623 K) = 10850 kg/m 3 , C P,f (623 K) = 294 J/kg/K.
No information is available for the sound speed c f and we, therefore, impose a value of c f = 6000 m/s which is an order of magnitude of the sound speed for steel. The EOS parameters are then computed using the relations: where the pressure is P = 155.0 10 5 P a.
By the way of a comparison, we also perform a simulation without fuel particle (i.e. α f = y f = z f = 0). The heating source term q f is then applied to the vapour and to the liquid: q l = q v = q f in the fixed domain [R 1 , R 2 ]. According to (48), the specific heat received by the mixture is then (y l q l + y v q v )dt = q f dt. In the simulation with fuel particles, the mass of fuel is conserved and the total heat received by the mixture is linearly increasing with time. On the other hand, in the simulation without fuel particle, the heat is received by a fixed volume that contains a non-constant mass of mixture. Roughly speaking, the mass of fluid in the domain [R 1 , R 2 ] decreases when the time increases because the temperature increases due to the heating of the fluid. It is thus not easy to predict the total amount of heat received by the mixture for the simulation without fuel particle. For the comparison, we impose equivalent initial heating source terms. For the simulation involving fuel particles, it is equal to: and for the simulation without fuel particle: where V is the volume between R 1 and R 2 . If we impose the two quantities to be equal, this yields q f ∼ 680/108.50 q f . We insist that the total amounts of heat injected in the domains during the whole simulations are different. Since for the simulation without fuel particle, the specific heat is injected in a fixed domain with a diminishing amount of mass (due to the increase of the temperature, the density decreases), the total amount of heat injected is indeed lower than for the simulation with the fuel particles.
The numerical scheme used to obtain the approximated solutions of system (48) is classical. Since the test case considered here is symmetric with respect to the axis r = 0, system (48) is written in axi-symmetric formulation. The overall scheme is based on a fractional step approach [37] in which the convective part and the source terms are treated successively. For the convective part, the numerical scheme is a finite volume scheme [11] where the numerical fluxes are approximated by a Rusanov scheme [34]. The relaxation source terms are then solved following the scheme described in [23]. The main difficulty concerning these source terms is to compute the equilibrium fractions. This computation follows the idea of [13,23] when the equilibrium is a liquid-steam-fuel equilibrium or a liquid-steam equilibrium. When the equilibrium involves the solid phase and only one phase among the water phases, the pressure-temperature equilibrium is solved using a classical Broyden method with Sherman-Morrison update of the inverse of the Jacobian [5,6,35]. Accounting for the heating source term is the last step of the algorithm; its is discretized using a semi-implicit Euler scheme that preserves the positivity of the energy fractions. In the latter the heating source terms q i are explicited and the system of ordinary differential equations associated to the heating source term is then solved analytically. In our simulation, the source terms q i are constant so that the scheme corresponds to an exact integration of the heating source terms for each time-step.
The overall scheme described above is then used to perform the two simulations for a mesh with 4000 cells and with a uniform radial mesh-size. In Figure 2 the values across the time are plotted at r = (R 1 − R 2 )/2, whereas different pressure and volume-fraction profiles along the radius r are plotted at five different instants in Figure  3. The Figure 4 shows the trajectory in time of each simulation in the pressure-temperature plane. The results obtained for the two configurations exhibit significant differences. First of all, since the time-scale λ has been chosen equal to zero, the thermodynamical equilibrium is achieved instantaneously. Hence, when the mixture fractions are in ]0, 1[ the two phases have the same pressure and temperature which are on the saturation curve. This can be verified with the trajectories of Figure 4 and with the temperature/volume fraction curves of Figure  2. When considering the time evolution of the results in the heating zone (see Figures 2 and 4), the simulation can be split into three periods which are: • the heating of the liquid, which corresponds to the beginning of the simulation, until the steam volume-fraction becomes positive; • the vaporization of the liquid, which corresponds to the time interval for which the volume-fraction is in ]0, 1[; • and the heating of the steam, which corresponds to the end of the simulation, when the steam volumefraction is equal to 1. In the heating zone, the maximum of the pressure magnitude arises during the second period. It can be observed that the simulation with fuel particles reaches a lower pressure which then decreases slower in time. The pressure wave that is then generated through the domain (see Figure 3) has a lower magnitude. It can also be seen from Figure 2 that in the heating zone the apparition of steam arises at almost the same time (at time 1.86 10 −5 s with fuel particles and 1.94 10 −5 s without fuel particle), whereas the complete vaporization of the liquid is achieved earlier without fuel particle (at time 4.227 10 −4 s with fuel particles and 3.25 10 −4 s without fuel particle). Once the heating zone only contains steam, the temperature rapidly increases. Hence, despite the lower amount of energy received by the fluid in the simulation without fuel particle, the associated results seem to correspond to a slightly more important level of severity of the steam explosion.    . Trajectories in the pressure-temperature plane along the simulation time at the radius r = (R 1 + R 2 )/2 (i.e. the center of the initial heating zone). The simulation without fuel particle is plotted using red circles, and the simulation with fuel particles is plotted using blue squares. The black plain line represents the saturation curve.

Conclusion
The homogeneous three-phase flow model proposed here relies on Newton's law, the first and second laws of thermodynamics. It accounts for the compressibility of the three phases and for the heat and the mass transfer between the phases through the thermodynamical disequilibrium between the phases (in terms of pressure, temperature and Gibbs enthalpy). Since no equilibrium assumption is done, the model should not be restricted to the representation of thermodynamical phenomenon at "large time-scale" and fine physical aspects of the thermodynamics might be caught. This feature seems mandatory when dealing with fast transient flows induced by the flashing of some liquid. From a mathematical point of view, the resulting model possesses very interesting properties, which allows to build efficient numerical schemes. In Section 4, an example of the simulations of the flashing of liquid water undergoing a violent heating has therefore been proposed to illustrate the capability of the model to be used in an industrial configuration. Even if the relaxation time-scale used in these simulations has been set to zero, the homogeneous model proposed in this paper allows to perform transient simulations involving strong pressure waves. Actually, some realistic models for the relaxation time-scale should be proposed on the basis of physical considerations in order to improve the simulations of Section 4.

Properties of the extensive mixture-entropy
In this section, we give the proofs that the intensive mixture-entropy (7) is such that: , with a small abuse of notation, The mixture entropy then reads: η(W ) = i η i (W i ).
Proof of property (i). The use of assumption (H 4 ) gives a straightforward proof of property (i).
Proof of property (iii). We use here assumption (H 6 ) which states that: Let us choose W ∈ (R + ) 9 and a > 0. Then we have:
Obviously, since η fulfills properties (i)-(iii) recalled in Appendix 5.1, the entropyη satisfies: (i) W →η(W ) is C 2 on H(M); (ii) W →η(W ) is concave; (iii) ∀a ∈ R * + , ∀W ∈ H(M),η(aW ) = aη(W ); The demonstration given here can be found in a more general form in [26]. The sketch of the proof is the following. We first exhibit the degeneracy manifold of the Hessian of the entropy η. Then we prove that its intersection with H(M) contains a single point. As a consequence, the degeneracy manifold of the Hessian of η also resumes to a single point, which proves that it is strictly concave.
If we derive the relation above with respect to a, we find that: and, thus: ∇ 2 Y (η(Y )) · X = 0, which for a = 1 leads to: ∇ 2 Y (η(Y )) · Y = 0. The degeneracy manifold of the Hessian of the entropy η at a point W is thus the set D η (W ) = {bW, b ∈ R}. This implies that the entropy η can not be strictly concave. This proves that the kernel of the Hessianη, the restriction of η to H(M), is restricted to a single point. As a consequence,η is strictly concave on H(M).

Concavity of the intensive entropies
We first demonstrate that the phasic intensive entropy s i , introduced in Section 1.2: is concave with respect to (V i /M i , E i /M i ). We choose two sets of variables We then have for any a ∈ [0, 1]: Thanks to the concavity of (V i , M i , E i ) → η i (V i , M i , E i ) (assumption (H 5 )), we get that: Hence, following the definition of the entropy s i , we have: The assumption that M # i = M * i = M i allows to obtain: and we obviously have: This property allows to define in a unique manner the shock waves for system (36). Indeed, when considering a Riemann problem for the convective part of system (36), it can be noticed that the fractions Y only vary through the contact wave U . Hence the fractions Y are constant across a shock wave. The property above then ensures that the specific entropy S is strictly concave through the shock waves.