Manuel Bernard
email: manuel.bernard@ifpen.fr

Stéphane Dellacherie
email: stephane.dellacherie@cea.fr

Gloria Faccanoni
email: faccanon@univ-tln.fr

Bérénice Grec
email: berenice.grec@parisdescartes.fr

Yohan Penel
email: yohan.penel@cerema.fr

Study of a low Mach nuclear core

Keywords:

Introduction

Several physical phenomena have to be taken into account when modelling a water nuclear reactor such as PWRs 1 or BWRs 2 (see [START_REF] Delhaye | Thermohydraulique des réacteurs[END_REF] for an introduction). In particular, the present work deals with the handling of high thermal dilation of the coolant uid induced by thermal transfers in nuclear cores (see Figure 1 for schematic pictures of PWR and BWR reactors). A natural approach is to represent the evolution of the ow by means of a system of PDEs similar to the 1 compressible Navier-Stokes equations coupled to the modelling of phase transition as it is the case in classic industrial codes [START_REF]TRACE V5.0 Theory Manual, Field Equations, Solution Methods and Physical Models[END_REF][START_REF] Bestion | The physical closure laws in the CATHARE code[END_REF][START_REF] Fillion | FLICA-OVAP: A new platform for core thermalhydraulic studies[END_REF].

In nominal and incidental situations as well as in some accidental situations studied in safety evaluations, the magnitude of the sound velocity is much higher than the one of the velocity of the coolant uid, which means that the Mach number of the ow is small. The discretization of compressible Navier-Stokes type systems at low Mach number may induce numerical issues directly related to the existence of fast acoustic waves (see for example [START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number[END_REF][START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit[END_REF][START_REF] Volpe | Performance of compressible ow codes at low Mach numbers[END_REF] when the convective part of the compressible Navier-Stokes system is discretized by means of a Godunov type scheme). Sereval numerical techniques based on the resolution of a Poisson equation have been proposed in the literature to extend incompressible methods to the (compressible) low Mach case. A pioneering work was that of Casulli & Greenspan [START_REF] Casulli | Pressure method for the numerical solution of transient, compressible uid ows[END_REF] where a nite-dierence scheme over a staggered grid is designed by impliciting terms involving the speed of sound in order to avoid restrictive stability conditions. In [START_REF] Colella | A projection method for low speed ows[END_REF], Colella & Pao made use of the Hodge decomposition to single out the incompressible part of the velocity eld.

Nevertheless, in a low Mach number regime, the acoustic phenomena can be neglected in energy balances although the ow is highly compressible because of the thermal dilation. Thus, to overcome the numerical diculties, S. Dellacherie proposed in [START_REF] Dellacherie | On a low Mach nuclear core model[END_REF] another model obtained by ltering out the acoustic waves in the compressible model. Let us underline that this approach was rst applied to model low Mach combustion phenomena [START_REF] Majda | Simplied equations for low Mach number combustion with strong heat release[END_REF][START_REF] Majda | The derivation and numerical solution of the equations for zero Mach number combustion[END_REF][START_REF] Sivashinsky | Hydrodynamic theory of ame propagation in an enclosed volume[END_REF], then astrophysical issues [START_REF] Almgren | Low Mach number modeling of type Ia supernovae. I. hydrodynamics[END_REF][START_REF] Almgren | Low Mach number modeling of type Ia supernovae. II. energy evolution[END_REF] and thermal dilation of the interface of bubbles at low Mach number [START_REF] Dellacherie | On a diphasic low Mach number system[END_REF]. This specic kind of models have been studied from a theoretical point of view. We mention for instance [START_REF] Embid | Well-posedness of the nonlinear equations for zero Mach number combustion[END_REF] for the well-posedness of a low Mach number system and [START_REF] Gonzalez-Santalo | An exact solution for ow transients in two-phase systems by the method of characteristics[END_REF][START_REF] Majda | Simplied equations for low Mach number combustion with strong heat release[END_REF][START_REF] Penel | Existence of global solutions to the 1D abstract bubble vibration model[END_REF][START_REF] Sivashinsky | Hydrodynamic theory of ame propagation in an enclosed volume[END_REF] for the derivation of monodimensional explicit solutions.

From a numerical point of view, 2D simulations have been performed in [START_REF] Almgren | Low Mach number modeling of type Ia supernovae. I. hydrodynamics[END_REF][START_REF] Almgren | Low Mach number modeling of type Ia supernovae. II. energy evolution[END_REF][START_REF] Lai | A projection method for combustion in the zero Mach number limit[END_REF] for astrophysical and combustion applications while a numerical study of the model established in [START_REF] Dellacherie | On a diphasic low Mach number system[END_REF] was proposed in [START_REF] Dellacherie | Numerical resolution of a potential diphasic low Mach number system[END_REF]. As for the low Mach number model derived in [START_REF] Dellacherie | On a low Mach nuclear core model[END_REF] and called the Low Mach Nuclear Core (Lmnc) model, it was discretized in [START_REF] Bernard | Study of low Mach nuclear core model for single-phase ow[END_REF] in the monodimensional (1D) case. Moreover, 1D unsteady analytical solutions were also given in [START_REF] Bernard | Study of low Mach nuclear core model for single-phase ow[END_REF] which allowed to assess the numerical schemes. Notice that 2D numerical results will be presented in [START_REF] Dellacherie | 2D numerical simulation of a low Mach nuclear core model with stiened gas using Freefem++[END_REF][START_REF] Dellacherie | Study of low Mach nuclear core model for two-phase ows with phase transition II: tabulated EOS[END_REF]. Despite relevant numerical results, the approach proposed in [START_REF] Bernard | Study of low Mach nuclear core model for single-phase ow[END_REF][START_REF] Dellacherie | On a low Mach nuclear core model[END_REF] was not satisfying since it was restricted to monophasic ows. Thus, we extend in this study the results stated in [START_REF] Bernard | Study of low Mach nuclear core model for single-phase ow[END_REF][START_REF] Dellacherie | On a low Mach nuclear core model[END_REF] by taking into account phase transition in the Lmnc model. If we neglect viscous eects, the Lmnc model proposed in this paper may be seen as the low Mach number limit of the Homogeneous Equilibrium Model (HEM) [START_REF] Clerc | Numerical Simulation of the Homogeneous Equilibrium Model for Two-Phase Flows[END_REF][START_REF] Goncalvès | Numerical study of cavitating ows with thermodynamic eect[END_REF][START_REF] Jaouen | Étude mathématique et numérique de stabilité pour des modeles hydrodynamiques avec transition de phase[END_REF][START_REF] Müller | The Riemann problem for the Euler equations with nonconvex and nonsmooth equation of state: construction of wave curves[END_REF][START_REF] Voss | Exact Riemann solution for the Euler equations with nonconvex and nonsmooth equation of state[END_REF] with source terms. Let us recall that the HEM model is the compressible Euler system in which the two phases are supposed to be at local kinematic, mechanic and thermodynamic equilibria.

A crucial step in the process is the modelling of uid properties through the equation of state (EOS). It is important from a physical point of view to match experimental data and from a mathematical point of view to close the system of PDEs. In the present work, this point is achieved by using the stiened gas EOS. A major result in this paper is the exhibition of 1D unsteady analytical solutions with phase transition (see Proposition 3.4). These solutions are of great importance: on the one hand they provide accurate estimates of heat transfers in a nuclear core in incidental and accidental situations, and on the other hand they enable to assess the robustness of the monodimensional numerical schemes presented in this article. In addition, regardless of the EOS used in the pure phases, when the thermodynamic pressure is constant and BWR is the steam void formation in the core of the latter. mechanic and thermodynamic equilibria, the mixture can always be considered a stiened gas (this point will be claried in the sequel): this important remark legitimates the study of the Lmnc model together with the stiened gas EOS.

Compared to the numerous low Mach number combustion models derived in the literature, the Lmnc model diers for several reasons. Some of them are due to the underlying uids: indeed, combustion issues are related to gas modelled by the ideal gas law and involved in miscible mixtures whereas our modelling comprises a more general equation of state and mixtures of immiscible uids. We must also mention that the system of PDEs is set in a bounded domain with nonperiodic boundary conditions whose inuence upon theoretical and numerical investigations is noticeable.

At last, we wish to underline that although this study is specic to dimension 1 (which is essential to obtain in particular the unsteady analytical solutions with phase change), it remains useful from an industrial point of view since many safety evaluations use a 1D modelling to describe the ow in each component of the nuclear reactor and thus within the nuclear core [START_REF] Bestion | The physical closure laws in the CATHARE code[END_REF][START_REF] Drouin | A consistent methodology for the derivation and calibration of a macroscopic turbulence model for ows in porous media[END_REF].

Nevertheless, the extension of this work to dimensions 2 and 3 is a natural and important perspective [START_REF] Dellacherie | 2D numerical simulation of a low Mach nuclear core model with stiened gas using Freefem++[END_REF]. This paper is organized as follows. In Section 1, the Lmnc model is recalled together with boundary/initial conditions and assumptions under which it is valid. We also study the existence of (more or less) equivalent formulations of the model that can be used depending on the variables we aim at focusing on. Section 2 is devoted to the modelling of phase transition and to the EOS that is necessary to close the system. In Section 3, we prove some theoretical results stated (without proof ) in [START_REF] Bernard | Study of low Mach nuclear core model for single-phase ow[END_REF] and we extend them to the multiphasic case. Exact and asymptotic solutions are thus exhibited. Numerical aspects are then investigated in Section 4. More precisely, we present some numerical schemes based on the method of characteristics proposed in [START_REF] Penel | An explicit stable numerical scheme for the 1D transport equation[END_REF]. We then prove that these schemes preserve the positivity of the density and of the temperature for any time step. Finally, these schemes are applied in Section 5 to various situations with occurrence of phase transition like a simplied scenario for a Loss of Flow Accident (LOFA) induced by a coolant pump trip event.

The low Mach nuclear core model

The Low Mach Nuclear Core (Lmnc) model introduced in [START_REF] Dellacherie | On a low Mach nuclear core model[END_REF] is obtained by ltering out the acoustic waves in a compressible Navier-Stokes type system. This is achieved through an asymptotic expansion with respect to the Mach number assumed to be very small in this framework.

One of the major consequences is the modication of the nature of the equations: the ltering out of the acoustic waves which are solutions of a hyperbolic equation in the compressible system introduces a new unknown (namely the dynamic pressure) which is solution to an elliptic equation in the Lmnc model. Another consequence is that we are able to compute explicit monodimensional unsteady solutions of the Lmnc model with or without phase transition 3 (see Section 3) and to construct 1D robust and accurate numerical schemes 4 (see Section 4).

In this section, we recall the Lmnc model and we present equivalent formulations for smooth solutions. Since we are interested in the 1D case in this paper, we do not extend results to 2D/3D. Nevertheless, this can easily be done (provided the boundary conditions are adapted).

Governing equations

The 1D nonconservative formulation of the Lmnc model [START_REF] Dellacherie | On a low Mach nuclear core model[END_REF] is written as

         ∂ y v = β(h, p 0 ) p 0 Φ(t, y), (1.1a)
ρ(∂ t h + v∂ y h) = Φ(t, y), (t, y) ∈ R + × [0, L] (1.1b) ∂ t (ρv) + ∂ y (ρv 2 + p) = F(v) -ρg, (1.1c)
where v and h denote respectively the velocity and the (internal) enthalpy of the uid. Pressure p 0 is a given constant see below. The density ρ = ρ(h, p 0 ) is related to the enthalpy by an equation of state (EOS) see Section 2. So does the dimensionless compressibility coecient β(h, p 0 ) which is dened by

β(h, p 0 ) def = - p 0 ρ 2 (h, p 0 ) • ∂ρ ∂h (h, p 0 ).
(1.

2)

The power density Φ(t, y) is a given function of time and space modelling the heating of the coolant uid due to the ssion reactions in the nuclear core. Finally, g is the gravity eld and F(v) models viscous eects: the classic internal friction in the uid, and also the friction on the uid due to technological devices in the nuclear core (e.g. the friction on the uid due to the fuel rods). In the sequel, we take

F(v) = ∂ y (µ∂ y v).
In this case, µ is a turbulent viscosity given by an homogenized turbulent model. Nevertheless, we explain in the sequel that the exact choice of F(v) is not important in the 1D case (this is no more the case in 2D/3D).

We must also emphasize that model (1.1) is characterized by two pressure elds, which is classic in low Mach number approaches: p 0 : this eld is involved in the equation of state and is named the thermodynamic pressure.

It is assumed to be constant throughout the present study and corresponds to an average pressure in the nuclear core; p(t, y): it appears in the momentum equation (1.1c) and is referred to as the dynamic pressure. It is similar to the pressure eld which appears in the incompressible Navier-Stokes system.

Notice that p 0 + p(t, y) is a 1st-order approximation of the classic compressible pressure in the nuclear core: the pressure p(t, y) is indeed a perturbation around the average pressure p 0 (this is due to the low Mach number hypothesis [START_REF] Dellacherie | On a low Mach nuclear core model[END_REF]). In the 1D case, equation (1.1c) decouples from the two other equations and may be considered a post-processing leading to the computation of p (this is why the expression of F(v) is not essential in 1D). Thus, equation (1.1c) will often be left apart in the sequel and equations (1.1a)-(1.1b) will often be referred to as the Lmnc model for the sake of simplicity.

Supplements

From now on, we suppose that:

Hypothesis 1.1. 1. Φ(t, y) is nonnegative for all (t, y) ∈ R + × [0, L].
2. p 0 is a positive constant.

The rst assumption characterizes the fact that we study a nuclear core where the coolant uid is heated. In the steam generator of a PWR type reactor (see Figure 1a) which could also be modelled by means of a Lmnc type model the uid of the primary circuit heats the uid of the secondary circuit by exchanging heat through a tube bundle: in that case, we would have Φ(t, y) ≤ 0 in the primary circuit and Φ(t, y) ≥ 0 in the secondary circuit.

The second assumption corresponds to real physical conditions even if it is not required in the setting of the equation of state (see Section 2.3).

Boundary conditions (BC) The uid is injected at the bottom of the core at a given enthalpy h e and at a given ow rate D e . We also impose the dynamic pressure p at the exit of the core (y = L). The inlet BC are h(t, 0) = h e (t),

(1.3a) (ρv)(t, 0) = D e (t), (1.3b)
while the outlet BC is p(t, L) = 0.

(1.3c)

The entrance velocity v e (t) to apply at y = 0 is deduced from the relation

v e (t) def = D e (t) ρ e (t)
where ρ e def = ρ(h e , p = p 0 ).

(1.4)

The fact that h e and D e depend on time enables to model transient regimes induced by accidental situations. For example, when D e (t) tends to zero, it models a main coolant pump trip event which is a Loss Of Flow Accident (LOFA) as at the beginning of the Fukushima accident in the reactors 1, 2 and 3.

We also assume in the sequel that: Hypothesis 1.2.

1. D e is nonnegative.

2. h e is such that ρ e is well-dened and positive.

The rst assumption corresponds to a nuclear power plant of PWR or BWR type: the ow is upward 5 . The second assumption means that the EOS ρ(h, p) is such that ρ(h e , p 0 ) can be computed. Moreover, we also suppose that ρ(h e , p 0 ) > 0 from a physical point of view.

We nally make the following modelling hypothesis:

Hypothesis 1.3. β is nonnegative.
Positivity assumptions about Φ, D e , β and ρ in Hypotheses 1.1, 1.2 and 1.3 ensure that the velocity v(t, y) remains nonnegative at any time and anywhere in the core. Otherwise, the system could become ill-posed (see Section 4.2 in [START_REF] Dellacherie | On a low Mach nuclear core model[END_REF] where this question is partially addressed).

Well-prepared initial conditions The model is nally closed by means of the initial condition h 0 (y) = h(0, y) satisfying the following hypothesis:

Hypothesis 1.4.

1. We impose the compatibility condition h 0 (y = 0) = h e (t = 0). 2. h 0 is such that ρ 0 is well-dened and positive.

The initial density deduced from the EOS directly satises the equality ρ 0 (y = 0) = ρ e (t = 0). Secondly, as system (1.1) consists of steady and unsteady equations, the initial velocity v 0 must satisfy equation (1.1a) for t = 0, which means v 0 (y) = β h 0 (y), p 0 p 0 Φ(0, y).

Hence, h 0 prescribes the initial velocity v 0 through the previous dierential equation together with the condition v 0 (y = 0) = v e (t = 0). The initial ow rate D 0 is thus given by D 0 (y) = ρ 0 (y)v 0 (y).

Such initial data h 0 and D 0 are said to be well-prepared. This will be implicitly assumed in the sequel. 5 The ow could be downward when the nuclear reactor is a material testing reactor.

Origin and dierent formulations of the model

The 1D Lmnc model (1.1) is written in [START_REF] Dellacherie | On a low Mach nuclear core model[END_REF] as

   ∂ y v = β(h, p 0 ) p 0 Φ(t, y), (1.5a) ρ(h, p 0 ) • (∂ t h + v∂ y h) = Φ(t, y).
(1.5b)

We recall that in the 1D case, equation (1.1c) is a post-processing of (1.5). It is important to note that the low Mach number model (1.5) is justied only under smoothness assumptions. To study the existence of weak solutions, it might be better to use a conservative formulation which is equivalent to (1.5) for smooth solutions. This conservative formulation is the following: Proposition 1.1. Under smoothness assumptions, system (1.5) is equivalent to

∂ t ρ + ∂ y (ρv) = 0, (1.6a)
∂ t (ρh) + ∂ y (ρhv) = Φ(t, y).

(1.6b) System (1.6) (coupled to equation (1.1c)) is the Lmnc model written in conservative variables.

Although (1.6) is more general than (1.5), system (1.5) is interesting as it emphasizes the fact that the ltering out of the acoustic waves turns the hyperbolic nature of the compressible Navier-Stokes system (related to the acoustic waves) to an elliptic constraint (upon the velocity) similar to the incompressible case.

Moreover, under smoothness assumptions and for a particular class of EOS, we can derive a semi-conservative formulation equivalent to (1.5), and which may be useful to derive ecient numerical schemes. Indeed, we have the following proposition.

Proposition 1.2. Under smoothness assumptions:

1. System (1.5) implies    ∂ y v = β(h, p 0 ) p 0 Φ(t, y), (1.7a)
∂ t ρ(h, p 0 )h + ∂ y ρ(h, p 0 )hv = Φ(t, y).

(1.7b)

2. For equations of state such that

∂ρ ∂h (h, p 0 ) = - ρ(h, p 0 ) h , (1.8)
systems (1.5) and (1.7) are equivalent.

Condition (1.8) upon ρ seems to be quite restrictive insofar as it does not enable to handle ideal gas (for which ∂ρ ∂h = -ρ h ). In the latter case, equations (1.7a) and (1.7b) are nothing but the same equation, which implies that we have to use formulations (1.5) or (1.6).

Remark 1.1. The equivalence between systems (1.5), (1.6) and (1.7) also holds in higher dimensions. Nevertheless, the momentum equation is strongly coupled to the other equations in 2D/3D and must be taken into account under conservative or nonconservative forms. Indeed, these forms are equivalent (as soon as the unknowns are smooth).

Proof of Proposition 1.1.

• ⇒ According to denition (1.2) of β, we have

∂ t ρ + ∂ y (ρv) = ∂ρ ∂h (∂ t h + v∂ y h) (1.5b) = Φ ρ +ρ ∂ y v (1.5a) = βΦ p 0 = 0 (1.9)
which gives (1.6a). We also obtain

∂ t (ρh) + ∂ y (ρhv) = h ∂ t ρ + ∂ y (ρv) + ρ (∂ t h + v∂ y h) = Φ,
using (1.5b) and (1.6a). We recover (1.6b).

• ⇐ Because of (1.6a), we deduce (1.5b) from (1.6b). Moreover, we deduce from (1.6a) and

(1.5b) that

∂ t ρ + ∂ y (ρv) = ∂ρ ∂h (∂ t h + v∂ y h) + ρ∂ y v = ∂ρ ∂h • Φ ρ + ρ∂ y v = 0,
which gives (1.5a) thanks to denition (1.2) of β.

Proof of Proposition 1.2. The rst point is a direct consequence of Proposition 1.1 since

∂ t ρ(h, p 0 )h + ∂ y ρ(h, p 0 )hv = ρ(h, p 0 ) • (∂ t h + v∂ y h) (1.5b) = Φ +h ∂ t ρ(h, p 0 ) + ∂ y ρ(h, p 0 )v
(1.9) = 0 .

To prove the second point, we just have to show that (1.7) implies (1.5) under condition (1.8).

On the one hand, since

∂ t ρ + ∂ y (ρv) = ∂ρ ∂h (h, p 0 )(∂ t h + v∂ y h) + ρ(h, p 0 )∂ y v,
by using (1.2) and (1.7a), we obtain

∂ t ρ + ∂ y (ρv) = ∂ρ ∂h (h, p 0 ) ∂ t h + v∂ y h - Φ ρ(h, p 0 )
.

(1.10)

On the other hand, (1.7b) leads to

ρ(h, p 0 ) • (∂ t h + v∂ y h) + h ∂ t ρ(h, p 0 ) + ∂ y ρ(h, p 0 )v = Φ(t, y)
that is to say

∂ t ρ + ∂ y (ρv) = - ρ(h, p 0 ) h ∂ t h + v∂ y h - Φ ρ(h, p 0 )
.

(1.11) Thus, by comparing (1.10) and (1.11), we obtain

∂ t h + v∂ y h - Φ ρ(h, p 0 ) = 0
under condition (1.8), which proves that (1.7) implies (1.5).

Equation of state for two-phase uids

For the system to be closed, an additional equation is required: the equation of state (EOS). It corresponds to the modelling of thermodynamic properties and consists of an algebraic relation between thermodynamic variables. Indeed, perturbations of the inlet velocity or of the power density may strongly modify the temperature in the uid and cause phase transition from liquid phase to vapor phase. At this modelling scale, the uid can thus be under liquid, vapor or mixture phases. The issue is here to construct an EOS that models all phases of a uid.

The model used in this study is based on the assumption of local mechanic and thermodynamic equilibria between phases. This means that the phases are assumed to move at the same velocity and that vaporisation, condensation and heat transfer processes are assumed to be instantaneous.

As a consequence, the two-phase ow can be considered a single-phase problem provided the EOS [START_REF] Clerc | Numerical Simulation of the Homogeneous Equilibrium Model for Two-Phase Flows[END_REF][START_REF] Goncalvès | Numerical study of cavitating ows with thermodynamic eect[END_REF][START_REF] Jaouen | Étude mathématique et numérique de stabilité pour des modeles hydrodynamiques avec transition de phase[END_REF][START_REF] Müller | The Riemann problem for the Euler equations with nonconvex and nonsmooth equation of state: construction of wave curves[END_REF][START_REF] Voss | Exact Riemann solution for the Euler equations with nonconvex and nonsmooth equation of state[END_REF] with source terms.

General thermodynamics

In classic thermodynamics, two variables are sucient to represent a thermodynamic state of a pure single-phase uid. This is done by means of an EOS which is a relation between the internal energy, the density and the entropy. In the literature, there exist numerous EOS specic to the uid and to the model which are under consideration. In the case of liquid-vapor phase transitions, the EOS must not only represent the behavior of each pure phase (liquid or vapor), but also model the rate of the mass transfer between one phase to the other.

Phase diagram of Figure 2b represent the coexistence curve p s (T ) which relates the pressure to the temperature when phase change occurs: the plane (T, p) is split by the coexistence curve into two regions in which one phase or another is stable. At any point on such a curve the two phases have equal Gibbs potentials and both phases can coexist.

In the (1/ρ, p) plane (see Figure 2a) this mixture where both phases coexist is called the saturation zone. The designation at saturation means that the steam is in equilibrium with the liquid phase. This region is bounded by two curves connected at the critical point (1/ρ c , p c ) which also belongs to the critical isotherm T = T c . Within the two-phase region, through any point passes an isotherm which is a straight line.

These curves can be obtained experimentally (see [START_REF] Lemmon | Thermophysical Properties of Fluid Systems[END_REF] for instance) and correspond to thermodynamic equilibria of temperature, pressure and Gibbs potential of the two phases. The Van der Waals law associated to the Maxwell construction is the most common example of this kind of EOS.

Nevertheless, it is very complicated to derive a unique EOS describing accurately both pure and mixture phases. To better handle pure phases and saturation curves, an idea consists in using two laws (one for each phase) so that each phase has its own thermodynamics. In the following section, we detail the general construction of the EOS in the mixture region given one EOS for each phase.

Construction of the EOS in the mixture

In this section, we explain how to specify the EOS in the mixture given an EOS for each pure phase.

Characterization of the two-phase media We consider each phase κ (κ = for the liquid phase and κ = g for the vapor phase) as a compressible uid governed by a given EOS: (ρ, ε) → η κ where ρ, ε and η κ denote respectively the specic density, the specic internal energy and the specic entropy of the uid. We assume that the function (τ, ε) → η κ (1/τ, ε) has a negative-denite Hessian matrix where τ def = 1/ρ is the specic volume [START_REF] Callen | Thermodynamics and an Introduction to Thermostatistics[END_REF].

We then dene classically for any phase κ the temperature T κ , the pressure p κ and the chemical potential g κ respectively by

T κ def = ∂η κ ∂ε ρ -1 , p κ def = -ρ 2 T κ ∂η κ ∂ρ ε , g κ def = ε + p κ ρ -T κ η κ .
Finally α denotes the volume fraction of vapor phase. This variable characterizes the volume of vapor in each unit volume: α = 1 means that this volume is completely lled by vapor; similarly, a full liquid volume corresponds to α = 0. Liquid and vapor are thus characterized by their thermodynamic properties.

The mixture density ρ and the mixture internal energy ε are dened by

ρ def = αρ g + (1 -α)ρ , (2.1a) ρε def = αρ g ε g + (1 -α)ρ ε , (2.1b)
where ρ g , ρ , ε g and ε denote respectively vapor/liquid densities and vapor/liquid internal energies. Recalling that the internal energy is connected to the enthalpy by the relation ρh = ρε + p, we can compute the mixture enthalpy h when the pressure is the same in both phases (which is the case in the Lmnc model where the thermodynamic pressure p is constant and equal to p 0 ). This leads to ρh = αρ g h

g + (1 -α)ρ h , (2.2)
where h g , h are respectively vapor/liquid enthalpies.

When taking phase transition into account, the two-phase mixture is constructed according to the second principle of thermodynamics. The key idea is that, when phases coexist (i.e. when 0 < α < 1), they have the same pressure, the same temperature and their chemical potentials are equal. The corresponding temperature, noted T s for temperature at saturation, is obtained by expliciting the equality of chemical potentials g (p, T s ) = g g (p, T s ). This implies a relation between T s and p (see for example [START_REF] Callen | Thermodynamics and an Introduction to Thermostatistics[END_REF][START_REF] Faccanoni | Étude d'un modèle n de changement de phase liquide-vapeur[END_REF][START_REF] Greiner | Thermodynamics and statistical mechanics[END_REF] for more details). In the sequel we choose to express the temperature in the mixture as a function of the pressure and we dene functions at

saturation ρ s κ and h s κ by p → ρ s κ def = ρ κ p, T s (p) and p → h s κ def = h κ p, T s (p)
. Consequently, all thermodynamic quantities can be expressed as functions of the enthalpy and the pressure as it will be seen below. The choice to focus on pressure relies on the fact that the pressure in the Lmnc model is supposed to be constant and equal to p 0 .

Remark 2.1. Notice that for most uids h s < h s g and ρ s > ρ s g . In fact, the dierence of enthalpies between the two phases in the saturated mixture is the latent heat of vaporisation

L s g (p) def = h s g (p) -h s (p)
and we have the Clapeyron's law linking the latent heat of vaporisation to the slope of the coexistence curve (see [START_REF] Callen | Thermodynamics and an Introduction to Thermostatistics[END_REF][START_REF] Greiner | Thermodynamics and statistical mechanics[END_REF] for more details):

L s g (p) = 1 ρ s g (p) -1 ρ s (p) T s (p) (T s ) (p)
.

Density of the two-phase media Given functions at saturation, we are now able to model density in pure and mixture phases. Using equations (2.1a) and (2.2), the density is written as a function of enthalpy h and pressure p as follows

ρ(h, p) =          ρ (h, p), if h ≤ h s (p), ρ m (h, p) = ρ s g ρ s (h s g -h s ) (p) [ρ s g h s g -ρ s h s ](p) -h • [ρ s g -ρ s ](p) , if h s (p) < h < h s g (p), ρ g (h, p), if h ≥ h s g (p). (2.3)
For the derivation of this formula, see Appendix A. Notice that we have

α(h, p) =          0, if h ≤ h s (p), ρ s (p)[h -h s (p)] [ρ s g h s g -ρ s h s ](p) -h • [ρ s g -ρ s ](p) , if h s (p) < h < h s g (p), 1, if h ≥ h s g (p).
(2.4)

Temperature of the two-phase media The temperature in the mixture T s is implicitly dened by the equation g (p, T s ) = g g (p, T s ) so that the temperature depends continuously on the enthalpy and on the pressure and reads

T (h, p) =      T (h, p), if h ≤ h s (p), T s (p), if h s (p) < h < h s g (p), T g (h, p), if h ≥ h s g (p).
(2.5)

We must emphasize that the function h → T (h, p = p 0 ) cannot be inverted in the mixture zone for a constant pressure (as it is the case in the Lmnc model). This remark prevents from working with equations on T instead of equations on h.

Compressibility coecient of two-phase media Computing the derivative of the density (2.3), we obtain the compressibility coecient β (previously dened by (1.2))

β(h, p) def = - p ρ 2 • ∂ρ ∂h p =          β (h, p), if h ≤ h s (p), β m (p) def = p • 1 ρ s g -1 ρ s h s g -h s (p), if h s (p) < h < h s g (p), β g (h, p), if h ≥ h s g (p). (2.6)
We notice that independently of the EOS in the pure phases, the compressibility coecient is constant in the mixture (since the pressure is constant in the Lmnc model). Moreover, it is generally discontinuous between pure and mixture phases.

The stiened gas EOS

Several EOS can be considered to describe thermodynamic properties of pure phases. In this paper (like in [START_REF] Bernard | Study of low Mach nuclear core model for single-phase ow[END_REF][START_REF] Dellacherie | On a low Mach nuclear core model[END_REF]) we use the stiened gas law. This EOS is the simplest prototype that contains the main physical properties of pure uids such as repulsive and attractive molecular eects, thereby facilitating the handling of thermodynamics through a simple analytical formulation. It is a generalization of the well-known ideal gas law (which is a commonly used EOS to describe the vapor phase), and it is an acceptable model for the liquid phase which is nearly incompressible (see Appendix C for water and steam parameters). Moreover, we will see that the EOS which models the mixture is a stiened gas EOS regardless of the EOS modelling the pure phases.

For a given pure phase, the stiened gas EOS is fully dened by the entropy η written as a function of the density ρ and the internal energy ε:

(ρ, ε) → η = c v [ln(ε -q -π/ρ) -(γ -1) ln ρ] + m.
(2.7)

The parameters c v > 0 (specic heat at constant volume), γ > 1 (adiabatic coecient), π (constant reference pressure), q (binding energy) and m (reference entropy, relevant only when phase transition is involved) are some constants describing thermodynamic properties of the phase. Note that the case of an ideal gas is recovered by setting π and q to zero. We refer to [START_REF] Meniko | The Riemann problem for uid ow of real materials[END_REF] for a more in-depth discussion on the physical basis for this EOS.

The classic denitions in thermodynamics provide the following expressions for the temperature T , the pressure p, the enthalpy h and the Gibbs potential g as functions of the density ρ and the internal energy ε:

p(ρ, ε) def = -T ρ 2 ∂η ∂ρ ε = (γ -1)(ε -q -π/ρ)ρ -π = (γ -1)(ε -q)ρ -γπ, T (ρ, ε) def = ∂η ∂ε ρ -1 = ε -q -π/ρ c v , h(ρ, ε) def = ε + p ρ = q + (ε -q -π/ρ)γ, g(ρ, ε) def = ε -T η + p ρ = q + (ε -q -π/ρ) γ - m c v -ln (ε -q -π/ρ)ρ 1-γ .
Physical considerations We underline that the various parameters of the stiened gas EOS cannot be chosen freely if a physically correct thermodynamic behavior is expected. Throughout this paper, we will consistently make the assumption that the parameters satisfy the following standard restrictions, which follow from thermodynamic stability theory.

For η to be well-dened, it is necessary to have ε -q -π/ρ > 0 and ρ > 0. Since c v > 0, the rst inequality implies T > 0. Because of γ > 0, we get h -q > 0. Moreover, we also have to satisfy ρ > 0: since γ > 1 and h -q > 0, this is satised when p + π > 0. We note that this hypothesis does not generally guarantee positivity of the pressure. This is consistent with the view that a stiened gas is obtained by shifting the zero point of an ideal gas pressure [START_REF] Meniko | The Riemann problem for uid ow of real materials[END_REF]. In particular, all derived thermodynamic quantities are well dened as long as p + π remains positive; see for instance [START_REF] Faccanoni | Modelling and simulation of liquid-vapor phase transition in compressible ows based on thermodynamical equilibrium[END_REF]. Hence there is no reason to discard negative-pressure solutions as unphysical. To summarize, the modelling hypotheses upon the EOS are

c v > 0, γ > 1 and p + π > 0.
(2.8)

The term (γ -1)(ε -q)ρ > 0 in the expression of p models repulsive eects that are present for any state (gas, liquid or solid) and is due to molecular vibrations. π leads to the stiened properties compared to an ideal gas: a large value of π implies near-incompressible behavior.

The product γπ > 0 represents the attractive molecular eect that guarantees the cohesion of matter in liquid or solid phases (hence π = 0 for a gas).

Temperature and enthalpy at saturation We assume that each phase κ is described by its own stiened gas EOS. To complete the results from Section 2.2, we have to express the temperature at saturation. As the temperature is constant in the mixture, we make a change of thermodynamic variables from (ρ, ε) to (p, T ) which can be made explicit for this kind of EOS. The variables are now given by

ρ κ (p, T ) = p + π κ (γ κ -1)c vκ T , ε κ (p, T ) = c vκ T p + π κ γ κ p + π κ + q κ , h κ (p, T ) = q κ + γ κ c vκ T, g κ (p, T ) = q κ + T c vκ γ κ -q κ -c vκ γ κ ln T + c vκ (γ κ -1) ln(p + π κ ) ,
where for the sake of simplicity we denoted

q κ def = m κ + c vκ γ κ ln c vκ + c vκ (γ κ -1) ln(γ κ -1) (2.9)
as in [START_REF] Métayer | Elaborating equations of state of a liquid and its vapor for two-phase ow models[END_REF][START_REF] Saurel | Modelling phase transition in metastable liquids: application to cavitating and ashing ows[END_REF]. We are now able to dene the temperature at saturation T s (p) of the mixture as the solution of the equation g (p, T s ) = g g (p, T s ) which yields

(c vg γ g -c v γ ) 1-ln T s (p) +c vg (γ g -1) ln(p+π g )-c v (γ -1) ln(p+π ) = q -q g T s (p) +q g -q . (2.10)
We suppose in the sequel that for κ ∈ { , g}, c vκ , γ κ and π κ satisfy the modelling hypothesis (2.8), and that q κ and q κ are such that T s (p) exists and is unique at least when p = p 0 (this is the case for the constants of Table 1 in Appendix C computed for liquid and steam water). Thus, we have in particular T s (p 0 ) > 0. We remark that if q = q g or if c vg γ g = c v γ , we can compute T s analytically. Otherwise, a Newton algorithm can be used to solve this nonlinear equation for any xed p. We then deduce the enthalpy at saturation for each phase

h s κ (p) = q κ + γ κ c vκ T s (p). (2.11)
Density The density is linked to the enthalpy by relation (2.3) where

ρ κ (h, p) = γ κ γ κ -1 p + π κ h -q κ , (2.12a)
ρ s κ (p) = p + π κ (γ κ -1)c vκ T s (p)
.

(

The density ρ s κ (p) denes the density at saturation for each phase κ ∈ { , g}.

Temperature The temperature satises relation (2.5) with

T κ (h, p) = h -q κ γ κ c vκ .
Compressibility coecient Relation (2.6) provides the expression of the compressibility coecient with

β(h, p) =                    β (p) = γ -1 γ p p + π , if h ≤ h s (p), β m (p) = p • 1 ρ s g -1 ρ s h s g -h s (p), if h s (p) < h < h s g (p), β g (p) = γ g -1 γ g p p + π g , if h ≥ h s g (p).
(2.13)

We notice that β κ is independent from h whereas β depends on h through the choice of the phase κ ∈ { , m, g}.

Binding energy We set

q(h, p) def =          q , if h ≤ h s (p), q m (p) def = ρ s g h s g -ρ s h s ρ s g -ρ s (p), if h s (p) < h < h s g (p), q g , if h ≥ h s g (p).
(2.14)

We notice that h s (p) > q , h s g (p) > q g and h s g (p) > h s (p) > q m (p).

(2.15)

The two rst estimates result from the denition (2.11) of h s κ while the last one is proved by a straightforward calculation as (h s g -h s )(ρ s g -ρ s ) < 0 (see Figure 3 and remark 2.1).

Graphs of density, temperature, compressibility coecient and speed of sound (whose expression is detailed in Sect. B) for liquid water and steam at p = 155 × 10 5 Pa with parameters of Table 1 (page 49) are pictured on Figure 3.

Conclusion: a unied stiened gas EOS Given relation (2.3), the density can be expressed by

ρ(h, p) = p/β(h, p) h -q(h, p) (2.16)
where β(h, p) and q(h, p) are given by (2.13) and (2.14). It is important to note that for any EOS dening the pure phases, ρ(h, p) in the mixture is always given by ρ m (h, p)

= p/βm(p) h-qm(p) .
Thus, since p is constant in the Lmnc model, the mixture can always be considered a stiened gas. Of course, this important property is mostly due to the local mechanic and thermodynamic equilibria hypothesis which gives ρ m (h, p) (see Appendix A).

Because of (2.16), PDE (1.1b) can be rewritten as

∂ t h + v∂ y h = β(h, p 0 ) p 0 h -q(h, p 0 ) Φ.
(2.17)

This formulation is the key point of the present study and will be used in the sequel instead of (1.1b).

Theoretical study

In this section we derive some analytical steady and unsteady solutions to system (1.1) together with BC (1.3) with stiened gas law so that equation (1.1b) is replaced by (2.17). For a single phase ow we obtain exact and asymptotic solutions for dierent power densities and inlet velocities. We then extend these calculations to two-phase ows with phase transition. We point out that our results generalize earlier works from Gonzalez-Santalo and Lahey [START_REF] Gonzalez-Santalo | An exact solution for ow transients in two-phase systems by the method of characteristics[END_REF]. In the latter paper, although the modelling of mixture is identical to what we present here, pure phases are considered incompressible. On the contrary, our work does not rely on any restriction in pure phases which allows for a physically more relevant modelling especially when pure gas phase appears. the enthalpies of the liquid and vapor phases at saturation are h s ≈ 1.627 × 10 6 J • K -1 and h s g ≈ 3.004 × 10 6 J • K -1 .

Validity of the model

As stated in 2.3, we must ensure of the positivity of h -q(h, p 0 ) for h solution to (2.17). We mention that in the framework of stiened gas satisfying (2.8), Hypotheses 1.2(2.) and 1.4(2.) reduce to: Hypothesis 3.1. Data h e and h 0 are such that

H def = min inf t≥0 h e (t), min y∈[0,L] h 0 (y) > q .
Indeed, if at instant t the inlet ow is in liquid phase (i.e. h e (t) < h s ) and satises Hypothesis 3.1, then h e (t) -q(h e (t), p 0 ) = h e (t) -q > 0. If the inlet ow is a mixture of liquid and steam (i.e. h e (t) ∈ [h s , h s g ]), then h e (t) -q(h e (t), p 0 ) = h e (t) -q m ≥ h s -q m > 0 according to (2.15). Likewise, if the inlet ow is in vapor phase (i.e. h e (t) > h s g ), then h e (t) -q(h e (t), p 0 ) = h e (t) -q g > h s g -q g > 0 according to (2.15). Thus, in each case, h e (t) -q(h e (t), p 0 ) > 0 which implies that ρ e (t) is well-dened and positive through (2.8) and (2.12a). The same proof applies to h 0 .

To go further, we need to dene some notations that will be useful for the whole section. To solve the transport equation (2.17), we make use of the method of characteristics. This method consists in constructing curves (called characteristic curves) along which PDE (2.17) reduces to an ordinary dierential equation (ODE). More precisely, let χ(τ ; t, y) be the position at time τ of a particle located in y at time t in a ow driven at velocity v. For t ≥ 0 and y ∈ (0, L), χ is thus solution to the parametrized ODE

   dχ dτ (τ ; t, y) = v τ, χ(τ ; t, y)), (3.1a)
χ(t; t, y) = y.

(3.1b)

The curve τ, χ(τ ; t, y) is the characteristic curve passing through the point (t, y). Properties of χ depend on the smoothness of the velocity eld v.

We denote in the sequel (see Figure 4) ξ(t, y) the foot of the characteristic curve, t * (t, y) the time at which the characteristic curve crosses the boundary y = 0 and y * (t) the location at time t of a particle initially placed at y = 0. In other words, ξ, t * and y * are dened by ξ(t, y) def = χ(0; t, y), χ t * (t, y); t, y = 0 and ξ t, y * (t) = 0.

(3.2)
We thus have the characterization

ξ(t, y) > 0 ⇐⇒ t * (t, y) < 0 ⇐⇒ y > y * (t).
This enables to prove the following result: Lemma 3.1. When Hypothesis 3.1 is satised, any smooth solution to PDE (2.17) with BC (1.3a) and well-prepared initial conditions satises h -q(h, p 0 ) > 0.

Proof of Lemma 3.1. We recall that q is dened through (2.14). For any t ≥ 0 and y ∈ (0, L), there are three possible situations:

If h(t, y) > h s g , then h(t, y) -q h(t, y), p 0 > h s g -q g > 0 because of (2.15);

If h(t, y) ∈ [h s , h s g ], then h(t, y) -q h(t, y), p 0 ≥ h s -q m > 0 because of (2.15); If h(t, y) < h s , then equation (2.17) reads (at least locally by a smoothness argument)

∂ t (h -q ) + v∂ y (h -q ) = β Φ p 0 (h -q ).
We infer that ĥ : (τ ; t, y) → h τ, χ(τ ; t, y) satises

    
∂ τ ĥ(τ ; t, y) -q = ĥ(τ ; t, y) -q β p 0 Φ τ, χ(τ ; t, y) , ĥ(t; t, y) = h(t, y).

Hence h(t, y) -q = ĥ(τ ; t, y) -q exp β p 0 t τ Φ σ, χ(σ; t, y) dσ

for any τ such that χ(τ ; t, y) ∈ (0, L).

This shows that h(t, y) -q is continuous and has the same sign as h e t * (t, y) -q (if ξ(t, y) ≤ 0),

or as h 0 ξ(t, y) -q (if ξ(t, y) ∈ [0, L]),
or as h s -q (if ξ(t, y) > L or if ξ(t, y) does not exist, i.e. when phase change occurs before and after).

All of them are positive thanks to Hypothesis 3.1 and (2.15).

Exact and asymptotic solutions for single-phase ow

In this section, we compute some analytical solutions of (1.1) supplemented with the stiened gas law for some particular cases (according to relevant values for Φ, h e and D e ) when a single phase κ ∈ { , m, g} is present. The compressibility coecient β and the coecient q are thus constant.

Since we focus on the 1D case and as β κ (h, p 0 ) = β κ in the case of the stiened gas law, we can compute the velocity v by a direct integration of equation (1.1a), which gives

v(t, y) = v e (t) + β κ p 0 y 0 Φ(t, z) dz, (3.3)
with v e dened by (1.4). This velocity is obviously nonnegative under Hypotheses 1.1, 1.2 and 1.3, so that it is compatible with the location of BC. As mentioned earlier, equation (1.1b) can be rewritten as (2.17). To compute the enthalpy, we apply the method of characteristics (see above) which gets simpler in the case of the stiened gas law.

The following results had rst been stated in [START_REF] Bernard | Study of low Mach nuclear core model for single-phase ow[END_REF]. The proofs are detailed below.

y L t ξ t 1 , y * ( t 1 ) = 0 t 1 y * (t 1 ) ξ ( t 1 , y 1 )
y 1 t * ξ ( t 1 , y 2 ) y 2
Figure 4: Sketch of the method of characteristics and denitions of ξ(t, y), t * (t, y) and y * (t).

Constant power density

h(t, y) =          q κ + h 0 ξ(t, y) -q κ e Φ0t ,
if ξ(t, y) ≥ 0, where H 0 (y) = 1/(h 0 (y) -q κ ), H 0 (y) = y/(h 0 (y) -q κ ), H e (t) = 1/(h e (t) -q κ ) and H e (t) = e -Φ0t /(h e (t) -q κ ).

q κ + h e t * (t, y) -q κ 1 + Φ0 y v e = h
Notice that Equation (3.5c) is a correction of Equation (9b) in [START_REF] Bernard | Study of low Mach nuclear core model for single-phase ow[END_REF].

Remark 3.1. As it has been stated in [START_REF] Bernard | Study of low Mach nuclear core model for single-phase ow[END_REF], if the inlet enthalpy is also constant and if the inlet velocity is nonzero, there is an asymptotic state which is reached in nite time. This time is equal to t ∞ = 1 Φ0 ln(1 + Φ0L ve ) and satises ξ(t ∞ , L) = 0 and y * (t ∞ ) = L. Hence, for t ≥ t ∞ , the solution (h, v, p)(t, y) is given by

               h ∞ (y) = h e + Φ 0 D e y, v ∞ (y) = v e + Φ0 y, p ∞ (y) = Φ0 [µ(y) -µ(L)] + gD e Φ0 ln   1 + Φ0L ve 1 + Φ0y ve   + Φ0 D e (L -y).
However, if the inlet velocity is v e = 0, then ξ(t, y) is always positive and no asymptotic state can be reached since the enthalpy increases continuously in time.

Remark 3.2. Proposition 3.1 applies for nonzero Φ 0 . When Φ 0 = 0, the same proof holds except the resolution of ODE (3.1) whose solution becomes χ(τ ; t, y) = y + v e (τ -t). We remark that solution (3.7) converges to this case when Φ0 → 0, which shows that the model is continuous with respect to Φ 0 . In particular, the enthalpy reads

h(t, y) = h 0 (y -v e t),
if y ≥ v e t, h e (t -y/v e ), otherwise, which was expected as h is a solution to a simple linear transport equation.

Proof of Proposition 3.1. Equation (3.3) becomes v(t, y) = v e + Φ0 y.

(3.6) Since v is linear, the Cauchy-Lipschitz theorem applied to ODE (3.1) ensures the existence of χ over some interval (depending on t and y). Moreover, χ is continuous with respect to (τ, t, y).

We then solve ODE (3.1) using expression (3.6). We obtain χ(τ ; t, y) = y + v e Φ0 e Φ0(τ -t) -v e Φ0 .

(3.7)

For xed t ≥ 0 and y ∈ (0, L), the requirement χ(τ ; t, y)

∈ (0, L) constrains the interval of existence τ ∈   max {0; t * (t, y)} , t + 1 Φ0 ln   1 + Φ0L ve 1 + Φ0y ve     . (3.8)
For (τ, t, y) satisfying (3.8), we note ĥ(τ ; t, y) def = h τ, χ(τ ; t, y) .

(3.9)

We deduce from equation (1.1b) rewritten under the form equation (2.17) that ĥ satises ∂ τ ĥ(τ ; t, y) -q κ = Φ0 ĥ(τ ; t, y) -q κ , ĥ(t; t, y) = h(t, y).

The solution of this linear rst order ODE is h(t, y) -q κ = ĥ(t; t, y) -q κ = ĥ(τ ; t, y) -q κ e Φ0(t-τ ) = h τ, χ(τ ; t, y) -q κ e Φ0(t-τ ) . (3.10)

Two cases must be investigated depending on the minimal time until which the characteristic curve remains in the domain or equivalently depending on the sign of ξ (see (3.8) and Figure 4): if ξ(t, y) ≥ 0, then the characteristic curve does not cross the boundary y = 0, which means that we can take τ = 0 in equation (3.10) and h(t, y) -q κ = h 0 χ(0; t, y) -q κ e Φ0t = h 0 ξ(t, y) -q κ e Φ0t ;

if ξ(t, y) < 0, the backward characteristic curve reaches the boundary at time t * (t, y) > 0 and h(t, y) -

q κ = [h e (t * ) -q κ ]e Φ0(t-t * ) = [h e (t * ) -q κ ] 1 + Φ0 y v e .
Noticing that βκ p0 (h

-q κ ) = 1 ρ leads to h(t, y) = q κ + [h e (t * ) -q κ ] + 1 ρ e (t * ) Φ 0 v e y = h e (t * ) + Φ 0 D e (t * ) y.
Proof of Corollary 3.1. The exact dynamic pressure p can be computed by integrating the momentum equation (1.1c) which is equivalent to

∂ y p = -∂ t (ρv) -∂ y (ρv 2 ) + ∂ y (µ∂ y v) -ρg.
Using the mass conservation law and observing that v given by (3.6) is independent of time, we obtain

∂ y p = -( Φ0 v + g)ρ + Φ0 ∂ y µ from which we deduce ∂ y p = - p 0 β κ g + Φ0 (v e + Φ0 y) h -q κ + Φ0 ∂ y µ.
Integrating between y and L, we get due to BC (1.3c)

p(t, y) = p 0 (g + Φ0 v e ) β κ L y 1 h(t, z) dz + p 0 Φ2 0 β κ L y z h(t, z) dz + Φ0 [µ(y) -µ(L)] (3.11)
where h def = h -q κ . Since h is dened piecewise by (3.4), we have to consider three cases:

If y * (t) > L, the curve τ → χ(τ ; t, y) lies in the green part of the graph on Figure 4, i.e.

above the curve τ → χ τ ; t, y * (t) . As y * (t) > L =⇒ χ(t, y) < 0 for all (t, y) ∈ R + × (0, L),

we can select the relevant value for h in (3.4), i.e. h(t, y) = h e t * (t, y) 1 + Φ0y ve .

To compute each integral in (3.11), we use the change of variables τ = t * (t, z), which yields

L y 1 h(t, z) dz = -v e t * (t,L) t * (t,y) 1 h e (τ ) dτ = v e H e (τ ) t * (t,y) t * (t,L) , L y z h(t, z) dz = v 2 e Φ0 t * (t,y) t * (t,L) e Φ0(t-τ ) -1 h e (τ ) dτ = v 2 e Φ0 e Φ0t H e (τ ) -H e (τ ) t * (t,y) t * (t,L)
.

We then infer (3.5a).

If y * (t) ≤ L and y ≥ y * (t) (=⇒ ξ(t, y) ≥ 0), we have h(t, y) = h 0 ξ(t, y) e Φ0t so that by means of the change of variable ζ = ξ(t, z) we specify the integrals in (3.11)

L y 1 H(t, z) dz = ξ(t,L) ξ(t,y) 1 h 0 (ζ) dζ = H 0 (ζ) ξ(t,L) ξ(t,y) , L y z h(z) dz = e Φ0t ξ(t,L) ξ(t,y) ζ h 0 (ζ) dζ + (e Φ0t -1) v e Φ0 ξ(t,L) ξ(t,y) 1 h 0 (ζ) dζ = e Φ0t H 0 (ζ) ξ(t,L) ξ(t,y) + (e Φ0t -1) v e Φ0 H 0 (ζ) ξ(t,L) ξ(t,y) .
Hence we deduce (3.5b).

The last corresponds to y < y * (t) ≤ L. The integration domain is split into two parts depending on the sign of ξ(t, z). More precisely

L y f (z) dz = L y * (t) f (z) dz + y * (t) y f (z) dz
for any function f . For integrals between y * (t) and L, we apply Formula (3.5b) with y = y * (t). For integrals between y and y * (t), we apply Formula (3.5a) with L replaced by y * (t). Summing all terms leads to (3.5c).

Varying power density (with y or t)

Let us generalize the results of Proposition 3.1 by taking Φ = Φ(y): Then ξ and t * dened by (3.2) are equal to ξ(t, y) = Θ -1 Θ(y) -t and t * (t, y) = t -Θ(y). The solution h of equation (1.1b) with BC (1.3a) is given by

h(t, y) =            q κ + v(y) h 0 ξ(t, y) -q κ v ξ(t, y) , if ξ(t, y) ≥ 0, q κ + v(y) h e t * (t, y) -q κ v e = h e t * (t, y) + 1 D e t * (t, y) y 0 Φ(z) dz, if ξ(t, y) < 0.
Let us note that v(y) = v 0 (y) since v(t, y) = v(y) and as v 0 is well-prepared (see 1.2).

We can also extend the results of Proposition 3.1 by taking Φ = Φ(t):

h(t, y) = q κ +    h 0 ξ(t, y) -q κ e Ψ(t) ,
if ξ(t, y) ≥ 0, h e t * (t, y) -q κ e Ψ(t)-Ψ(t * ) , if ξ(t, y) < 0.

Remark 3.3. In the general case, i.e. with time dependence for h e and D e and time/space dependence for Φ, it does not seem possible to compute an exact solution. However, if h e (t), D e (t) and Φ(t, y) have a nite limit as t → +∞ denoted by h ∞ e , D ∞ e , Φ ∞ (y), then there exists a steady solution for the enthalpy

h ∞ (y) = h ∞ e + 1 D ∞ e y 0 Φ ∞ (z) dz
and all other quantities are deduced from h ∞ . We recover the result of [START_REF] Dellacherie | On a low Mach nuclear core model[END_REF].

Up to now, there is no theoretical result in the general case about the convergence of the unsteady solution to this steady state. We just notice that solutions from Props. 3.2 and 3.3 actually do.

Proof of Proposition 3.2. Since the inlet velocity v e is a positive constant and as Φ does not depend on t, the velocity v(t, y) = v(y) = v 0 (y) is a positive function independent of time due to (3.3). Function Θ introduced in the statement of Proposition 3.2 is thus invertible. In a similar way as for Proposition 3.1, we solve the characteristic ODE (3.1) which is equivalent to

1 = χ v(χ) = Θ (χ)χ = [Θ(χ)] .
For the sake of simplicity, we set χ = dχ dτ and χ = d 2 χ dτ 2 . This leads to Θ χ(τ ; t, y) -Θ(y) = Θ χ(τ ; t, y) -Θ χ(t; t, y) = τ -t and χ(τ ; t, y) = Θ -1 Θ(y) + τ -t .

(3.12)

We can ensure that χ(τ ; t, y) ∈ (0, L) provided τ ∈ max 0, t * (t, y) , t + Θ(L) -Θ(y) . Keeping the same notation (3.9) for ĥ, equation (2.17) becomes    ∂ τ ĥ(τ ; t, y) -q κ = β κ p 0 ĥ(τ ; t, y) -q κ Φ χ(τ ; t, y) , ĥ(t; t, y) = h(t, y).

(3.13)

We dierentiate ODE (3.1) to obtain χ (τ ; t, y) = χ (τ ; t, y) dv dy χ(τ ; t, y) = χ (τ ; t, y) β κ p 0 Φ χ(τ ; t, y) .

The positivity of v implies that χ (τ ; t, y) > 0 and

β κ p 0 Φ χ(τ ; t, y) = χ (τ ; t, y) χ (τ ; t, y) = ln χ (τ ; t, y) = ln v χ(τ ; t,

y) .

Inserting this relation in equation (3.13), we have ∂ τ ( ĥ -q κ ) = ( ĥ -q κ )∂ τ ln v χ(τ ; t, y) .

Hence h(t, y) -q κ = v(y) h τ, χ(τ ; t, y) -q κ v χ(τ ; t, y) .

Given expression (3.12) for χ, we nally obtain

h(t, y) = q κ +          v(y) h 0 ξ(t, y) -q κ v ξ(t, y) , if ξ(t, y) ≥ 0, v(y) h e t * (t, y) -q κ v e , otherwise.
Proof of Proposition 3.3. As previously, the key point is the integration of the characteristic ODE (3.1) which reads

d dτ χ(τ )e -Ψ(τ ) = v e (τ )e -Ψ(τ ) .
This leads to χ(τ ; t, y) = ye Ψ(τ )-Ψ(t) + τ t v e (s)e Ψ(τ )-Ψ(s) ds.

Consequently h(t, y) -q κ = h τ, χ(τ ; t, y) -q κ e Ψ(t)-Ψ(τ ) .

The main issue is then to determine the interval for τ such that χ(τ ; t, y) ∈ (0, L). The equality χ(τ ; t, y) = 0 can be rewritten as y = t τ v e (s)e Ψ(t)-Ψ(s) ds.

As the right hand side vanishes for τ = t, two cases may occur: either y is greater than the right hand side for τ = 0 (which would imply that the left bound of the interval is 0), or there exists τ = t * > 0 such that the previous equality holds. In the former case, we obtain h(t, y) -q κ = h 0 χ(0; t, y) -q κ e Ψ(t) , while in the latter case h(t, y) -q κ = h e t * (t, y) -q κ e Ψ(t)-Ψ(t * ) .

Exact and asymptotic solutions for two-phase ow with phase transition

In the case of a two-phase ow with phase transition and if we can determine the position of each phase, we can deduce the exact solution using the single-phase ow results given in Section 3.2.

In the particular case where all parameters and boundary/initial data are constant, the result is the following. Let us assume that h e (t) ≡ h e > q , v e (t) ≡ v e > 0, Φ(t, y) ≡ Φ 0 and h 0 (y) ≡ h 0 = h e > q . Let Φκ def = β κ Φ 0 /p 0 , where κ is respectively , m, g in the liquid, mixture and gas phases. We suppose that the initial and boundary data correspond to the liquid phase, i.e. h e (= h 0 ) < h s . We set

y s def = D e Φ 0 (h s -h e ), y s g def = D e Φ 0 (h s g -h e ).
Let us dene three curves in R + × R + as pictured on Figure 5b t Let us also dene

t s def = t (y s ) = 1 Φ ln h s -q h 0 -q , t s g def = t m (y s g ) = t s + 1 Φm ln h s g -q m h s -q m .
Then the spatiotemporal domain R + ×R + consists of three regions corresponding to liquid, mixture and vapor phases as follows (see Figure 5a):

L = (t, y) ∈ R + × R + t ≤ t s or y ≤ y s , M = (t, y) ∈ (t s , +∞) × (y s , +∞) t ≤ t s g or y ≤ y s g , G = (t, y) ∈ (t s g , +∞) × (y s g , +∞) ;
the solution v of equation (1.1a) is given by

v(t, y) =      v e + Φ y, if (t, y) ∈ L, v e + Φ y s + Φm (y -y s ), if (t, y) ∈ M,
v e + Φ y s + Φm (y s g -y s ) + Φg (y -y s g ), if (t, y) ∈ G, and the solution h of equation (1.1b) is given by (see Figure 5b)

h(t, y) =          q + (h 0 -q )e Φ t ,
if (t, y) ∈ L and t < t (y), q m + (h s -q m )e Φm(t-t s ) , if (t, y) ∈ M and t < t m (y), q g + (h s g -q g )e Φg(t-t s g ) , if (t, y) ∈ G and t < t g (y),

h e + Φ0 De y,
otherwise.

Remark 3.4. We notice that constant Φm is the Zuber's characteristic reaction frequency dened in [START_REF] Zuber | Flow excursions and oscillations in boiling, two-phase ow systems with heat addition[END_REF] and used in [START_REF] Gonzalez-Santalo | An exact solution for ow transients in two-phase systems by the method of characteristics[END_REF]. Contrary to the latter paper, we do not assume the pure phases to be incompressible. Thus we can extend the Zuber's frequency to liquid and gas phases as a piecewise constant function.

Remark 3.5. In any case, we observe that a steady state is reached. It is given by

h ∞ (y) = h e + Φ 0 D e y
and the velocity is piecewise linear. Moreover, we can determine the time t ∞ at which the steady state (it is thus an asymptotic state in that case) is reached

t ∞ def =      t (L), if y s > L, t m (L), if y s ≤ L ≤ y s g , t g (L), if y s g < L.
The assumptions of Proposition 3.4 may seem restrictive but it can be trivially extended to several cases: to other constant initial and boundary conditions (i.e. mixture or vapor); if Φ = Φ 0 < 0 such that v remains positive, the enthalpy is then still monotone.

If Φ and/or boundary-initial conditions are not constant anymore, the enthalpy is no longer monotone. It is thus dicult to determine regions L, M and G. However, if this can be achieved, we can similarly apply the monophasic results in each region to compute the exact solution. In any case, the following result holds (similar to Remark 3.3):

Remark 3.6. For any inlet velocity v e , inlet ow rate D e and power density Φ which have nite limits in time, let us denote

h ∞ e , D ∞ e , Φ ∞ (y) = lim t→+∞ h e (t)
, D e (t), Φ(t, y) . Then, there exists a steady solution h ∞ (y), v ∞ (y), p ∞ (y) given by

h ∞ (y) = h ∞ e + 1 D ∞ e y 0 Φ ∞ (z) dz, v ∞ (y) = D ∞ e ρ h ∞ (y), p 0 , p ∞ (y) = g L y ρ h ∞ (z), p 0 dz -µ β h ∞ (z), p 0 Φ ∞ (z) p 0 - (D ∞ e ) 2 ρ h ∞ (z), p 0 z=L z=y ,
with ρ(h ∞ , p 0 ) given by (2.3).

Proof of Proposition 3.4. Since ρ e , v e , h 0 and Φ 0 are constant and correspond to the liquid phase, equations (1.1a)-(2.17) are written over a time interval [0, T ) for some T > 0 to be specied later (which corresponds to the rst time another phase appears)

               ∂ y v = Φ , ∂ t h + v∂ y h = Φ (h -q ), v(t, 0) = v e , h(t, 0) = h e , h(0, x) = h 0 = h e .
We can apply Proposition 3.1 to this system which leads to

v(y) = v e + Φ y and h(t, y) =      q + (h 0 -q )e Φ t , if t < t (y), h e + Φ 0 D e y, otherwise,
where the curve t = t (y) corresponds to the characteristic curve coming from (t = 0, y = 0), that is ξ (t, y) = 0 (see Proposition 3.1 for the denition of ξ ). The enthalpy h(t, •) is a monotoneincreasing function consisting (spatially) of a linear part and a constant part at each time. Two situations may occur:

either h e + Φ0 De L ≤ h s : the uid remains liquid indenitely (T = +∞) and the enthalpy is equal to h e + Φ0

De y everywhere as soon as t ≥ t (L); or h e + Φ0

De L > h s : a mixture phase appears. In the latter case, there exists t > 0 and y ∈ (0, L) such that h(t, y) = h s . We then dene T = t s as the solution of h(t s , L) = h s and y s as the smallest y such that h(t s , y) = h s , i.e.

q + (h 0 -q )e Φ t s = h s , h e + Φ 0 D e y s = h s .
For t > t s , the uid is in liquid phase for y < y s and in mixture phase for y ≥ y s and t small enough, i.e. We adapt Proposition 3.1 to the current spatiotemporal domain so that we obtain

               ∂ y v = Φm , ∂ t h + v∂ y h = Φm (h -q m ), v(t, y s ) = v e + Φ y s , h(t, y s ) = h s , h(t s , y) = h s . y 2 y N -1 h n i v n i p n i y i 0 y 1 h n 1 =he(t n ) v n 1 =ve(t n ) L y N p n N = p 0
v(y) = v e + Φ y s + Φm (y -y s ) and h(t, y) =      q m + (h s -q m )e Φm(t-t s ) , if t < t m (y), h e + Φ D e y,
otherwise.

The curve t = t m (y) corresponds to the characteristic curve passing through (t s , y s ), i.e.

t m (y) = 1 Φm ln v e + ( Φ -Φm )y s + Φm y v e + Φ y s + t s .
This can be rewritten in a simpler form. To this end, we use the denition of y s , the relation

1 ρe = β p0 (h e -q
) and the continuity of ρ at y s (implying

β m (h s -q m ) = β (h s -q )) which leads to t m (y) = 1 Φm ln h e -q m + Φ0y De h s -q m + t s .
If needed, we apply the same procedure to nd (t s g , y s g ) such that h s g is reached and the proof follows.

Numerical scheme

The main advantage of dimension 1 is that equations in (1.1) can be decoupled in an explicit numerical procedure which means that it suces to compute h from equation (1.1b) i.e. (2.17) for the stiened gas law in order to deduce all other variables. To solve the transport equation with source term (2.17), the numerical method of characteristics (MOC) seems suitable insofar as this method was essential in the theoretical part of the study. Moreover, this method has interesting properties such as unconditional stability and positivity-preserving for variable h -q (see 4.3). Accuracy can be improved using 2nd-order results from [START_REF] Penel | An explicit stable numerical scheme for the 1D transport equation[END_REF].

Given ∆y > 0 and ∆t > 0, we consider a uniform Cartesian grid { y i = i∆y } 1≤i≤N such that y 1 = 0 and y N = L (see Figure 6) as well as a time discretization { t n = n∆t } n≥0 . Unknowns are collocated at the nodes of the mesh. We set the initial values v 0 i = v 0 (y i ) and h 0 i = h 0 (y i ) for i = 1, . . . , N .

We must emphasize that as the pressure variable p is supposed to be constant and equal to p 0 in the Lmnc model, most parameters are constant throughout the study (such as T s , β κ or q κ ).

That is why references to the dependence on the pressure may be dropped in the sequel for the sake of simplicity.

Key idea of the scheme

For details about numerical methods of characteristics, the reader may refer to [START_REF] Acrivos | Method of characteristics technique. Application to heat and mass transfer problems[END_REF][START_REF] Durran | Numerical methods for uid dynamics[END_REF][START_REF] Penel | An explicit stable numerical scheme for the 1D transport equation[END_REF] and references therein. This method amounts to tracking particles along the ow and to solve ODEs along the characteristic curves. More precisely, it consists in locating accurately the position of particles and updating the unknown according to the source term as described from a theoretical point of view in Section 3.2. Actually, at time t n , we aim at approximating the solution of

d dτ ĥn+1 i (τ ) = β ĥn+1 i (τ ) Φ τ, χ(τ ; t n+1 , y i ) p 0 ĥn+1 i (τ ) -q ĥn+1 i (τ ) , (4.1a)
where τ → ĥn+1 i (τ ) def = h τ, χ(τ ; t n+1 , y i ) and the characteristic ow χ satises

   d dτ χ(τ ; t n+1 , y i ) = v τ, χ(τ ; t n+1 , y i ) , τ ≤ t n+1 , χ(t n+1 ; t n+1 , y i ) = y i , (4.1b)
with v an approximate solution to (1.1a). The reader may refer to [START_REF] Penel | An explicit stable numerical scheme for the 1D transport equation[END_REF] for further details about the resolution of (4.1b).

The classic Euler-type MOC method applied to ODE (4.1a) over the interval [t n ;

t n+1 ] is writ- ten h(t n+1 , y i ) = ĥn+1 i (t n+1 ) = ĥn+1 i (t n ) + ∆t β ĥn+1 i (t n ) Φ t n , χ(t n ; t n+1 , y i ) p 0 ĥn+1 i (t n ) -q ĥn+1 i (t n ) + O(∆t 2 ).
Hence the scheme can be written

h n+1 i = hn+1,n i + ∆t β hn+1,n i Φ t n , ξ n i p 0 hn+1,n i -q hn+1,n i , (4.2)
where ξ n i and hn+1,n i are some approximations of χ(t n ; t n+1 , y i ) and ĥn+1 i (t n ) respectively. In the sequel, scheme (4.2) will be referred to as the MOC scheme.

To reach higher order, we rewrite ODE (4.1a) using the facts that β(h) and h -q(h) are positive under Hypotheses 1.3 and 3.1. We set R( ĥ) def = ĥ

H dh β(h) • h -q(h)
,

where H is dened within Hypothesis 3.1. Then, equation (4.1a) reads

R ĥn+1 i d dτ ĥn+1 i (τ ) = Φ τ, χ(τ ; t n+1 , y i ) p 0 (4.3)
and thus can be integrated explicitly between t n and t n+1

R ĥn+1 i (t n+1 ) -R ĥn+1 i (t n ) = 1 p 0 t n+1 t n Φ τ, χ(τ ; t n+1 , y i ) dτ.
As Φ is a datum, the right hand side can be expanded at any order or exactly computed. Using the trapezoidal rule, the scheme reads

h n+1 i = R -1 R hn+1,n i + ∆t p 0 Φ(t n , ξ n i ) + Φ(t n+1 , y i ) 2 . (4.4)
Notice that we can give explicit expressions for R and R -1 :

R(h) =          1 β ln h-q H-q , if h ≤ h s , R s + 1 βm ln h-qm h s -qm , if h s < h < h s g , R s g + 1 βg ln h-qg h s g -qg , if h ≥ h s g , R -1 (r) =      q + (H -q )e β r , if r ≤ R s , q m + (h s -q m )e βm(r-R s ) , if R s < r < R s g , q g + (h s g -q g )e βg(r-R s g ) , if r ≥ R s g , where R s def = 1 β ln h s -q H -q , R s g def = R s + 1 β m ln h s g -q m h s -q m .
Strategy (4.4) is named INTMOC. We recall that h s κ > q κ and h s g > h s > q m see (2.15) and that H > q under Hypothesis 3.1. Thus, R s and R s g are well-dened. Strategies (4.2) and (4.4) being set up, it remains to specify how to compute ξ n i and hn+1,n i .

These computations are the core of numerical methods of characteristics and are detailed in the next section.

Description of the scheme

Given the numerical solutions (h n i , v n i , p n i ), the overall process at step n + 1 consists in computing successively h n+1 i , v n+1 i and p n+1 i as follows.

• Enthalpy. For the boundary condition (i = 1) we impose h n+1 1 = h e (t n+1 ). Then h n+1 i is determined in two steps:

x Solve ODE (4.1b) over the interval [t n , t n+1 ]. The approximation ξ n i of χ(t n ; t n+1 , y i )
is computed either at order 1 or 2:

(i) at order 1 in time, we have χ(t n ; t n+1 , y i ) ≈ y i -∆t • v(t n , y i ) so that we set

ξ n i = y i -∆t • v n i ; (4.5)
(ii) at order 2 in time (see [START_REF] Penel | An explicit stable numerical scheme for the 1D transport equation[END_REF][START_REF] Penel | Existence of global solutions to the 1D abstract bubble vibration model[END_REF] for more details), we have

χ(t n ; t n+1 , y i ) ≈ y i -∆t • v(t n , y i ) - 1 2 ∆t 2 ∂ t v(t n , y i ) -v(t n , y i )∂ y v(t n , y i ) .
Due to (1.1a) and with a standard 1st-order nite-dierence discretization for ∂ t v, we set

ξ n i = y i -∆t 3 2 v n i - 1 2 v n-1 i + ∆t 2 2 β(h n i ) p 0 v n i Φ(t n , y i ). (4.6)
y Update the enthalpy: if ξ n i > 0 (see Figure 7a), let j be the index such that ξ n i ∈ [y j , y j+1 ) and θ n ij def =(y j+1 -ξ n i )/∆y. As ξ n i is generally not a mesh node, we use an interpo- lation procedure to evaluate hn+1,n i :

(i) at order 1 hn+1,n i = θ n ij h n j + (1 -θ n ij )h n j+1 ; (4.7)
(ii) at higher order (see [START_REF] Penel | An explicit stable numerical scheme for the 1D transport equation[END_REF] for more details)

hn+1,n i = λ n i h - j + (1 -λ n i )h + j (4.8)
where

λ n i def =          1+θ n ij 3 , if P + j (θ n ij ) ≥ 0 and P - j (θ n ij ) ≥ 0, 0, if P + j (θ n ij ) ≥ 0 and P - j (θ n ij ) < 0, 1, if P + j (θ n ij ) < 0 and P - j (θ n ij ) ≥ 0, θ n ij , otherwise, h - j def =          h n j , if P + j (θ n ij ) < 0 and P - j (θ n ij ) < 0, (θ n ij ) 2 2 h n j-1 -2h n j + h n j+1 - θ n ij 2 h n j-1 -4h n j + 3h n j+1 + h n j+1 , otherwise, h + j def =          h n j+1 , if P + j (θ n ij ) < 0 and P - j (θ n ij ) < 0, (θ n ij ) 2 2 h n j+2 -2h n j+1 + h n j - θ n ij 2 h n j+2 -h n j + h n j+1 , otherwise,
and This procedure has been designed in [START_REF] Penel | An explicit stable numerical scheme for the 1D transport equation[END_REF] in order to ensure the maximum principle by means of a variable stencil (and an adaptive order). Even if there is no maximum principle associated to equation (2.17), this scheme preserves the property hn i -q( hn i ) > 0 (see Proposition 4.1).

P ± j (θ) def =(θ -δ ± j )(θ -δ ± j+1 ) with δ - j def = 2(h n j+1 -h n j ) h n j-1 -2h n j + h n j+1 , δ + j def = 2(h n j+1 -h n j ) h n j -2h n j+1 + h n j+2 , δ - j+1 def = h n j-1 -4h n j + 3h n j+1 h n j-1 -2h n j + h n j+1 , δ + j+1 def = h n j+2 -h n j h n j -2h n j+1 + h n j+2 . y t y j-1 y j y j+1 t n+1 y i t n ξ n i (a) ξ n i > 0 y t 0 t n+1 y i t n ξ n i t * i (b) ξ n i ≤ 0
We then update h n+1 i by formulae (4.2) or (4.4).

if ξ n i ≤ 0 (see Figure 7b), we compute the time t * i at which the characteristic curve τ → χ(τ ; t n+1 , y i ) crosses the inow boundary. There we have h(t * i , 0) = h e (t * i ). Using a rst order approximation in time, we set t * i = t n+1 -y i /v n i and we compute the updated enthalpy similarly to what is detailed above:

(i) by integrating ODE (4.1a) over [t * i , t n+1 ] (MOC strategy)

h n+1 i = h e (t * i ) + (t n+1 -t * i ) β h e (t * i ) Φ(t * , 0) p 0 h e (t * i ) -q h e (
h n+1 i = R -1 R h e (t * i ) + t n+1 -t * i p 0 Φ(t * i , 0) + Φ(t n+1 , y i ) 2 . (4.9b)
The boundary y = 0 is the only one we need to care about since characteristic curves cannot exit from the domain at y = L (we assumed that v e > 0 and Φ ≥ 0 which implies that v > 0).

• Velocity. For the boundary condition (i = 1), we set v n+1 = v e (t n+1 ). Then, we integrate equation (1.1a) over [y i , y i+1 ]. Depending on the ability to compute the primitive function of Φ (and as β is piecewise constant), the velocity eld can be computed directly

v n+1 i = v n+1 i-1 + 1 p 0 yi yi-1 β h(t n+1 , z) Φ(t n+1 , z) dz, (4.10)
or approximated for example by the following upwind approach

v n+1 i = v n+1 i-1 + ∆y p 0 β(h n+1 i-1 )Φ(t n+1 , y i-1 ) (4.11)
since h is transported by v n i ≥ 0. However, since the coecient β is discontinuous at phase change points (see Figure 3c), we have to adapt the previous algorithm in cells where the uid changes from a phase to the mixture. It is reasonable to suppose that at most a pure phase and a mixture are present within a single cell (never liquid, mixture and steam simultaneously). Then, if h s κ ∈ (h n+1 i-1 , h n+1 i ), let y * be the linear approximation of y s κ , i.e.

y * def = y i-1 + ∆y h s κ -h n+1 i-1 h n+1 i -h n+1 i-1 . Hence yi yi-1 β h(t n+1 , z) Φ(t n+1 , z) dz = y * yi-1 β h(t n+1 , z) Φ(t n+1 , z) dz + yi y * β h(t n+1 , z) Φ(t n+1 , z) dz (4.10')
or when the primitive function of Φ is not known

yi yi-1 β h(t n+1 , z) Φ(t n+1 , z) dz ≈ (y * -y i-1 )β(h n+1 i-1 )Φ(t n+1 , y i-1 ) + (y i -y * )β(h n+1 i )Φ(t n+1 , y i ). (4.11')
• Pressure. For the boundary condition (i = N ), we set p n N = p 0 . Then we rewrite equation (1.1c) in the following equivalent form

-∂ y p = ∂ t (ρ(h)v) + ∂ y (ρ(h)v 2 ) -∂ y (µ∂ y v) + ρ(h)g = ρ(h)∂ t v + ρ(h)v∂ y v -∂ y (µ∂ y v) + ρ(h)g. Using (1.1a) it becomes -∂ y p = ρ(h)∂ t v + ρ(h)v β(h)Φ p 0 -∂ y µ β(h)Φ p 0 + ρ(h)g. Let us note ρ n+1 i = ρ(h n+1 i ) and β n+1 i = β(h n+1 i ). Integrating this equation over [y i-1 , y i ],
we obtain with well-prepared initial conditions ensure the positivity of h n i -q(h n i ) for any ∆t > 0.

p n+1 i-1 = p n+1 i + ∆y 2 ρ n+1 i + ρ n+1 i-1 g + ρ n+1 i v n+1 i -v n i ∆t + ρ n+1 i-1 v n+1 i-1 -v n i-1 ∆t + ρ n+1 i v n+1 i β n+1 i p 0 Φ(t n+1 , y i ) + ρ n+1 i-1 v n+1 i-1 β n+1 i-1 p 0 Φ(t n+1 , y i-1 ) -µ β n+1 i p 0 Φ(t n+1 , y i ) - β n+1 i-1 p 0 Φ(t n+1 , y i-1 ) , i ∈ {2, . . . , N }.
Proof. For the sake of clarity, we set κ(h) the phase corresponding to the value of h, namely

κ(h) = if h < h s , m if h ∈ [h s , h s g ] and g if h > h s g .
Let us show by induction the inequalities

h n i ≥ H and h n i > q κ(h n i ) . (4.13)
We recall that H is dened within Hyp. 3.1. Given the physical framework, we assume without loss of generality that the uid is under the liquid phase initially and at the entry, which means h e (t) < h s and h 0 (y) < h s for all t ≥ 0 and y ∈ (0, L). Estimates (4.13) is satised for n = 0 by hypothesis. If it is satised at iteration n, then hn+1,n i ≥ H according to the maximum principle veried by the interpolation processes (4.7) (by convexity) and (4.8) (by construction: see [START_REF] Penel | An explicit stable numerical scheme for the 1D transport equation[END_REF]).

Then, if κ( hn+1,n i ) ∈ {m, g}, we directly have hn+1,

n i > q κ( hn+1,n i ) due to (2.15). If κ( hn+1,n i ) = , then hn+1,n i ≥ H > q due to Hyp. 3.1.
To summarize, we have hn+1,

n i > q κ( hn+1,n i ) .
To complete the proof, we remark that no matter what scheme is used to compute h n+1 i from hn+1,n i , the cases κ(h n+1 i ) ∈ {m, g} are trivially handled due to (2.15). It thus remains to deal with κ(

h n+1 i ) = . If h n+1 i is computed via the MOC scheme (4.2), as hn+1,n i > q κ( hn+1,n i ) , we remark that h n+1 i > hn+1,n i , which implies that κ( hn+1,n i ) = . Thus, (4.2) yields h n+1 i -q = ( hn+1,n i -q ) 1 + ∆t β Φ(t n , ξ n i ) p 0 > 0.
h n+1 i -q = ( hn+1,n i -q ) exp ∆t β Φ(t n , ξ n i ) + Φ(t n+1 , y i ) 2p 0 .
Hence the conclusion, which also holds in the exiting characteristic curve case via schemes (4.9a) and (4.9b). Both schemes are thus compatible with our modelling of thermodynamics.

Both schemes (MOC and INTMOC) are explicit and unconditionally stable which are standard features for numerical methods of characteristics. The handling of boundary conditions is achieved in the present study at order 1. Although the time step can be chosen independently from the mesh size, it must be small enough to provide accurate transient results.

Numerical examples

In this section, we focus on data sets for which phase transition occurs (for the simulations of single-phase ows, we refer to [START_REF] Bernard | Study of low Mach nuclear core model for single-phase ow[END_REF]).

To simulate some scenarii in a PWR, parameters are set as follows:

Geometry of the reactor: L = 4.2 m;

Discretization parameters: N = 100 mesh nodes (∆y ≈ 4.3 cm) and ∆t = 0.01 s;

Reference value for the pressure, the power density and the velocity: p 0 = 155 × 10 5 Pa, Φ 0 = 170 × 10 6 W • m -3 , ṽ = 0.5 m • s -1 ;

Initial data: h 0 (y) = h e , v 0 (y) = v e + y 0 β(h 0 (z))Φ(0, z)/p 0 dz (well-prepared insofar as they satisfy the divergence constraint (1.1a));

Constant viscosity: µ 0 = 8.4

• 10 -5 kg • m -1 • s -1 .
To perform simulations, we have to provide physical values for liquid and vapor stiened gas parameters. This is achieved by using the tables of water and steam [START_REF] Lemmon | Thermophysical Properties of Fluid Systems[END_REF]

Two-phase ow with phase transition

In the rst test, we investigate the ability of our model to deal with two-phase ows with phase transition. We consider the case when the inlet density ρ e , the inlet velocity v e , the initial condition h 0 and the power density Φ are constant, so that we can apply Proposition 3.4 to compute exact transient and asymptotic solutions. The boundary conditions are ρ e (t) = 750 kg • m -3 and v e (t) = ṽ, thus h 0 (y) = h e (t) = h (ρ e ) ≈ 1.190 × 10 6 J • K -1 . The power density is set constant in space and time and equal to Φ 0 . With these parameters, the domain is initially lled with liquid. Then according to Proposition 3.4 mixture appears at time t s 1.769 s for y > y s 0.964 m and pure vapor appears at time t s g 2.929 s for y > y s g 4.002 m. The asymptotic state is reached at t g (L) 2.957 s.

Figure 8 displays numerical results for the enthalpy and the velocity at instants t = 2.1 s, t = 2.8 s and t = 3.5 s. At this last time the solution has already reached the asymptotic regime. Figure 9 displays the mass fraction and the Mach number computed from the enthalpy and the velocity. We mention that the foot of the characteristic curve ξ n i is always computed by means of the 2nd-order expression (4.6) and the velocity eld by means of formulae (4.10) and (4.10').

We observe that all numerical results match the behaviour of the exact solution including phase transition. Nevertheless as it can be noticed at time t = 2.8 s, there is some discrepancy close to the singularity (interface between information coming from the boundary and the inner domain). This can easily be accounted for. Indeed, the scheme relies on an interpolation process which assumes some smoothness for the solution. As it is not smooth at the singularity, the scheme provides a smoothened version of the exact solution. It must be underlined that the higher the degree of interpolation, the more accurate the solution. Moreover, there is a clear advantage of INTMOC versions of the scheme over MOC which is due to the better approximation of the right hand side in the former case. For these reasons, from now on, we will use only the INTMOC_2 scheme which is more accurate.

Figure 9 shows that the Mach number remains lower than 0.05 so that the low Mach number hypothesis of the model is valid in this conguration. Observe that during the transient state, the Mach number is greater than at the asymptotic state.

From a physical point of view, it is worth emphasizing that our model predicts the appearance of some vapor at the top of the reactor for this data set. This is due to a too small inow velocity.

Non monotone Φ

We perform a simulation for a steady power density piecewise constant in space

Φ(t, y) = Φ 0 , if y ≤ L/2, 0, if y > L/2.
This power density models the situation where the control rods are blocked at the middle of the reactor. The boundary conditions are ρ e (t) = 750 kg • m -3 and v e (t) = ṽ, thus h 0 (y) = h e (t) = h (ρ e ) ≈ 1.190 × 10 6 J • K -1 . At the initial state the core thus is lled with liquid, so that the void fraction is 0.

We emphasize that due to the steady boundary condition for the pressure and to the fact that the dynamic pressure is decoupled from the other variables in 1D, the right part [L/2, L] of the domain does not inuence the left part [0, L/2]. This can also be seen in the statement of Proposition 3.4 where the process (for h and v) evolves from y = 0 to y = L. Hence, we can apply Proposition 3.4 over [0, L/2] since Φ is constant within. However, given the solution over [0, L/2] (and more specically at y = L/2), we cannot apply Proposition 3.4 over [L/2, L] with these inlet boundary conditions at y = L/2 as they are unsteady until the asymptotic state is reached. Besides, we know the global asymptotic solution which is given by Proposition 3.6.

We plot on Figure 11 the asymptotic solution from Proposition 3.6 and the numerical one computed by means of the INTMOC_2 scheme for the enthalpy, the mass fraction, the dynamical pressure, the density and the temperature at the instants t = 0.0 s, t = 1.0 s, t = 2.0 s, t = 3.0 s and t = 4.0 s. The velocity is given at instants t < t s , t = 1.8 s and t = 2.0 s. At t = t s 1.768 s mixture appears over [y s , L/2] where y s 0.964 m. This mixture is transported to ll the domain y ≥ y s . As long as the domain is lled with liquid, the velocity does not depend on time. When mixture appears, β is piecewise constant in the domain so that, according to Proposition 3.4, the velocity increases up to the asymptotic value. The dynamical pressure is computed by scheme (4.12) and decreases towards the asymptotic state.

A simplied scenario for a Loss of Flow Accident

Our model is tested here on an accidental transient regime: a main coolant pump trip which is a Loss Of Flow Accident (LOFA) as at the beginning of the Fukushima accident in reactors 1, 2 and 3. To simulate this scenario, the domain is lled at rst with liquid water and we have a normal behaviour of the reactor: the pumps work normally and control rods are in upper position. We impose at the entrance the velocity v e = 10ṽ and the power density is equal to Φ 0 .

These constants are chosen so as to prevent the appearance of mixture. The asymptotic state is reached at t 0.8 s.

At time t 1 = 1.5 s most of the pumps stop, so that the entrance velocity decreases. At time t 2 the security system makes control rods drop into the core decreasing abruptly the power density.

However, there is still some residual power density (here set to 7%Φ 0 ). This instant t 2 2.85 s corresponds to the reaction time (0.3 s) of the safety system (drop of the control rods) after mixture appears (t 2.55 s) and is detected. At time t 3 the security pumps are turned on, thus the inow is re-established. We then compute the solution until the asymptotic state is reached.

Functions v e (t) and Φ(t) can be modelled as follows (see Figure 12):

v e (t) =      10ṽ, if 0 ≤ t < t 1 , 0.2ṽ, if t 1 ≤ t < t 3 , 10ṽ, if t ≥ t 3 , Φ(t) = Φ 0 , if 0 ≤ t < t 2 , 7%Φ 0 , if t ≥ t 2 .
In Figure 13 we report the behaviour of the mass fraction and temperature computed at dierent instants. Thanks to our model we can predict the appearance of some steam in the core depending of the value of t 3 : Case A: for t 3 = 40 s, the asymptotic state corresponding to Φ = 7%Φ 0 and v e = 0.2ṽ is established. In this case, due to the residual power density the uid is completely vaporized at the top of the domain during the transition (see for example Figure 13a at time t = 30 s) even though there are only liquid and mixture phases in the corresponding asymptotic state.

Case B: for t 3 = 20 s, the pumps are re-started soon enough so that the appearance of pure steam is avoided.

Case C: for t 3 = 4 s, even though the pumps are re-started almost instantaneously, some mixture appears in a large part of the domain.

In all cases when the pumps are re-started, the uid comes back to the liquid phase. As expected, the sooner the pumps are re-started, the safer the situation.

Conclusion & Perspectives

We proposed in this paper simulations of a low Mach number model named Lmnc for uid ows in nuclear reactor cores coupled to an adaptive stiened gas equation of state (EOS) which varies according to the phase of the coolant uid (which can be pure liquid water, pure steam or mixture of these two phases). The monodimensional (1D) numerical strategy we worked out leads to accurate and relevant qualitative results in accordance with what was expected. These results were obtained with an optimal computational cost since the velocity and the dynamic pressure are directly integrated in the 1D case.

More precisely, compared to a previous study [START_REF] Bernard | Study of low Mach nuclear core model for single-phase ow[END_REF], the present work enables to deal with more realistic situations insofar as the model allows for phase transition and, thus, for accidental scenarii such as a Loss of Flow Accident (LOFA) induced by a coolant pump trip event. Nevertheless, the method used to determine parameters involved in the stiened gas EOS seems to be restrictive given the range of temperature that must be considered. That is why another strategy will be investigated in [START_REF] Dellacherie | Study of low Mach nuclear core model for two-phase ows with phase transition II: tabulated EOS[END_REF] using tabulated laws instead of the stiened gas EOS.

From a numerical point of view, a variant of the classic method of characteristics was proposed.

It relies on the EOS and provided second-order accurate results. This algorithm was compared to the unsteady analytical solutions derived in this study for dierent data sets including congurations where phase transition occurs.

Moreover, we recall that the thermodynamic pressure is constant in the Lmnc model. Hence, whatever the EOS used to model each phase and no matter what the dimension of the domain (1D, 2D or 3D), the computational cost related to the EOS will be much lower in the discretization of the Lmnc model than in any strategy for the compressible model.

The adaptation of the proposed algorithm to dimensions 2 and 3 will be carried out in future works. However, it must be underlined that although the Lmnc model may be a useful tool for safety studies, it is not designed to handle all potential situations especially when the Mach number cannot be considered small anymore. This is why it is important to study the possibility to couple the Lmnc model to the compressible system from which it is derived (this approach will be investigated in future works). Likewise, couplings with systems dedicated to other circuits in the reactor must be considered.

A. Derivation of the mixture density formula

The mechanic and thermodynamic equilibria hypothesis in the mixture at saturation implies that the pressures are identical in each phase. This allows to obtain (2.2) by using (2.1b). This results in the equality for the mixture density ρ m ρ m = αρ s g h s g + (1 -α)ρ s h s h = αρ s g + (1 -α)ρ s where ρ s κ and h s κ are the density and the enthalpy of the phase κ ∈ { , g} at saturation, i.e.

ρ s κ = ρ κ p, T s (p) and h s κ = h κ p, T s (p) where T s (p) is the temperature T solution of g (p, T ) = g g (p, T )). We deduce α is prescribed by α(h, p) = ρ s h -h s (ρ s g h s g -ρ s h s ) -h • (ρ s g -ρ s )

.

Hence we obtain ρ m (h, p) = ρ s + α(h, p) • (ρ s g -ρ s ) = ρ s g ρ s (h s g -h s ) ρ s g h s g -ρ s h s -h • (ρ s g -ρ s ) , which gives the expression (2.3) of ρ m as a function of (h, p).

B. Computation of the speed of sound for two-phase media

We detail in this section how to express the speed of sound in the diphasic case. This enables in particular to compute the Mach number in any case and to highlight that it remains small.

B.1. General case

In [START_REF] Allaire | A strictly hyperbolic equilibrium phase transition model[END_REF][START_REF] Faccanoni | Étude d'un modèle n de changement de phase liquide-vapeur[END_REF][START_REF] Faccanoni | Modelling and simulation of liquid-vapor phase transition in compressible ows based on thermodynamical equilibrium[END_REF], it has been proved that the speed of sound of the mixture is always positive. To compute explicitly the speed of sound as a function of the enthalpy and the pressure, we start from the usual thermodynamic relation T dη = dε -p ρ 2 dρ.

On the one hand, the relation h = ε + p ρ coupled to the previous statement leads to d(ρε) = (ρT )dη + hdρ. In the mixture, using (2.2) and noticing that ρ κ h κ = ρ s κ h s κ (since T = T s (p) in the mixture), which only depends on p, we can write The speed of sound is discontinuous across the saturation curve, and it is smaller in the mixture than in any pure phase. The minimum value is reached for void ratios close to zero (see [START_REF] Faccanoni | Étude d'un modèle n de changement de phase liquide-vapeur[END_REF] and Figure 3d).

B.2. Stiened gas EOS

To compute the speed of sound c κ in a pure phase (κ ∈ { , g}) using (B.2a), we need to express the enthalpy h as a function of the density ρ κ and the pressure p. Inverting (2.12a), we obtain h κ (ρ, p) = q κ + γ κ (p + π κ ) ρ(γ κ -1) so that (B.2a) becomes c κ (ρ, p) = γ κ (p + π κ ) ρ .

In the mixture, the speed of sound c m satises equation (B.2b) which is written as c m (h, p) 2 = h -q m -α(ρ s g ) + (1 -α)(ρ s ) q m + α(ρ s g h s g ) + (1 -α)(ρ s h s ) -1

.

In this expression, we introduced the following notations

C. Approximation of the water EOS with the stiened gas law

In this appendix, we describe the method introduced in [START_REF] Métayer | Elaborating equations of state of a liquid and its vapor for two-phase ow models[END_REF][START_REF] Métayer | Modelling evaporation fronts with reactive Riemann solvers[END_REF][START_REF] Saurel | Modelling phase transition in metastable liquids: application to cavitating and ashing ows[END_REF] to choose parameters describing each pure phase of water with a stiened gas law to best t with saturated curves (see Figure 2a).

Then, we determine relevant values for liquid water and steam.

More precisely, we suppose that each phase is described by its own stiened gas EOS (2.7), so that we have to compute ve constants for each phase: the specic heat at constant volume c vκ , the adiabatic coecient γ κ , the reference pressure π κ , the binding energy q κ and the reference entropy q κ dened by (2.9). We recall that κ = refers to the liquid phase and κ = g to the vapor one. These parameters will be determined using experimental values at saturation.

The stiened gas law consists of the following expressions for the density, the enthalpy and the Gibbs potential as functions of the temperature T and the pressure p: ρ κ (T, p) = p + π κ (γ κ -1)c vκ T ,

(C.1a)
h κ (T, p) = q κ + γ κ c vκ T, (independent of p)

(C.1b)
g κ (T, p) = q κ + T c vκ γ κ -q κ -c vκ γ κ ln T + c vκ (γ κ -1) ln(p + π κ ) .

(C.1c)

For a given temperature T , the experimental data needed for the computation are the pressure at saturation p s exp (T ), the enthalpies at saturation h s κ,exp (T ) and the densities at saturation ρ s κ,exp (T ). We can nd these data in [START_REF] Lemmon | Thermophysical Properties of Fluid Systems[END_REF] for several uids. Let T 0 and T 1 two reference temperatures (recorded in [START_REF] Lemmon | Thermophysical Properties of Fluid Systems[END_REF]) between the triple point and the critical point (see Figure 2b).

•

Step I: computation of γ κ c vκ . From the analytical expression of the enthalpy (C.1b), we have h κ (T ) = γ κ c vκ . This yields an averaged value of the product c pκ def = γ κ c vκ by means of a linear approximation of the experimental enthalpies between reference states T 0 and T 1 as c pκ = h s κ,exp (T 1 ) -h s κ,exp (T 0 ) T 1 -T 0 .

(C.2)

•

Step II: computation of q κ . Equation (C.1b) applied to the reference state T 0 provides the approximation of the binding energy q κ : q κ = h s κ,exp (T 0 ) -c pκ T 0 .

•

Step III: computation of π κ . At saturation, there is an algebraic relation between the temperature and the pressure, e.g. we can write p = p s (T ). Thus, using (C.1a), the specic densities at saturation are expressed as By convention, we take q = 0 J • K -1 and determine the coecient q g using experimental values of the pressure at saturation q g = A ln p s exp (T 0 ) + π g -B ln p s exp (T 0 ) + π -C(ln T 0 -1) + D T 0 .

(

  which is the case in the Lmnc model) and when the phase change is modelled by assuming local .: 330 • C) Water coolant (input temp.: 290 • C) (a) Scheme of the primary circuit of a PWR. .: 290 • C) Water coolant (input temp.: 280 • C) Liquid Vapor (b) Scheme of the primary circuit of a BWR.

Figure 1 :

 1 Figure 1: Scheme of nuclear reactors whose coolant is water: the major dierence between PWR

Figure 2 :

 2 Figure 2: Saturation and coexistence curves.

  Density as a function of the enthalpy: the densities of the liquid phase and the vapor phase at saturation are ρ s ≈ 632.663 kg • m -3 and ρ s g ≈ 52.937 kg • m -3 . Temperature as a function of the enthalpy: the temperature at saturation is T s ≈ 654 K. Compressibility coecient as a function of the enthalpy: the compressibility coecients of the liquid phase, the mixture at saturation and the vapor phase are β ≈ 0.008768, βm ≈ 0.194852 and βg ≈ 0.300699.

  c m (h s g ) c (h s ) h s c m (h s ) (d)Speed of sound as a function of the enthalpy: the speed of sound of the liquid phase and the vapor phase at saturation are c (h s ) ≈ 1942 m • s -1 and cg(h s g ) ≈ 647 m • s -1 ; in the mixture at saturation the speed of sound is in the range (187 m • s -1 , 579 m • s -1 ).

Figure 3 :

 3 Figure 3: EOS with phase transition for parameters of Table 1 at page 49 with p 0 = 155 × 10 5 Pa:the enthalpies of the liquid and vapor phases at saturation are h s ≈ 1.627 × 10 6 J • K -1 and h s g ≈ 3.004 × 10 6 J • K -1 .

Proposition 3 . 2 .

 32 Let us assume that the power density Φ depends only on space; h e and D e satisfy Hypotheses 1.2 and 3.1 and are such that v e = D e /ρ(h e , p 0 ) is independent of time;

  Liquid (L), mixture (M) and vapor (G) regions of the spatiotemporal domain R + × R + for the velocity. Denition of four regions of the spatiotemporal domain R + × R + for the enthalpy.

Figure 5 :

 5 Figure 5: Denition of regions for Proposition 3.4

Figure 6 :

 6 Figure 6: Grid and boundary conditions.

Figure 7 :

 7 Figure 7: Numerical method of characteristics.

1

 1

(4. 12 ) 4 . 3 .

 1243 Positivity-preserving property We have the following result which is the discrete version of Lemma 3.1: Proposition 4.1. When Hypothesis 3.1 is satised, strategies (4.2)& (4.9a) and (4.4)& (4.9b)

Figure 10

 10 Figure 10 displays the density and the temperature computed from the enthalpy. In those gures we compare the exact and asymptotic solutions to the dierent versions of the numerical scheme, namely: by scheme (4.2)&(4.9a) with linear interpolation (4.7) for hn+1,n i (called MOC_1), by scheme (4.2)&(4.9a) with high order interpolation (4.8) for hn+1,n i (called MOC_2), by scheme (4.4)&(4.9b) with linear interpolation (4.7) for hn+1,n i (called INTMOC_1), by scheme (4.4)&(4.9b) with high order interpolation (4.8) for hn+1,n i

Figure 8 :

 8 Figure 8: Numerical enthalpy (left) and velocity (right) for test 5.1: Two-phase ow with phase transition. We compare the four numerical solutions to the analytical and asymptotic ones. The horizontal dotted lines correspond to = h s and h = h s g .

Figure 11 :

 11 Figure 11: Numerical results of the test 5.2: Non monotone Φ.

  t

Figure 12 :

 12 Figure 12: Power density and entrance velocity for the test 5.3: A simplied scenario for a Loss of Flow Accident.

Figure 13 :

 13 Figure 13: Numerical mass fraction (left) and temperature (right) of the test 5.3: A simplied scenario for a Loss of Flow Accident.

1 ,(B. 1 )

 11 since ρε = ρh -p.In the pure phase κ ∈ { , g}, the volume fraction α is 0 (i.e. κ = ) or 1 (i.e. κ = g), so that ρ = ρ κ , h = h κ and equation (B.

(B. 3 )

 3 Remark B.1.

  -1)c vκ T s (p) 1 -(p + π κ ) (T s ) T s (p) , -1)c vκ T s (p) q κ + γ κ c vκ T s (p) -(p + π κ )q κ (T s ) T s (p) , (T s ) (p) = (γ g -1)c vg p + π g -(γ -1)c v p + π (c vg γ g -c v γ ) + q g -q T s (p) ) -h s (p) T s (p) = T s (p) p β m (p).

••

  ρ s κ (T ) = ρ κ T, p s (T ) = p s (T ) + π κ (c pκ -c vκ )T , (C.3)which evaluated at T 0 and T 1 yieldsρ s κ (T 1 ) ρ s κ (T 0 ) = (p s (T 1 ) + π κ ) T 0 (p s (T 0 ) + π κ ) T 1 ,so that we can determine an averaged value of the coecient π κ using experimental values of the density and the pressure at saturationπ κ = T 0 ρ s κ,exp (T 0 )p s exp (T 1 ) -T 1 ρ s κ,exp (T 1 )p s exp (T 0 ) T 1 ρ s κ,exp (T 1 ) -T 0 ρ s κ,exp (T 0 ).Step IV: computation of c vκ . Equation (C.3) applied to the reference state T 0 provides the approximation of c vκc vκ = c pκ -p s exp (T 0 ) + π κ T 0 ρ s κ,exp (T 0 ).Step V: computation of γ κ . Using (C.2) we deduce the approximation of γ κ γ κ = c pκ c vκ .•Step VI: computation of q κ . At thermodynamic equilibrium, the two Gibbs potentials are equal. Using (C.1c), this impliesA ln p s (T ) + π g -B ln p s (T ) + π -C(ln T -1) + D T + q -q g = 0 with A def = c pg -c vg , B def = c p -c v , C def = c pg -c p , D def = q g -q .

  Proposition 3.1. Let us assume that the power density Φ = Φ 0 > 0 is constant in time and space; h e and D e satisfy Hypotheses 1.2 and 3.1 and are such that v e = D

	ξ(t, y) = y +	v e Φ0	e -Φ0t -	v e Φ0	,
	t * (t, y) = t -	1 Φ0	ln 1 +	Φ0 v e	y = -	1 Φ0	ln 1 +	Φ0 v e	ξ(t, y) ,

e /ρ(h e , p 0 ) is independent of time; h 0 veries Hypotheses 1.4 and 3.1. Let us denote Φ0 def = β κ Φ 0 /p 0 . Then ξ(t, y) and t * (t, y) dened by (3.2) are equal to and the solution h of equation (1.1b) supplemented with BC (1.3a) is given by

  Corollary 3.1. Under the hypotheses of Prop. 3.1, we have y * (t) = (e Φ0t -1) ve Φ0 and the solution p of equation (1.1c) together with BC (1.3c) is given by if y

	(3.4) ; (3.5b) 0 e Φ0t H 0 ξ(t, L) -H 0 (0) g + Φ0 v e e Φ0t H 0 ξ(t, L) -H 0 ξ(t, y) p(t, y) = Φ0 [µ(y) -µ(L)]+ v e p 0 β κ p 0 β κ + Φ2 0 e Φ0t H 0 ξ(t, L) -H 0 ξ(t, y) otherwise p(t, y) = Φ0 [µ(y) -µ(L)]+ g H e t p(t, y) = Φ0 [µ(y) -µ(L)]+ p 0 β κ g + Φ0 v e e Φ0t H 0 ξ(t, L) -H 0 (0) + Φ2

e t * (t, y) + Φ 0 y D e t * (t, y) , if ξ(t, y) < 0. * (t) > L, then * (t, y) -H e t * (t, L) + Φ0 v e e Φ0t H e t * (t, y) -H e t * (t, L) ; (3.5a) if y * (t) ≤ L and y ≥ y * (t) then + v e g H e t * (t, y) -H e (0) + Φ0 v e e Φ0t H e t * (t, y) -H e (0) ; (3.5c)

  Then the solution h of equation (1.1b) with BC (1.3a) is given by

	Proposition 3.3. Let us assume that
	the power density Φ depends only on time;
	h e and D e satisfy Hypotheses 1.2 and 3.1;
	h 0 veries Hypotheses 1.4 and 3.1.
	Let us dene	Ψ(t) def =	β κ p 0	0	t	Φ(s) ds.
	We thus have

ξ(t, y) = ye -Ψ(t) -t 0 v e (s)e -Ψ(s) ds and t * (t, y) is the solution of the equation (upon t) y = t t * v e (s)e Ψ(t)-Ψ(s) ds.

  as described in the appendix C. Parameters involved in the stiened gas law for the following simulations are those in Table1, page 49. For this data set, variables at saturation are h s ≈ 1.627 × 10 6 J • K -1 , h s g ≈ 3.004 × 10 6 J • K -1 , ρ s ≈ 632.663 kg • m -3 and ρ s g ≈ 52.937 kg • m -3 .

	In Test 5.1, a two-phase ow with phase transition is considered and enables to compare the
	numerical schemes presented in Section 4, while all other tests only involve scheme (4.4)&(4.9b)
	with high order interpolation (4.8) (noted INTMOC_2) to assess the robustness of the scheme and
	the relevance of the model.

This is not the case for the 1D compressible system from which the Lmnc is derived.

The existence of fast acoustic waves in the compressible system induces numerical diculties see[START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number[END_REF][START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit[END_REF] for example which cannot arise in the Lmnc model since the acoustic waves have been ltered out to obtain this low Mach number model.

 1: Liquid water and steam, parameters computed in [START_REF] Métayer | Elaborating equations of state of a liquid and its vapor for two-phase ow models[END_REF][START_REF] Métayer | Modelling evaporation fronts with reactive Riemann solvers[END_REF].

Reference states T 0 and T 1 must now be specied. First, T 0 is a reference state chosen in order to have the best t between the theoretical pressure p s (T ) and the experimental one p s exp (T ).

Secondly, this strategy for the determination of parameters for a stiened gas EOS is accurate if the two reference states are suciently close. Practically (following [START_REF] Métayer | Elaborating equations of state of a liquid and its vapor for two-phase ow models[END_REF][START_REF] Métayer | Modelling evaporation fronts with reactive Riemann solvers[END_REF]), we took for the liquid phase T 0 = 298 K and T 1 = 473 K (steps I and II) and T 0 = 439 K and T 1 = 588 K (steps III to VI). As for the steam, we set T 0 = 298 K and T 1 = 473 K in any step.

The values for liquid water and steam are given in Table 1. As underlined in [START_REF] Faccanoni | Étude d'un modèle n de changement de phase liquide-vapeur[END_REF][START_REF] Métayer | Elaborating equations of state of a liquid and its vapor for two-phase ow models[END_REF][START_REF] Métayer | Modelling evaporation fronts with reactive Riemann solvers[END_REF], these parameter values yield reasonable approximations over a temperature range from 298 K to 473 K.

Nevertheless near the critical point, there are some restrictions due to the nonlinearity of the enthalpy with respect to the temperature. To circumvent this limitation, we shall propose in [START_REF] Dellacherie | Study of low Mach nuclear core model for two-phase ows with phase transition II: tabulated EOS[END_REF] a new strategy to better approximate the thermodynamic quantities involved in the Lmnc model.