A homogeneous relaxation low Mach number model

. In the present paper, we investigate a new homogeneous relaxation model describing the behaviour of a two-phase ﬂuid ﬂow in a low Mach number regime, which can be obtained as a low Mach number approximation of the well-known HRM. For this speciﬁc model, we derive an equation of state to describe the thermodynamics of the two-phase ﬂuid. We prove some theoretical properties satisﬁed by the solutions of the model

A hierarchy of relaxation two-phase flow models can be found in [14,26,42,43]. In this hierarchy, two families can be considered: the inhomogeneous flows [42], with different velocities for the two phases formulated using two momentum equations and velocity relaxation; and the homogeneous models [43], where there is no relative velocity between the two phases.

Homogeneous models
In homogeneous models, a single velocity is considered to describe the flow, and the mixture is treated as a single fluid. These models require three partial differential equations which govern the evolution of the total mass, the global momentum and the total energy of the whole mixture and some transport equations as well as a number of externally supplied relations to specify the interaction between the two phases.
The most general homogeneous model is the 6-equation model presented for instance in [22,31,32,46,48] where the flow is described using three fractions (mass, volume and energy) and the phase interactions are accounted by source terms of the form ( eq − )/ , where eq is the fraction at equilibrium and the exchange between the two phases occurs with a given characteristic time .
A homogeneous model is relaxed when at least one equilibrium between phases is assumed. Relaxed models will be denoted HRM for Homogeneous Relaxation Model. Among these models, we can consider the models where only the equality of pressures is taken into account [1], or only the equality of temperatures [1,36], or equality of both pressures and temperatures [26,43], or in the case of liquid and vapour water, equality of pressures and the saturation of the water vapour phase [4,12,15,20,23].
The Homogeneous Equilibrium Model (acronym HEM) is the simplest of the homogeneous models. It assumes that the two phases are at thermodynamic equilibrium. In this case, a set of three equations in order to account for total mass, total momentum and total energy balances is sufficient to describe the flow.
In this paper we focus on two-phase flows (liquid and vapour phases of the same fluid), described by a fourequation model with equality of both pressures and temperatures of the phases and supplemented with a source term relaxing the mass fraction to an equilibrium (saturation) mass fraction.

Low Mach number models
In some applications, like in nuclear cores, convection is characterised by a low Mach number flow, where the convective velocities are much slower than the speed of sound in the fluid, typically by one to two orders of magnitude. In fact, when dealing with homogeneous models for liquid-gas flows, the speed of sound in the mixture can be much lower than the one in the pure phases, leading to a potentially "high" Mach number, even when the convective velocity of the flow is "small". Nevertheless, this effect is not so significant in the high pressure/high temperature regime, making low Mach number models legitimate in nuclear core applications. In practical terms, the speed of sound in the HEM is always smaller than the speed of sound in the HRM [42], leading to a higher Mach number. Since the low Mach number hypothesis has already been checked for the HEM in the nuclear core setting [8], it is also valid for the HRM.
The disparity between the time scales of convective motions and sound waves is a major computational challenge when convection is the phenomenon of interest [47]. Following sound waves explicitly, as in a traditional compressible approach, introduces a much shorter time scale than the one driven by the convective motions in the computation, making it difficult to simulate expected time scales, which are long with respect to the convective time scale. This made the development of so-called low Mach number models attractive, as sound waves are filtered out (e.g. [2,3,21,44]).
Low Mach number formulations replace the compressible flow equations by a constrained system of partial differential equations with a similar structure to the incompressible Navier-Stokes equations. The equations of hydrodynamics are reformulated in order to analytically remove acoustic waves while keeping local compressibility effects due to heat release and phase changes. Because low Mach number models do not track the propagation of acoustic waves, one can use a time step based on the fluid velocity rather than the speed of sound and thus often gain an order of magnitude or more in computational efficiency over a traditional compressible approach.
The fundamental approximation made in the the low Mach number equations is that the compressible pressure can be approximated by a reference pressure in the equation of state. Mathematically, this leads to the addition of a constraint on the velocity field compared to standard hyperbolic evolution equations. The pressure is thus decomposed as ( , ) = * ( ) + ( , ), where * is the reference state pressure (or "background" or "ambient thermodynamic" pressure) and is the perturbed pressure (often called "dynamic pressure"). For low Mach number flows, an asymptotic analysis shows that ( , )/ ( , ) = (ℳ 2 ), where ℳ is the Mach number. We can then approximate ( , ) by * ( ) in the computation of thermodynamics quantities.
A hierarchy of low Mach number models for nuclear reactors In the context of pressurised water reactor cores, an asymptotic low Mach number model for the HEM system has been derived and investigated in a series of papers [8,16,19]. The fluid is described by a single equation of state taking into account the phase transition by supposing that, when both vapour and liquid phases are present, they have the same pressure, temperature and chemical potential. It incorporates large compressibility effects due to the vaporisation and thermal processes with a spatially constant background pressure * = 155 bar.
In the present paper we are interested in studying an asymptotic low Mach number model for a 4-equation HRM: we assume that both the vapour and liquid phases have the same pressure and temperature but different chemical potentials. In the following the low Mach asymptotic expansion of the HEM described in [8] will be referred to as 3-Lmnc model and the low Mach asymptotic expansion of the 4-equation HRM as 4-Lmnc model, where Lmnc stands for Low Mach Nuclear Core.

Content of the paper
The organisation of the paper is as follows. We first present the 4-Lmnc model in Section 2: the system of equations with boundaries and initial conditions. In Section 3, we derive the isothermal isobar equation of state to close the system based on a Noble-Abel stiffened gas law for each phase. Then, we describe the relaxation term in Section 4.1 and the different regimes that can be considered. Section 4 is devoted to the investigation of the non-instantaneous relaxation regime: we prove some properties of the system (maximum principle, positivity of the source terms, analytical steady solution) and we introduce a well-balanced numerical scheme mimicking these theoretical properties.
In Section 5, we study the instantaneous relaxation regime towards the 3-Lmnc model: we first recall the equations governing the 3-Lmnc model and its closure law (isothermal, isobar, iso-chemical potential equation of state). We then prove the formal convergence of the 4-Lmnc model towards the 3-Lmnc model. An asymptoticpreserving numerical scheme is then introduced to handle numerical simulations in the stiff relaxation regime.
Let us notice that for the sake of simplicity, the analysis of the present paper is restricted to dimension 1.

Governing equations
As for the 3-Lmnc model [16], starting from the 4-equation HRM, a standard asymptotic analysis around the reference pressure * yields the following system of conservative equations (called 4-lmnc model) where is the total specific density, the velocity field, ℎ the total specific enthalpy. The power density Φ( , ) ≥ 0 models the heating due to fission reactions, and might be varying in time and space. ∈ [0; 1] is the mass fraction: if = 1 the fluid is in vapour phase, if = 0 is in liquid phase and when 0 < < 1 the fluid is a liquid-vapour mixture. Finally ℛ is a relaxation term accounting for interactions between phases.
To close the system, we have to provide a closure relation between the thermodynamic variables , ℎ and , called "equation of state", modelling the thermodynamic properties of the fluid. The fluid can be a pure phase (liquid or vapour) of the same fluid (e.g. water and steam) or a mixture of both phases.
The fundamental change with respect to the fully compressible model (the HRM model) lies in the replacement of the full pressure field in the equation of state by the constant reference pressure * > 0 ( * = 155 bar for a Pressurised Water Reactor (PWR)), and in the momentum equation by the dynamic pressure .
Under smoothness assumptions we can derive a non-conservative formulation equivalent to (2.1), which is the one we shall focus on for analysis and derivation of numerical schemes.
Let us define = 1/ the specific volume. System (2.1) can be written in a nonconservative form Observe that in dimension 1, this equation allows to integrate directly the velocity (with boundary conditions), whereas in higher dimensions, this is a divergence constraint (as in the Navier-Stokes system for example). The (unknown) dynamic pressure is only involved in the fourth equation of System (2.2), and can be computed as a post-processing after the computation of ℎ, and . Because of this decoupling, this last equation will be left apart in the following. Finally, we focus on the following system

Boundary and initial conditions
The model is set in some bounded domain Ω = (0, ), which may represent the nuclear core. The fluid is injected at the bottom with a given enthalpy ℎ > 0, mass fraction ∈ [0; 1] and at a given flow rate ( ) > 0: we assume the flow to be upward (which corresponds to a nuclear power plant of PWR or BWR type). The boundary conditions are thus written The system is supplemented with initial conditions:

Closure relation and relaxation source term
As already mentioned, the system requires appropriate equations of state (referred to as the EoS in the following) in order to describe the "pure vapour" phase, the "pure liquid" phase but also the "mixture" phase. They are specified in Section 3, and the relaxation source term ℛ in Section 4.1.

Equation of state for the 4-lmnc model
The equation of state (EoS) corresponds to the modelling of thermodynamic properties of the fluid and consists of an algebraic relation between thermodynamic variables.
In classic thermodynamics, the thermodynamic state of a pure single-phase fluid is represented by means of a relation between the internal energy , the specific volume and the entropy (see e.g. [13,29]). As a preliminary study, we chose a simple analytical form capable of capturing the essential physics of a pure phase, which is the Noble-Abel stiffened gas equation of state (NASG EoS).

Equation of state for a pure phase
For a given fluid, the Noble-Abel stiffened gas EoS is fully defined by the relation (3.1) The constant Noble-Abel stiffened gas parameters describing thermodynamic properties of the fluid are the following: is the specific heat at constant volume, -> 1 is the adiabatic index, which is a non-dimensional coefficient, -− [Pa] is the minimal admissible pressure, -[m 3 · kg −1 ] is the covolume, the minimal admissible volume, -[J · kg −1 ] is a reference enthalpy, -[J · K −1 · kg −1 ] is a reference entropy (relevant when phase transition is involved).
The Noble-Abel stiffened gas EoS, proposed in [40], is an extension of the classic stiffened gas EoS [38], which is recovered when = 0. The ideal gas EoS is recovered when = = = 0.
Thanks to the Gibbs relation d = d + d , the classic definitions in thermodynamics provide the following expressions of the temperature , the pressure , the enthalpy ℎ and the Gibbs potential as functions of the specific volume and the internal energy : The definition of the entropy requires − − ( − ) > 0 and − >0, which is equivalent to + > 0 (we refer to [49] for a more in-depth discussion on the physical basis for this EoS). Making a change of thermodynamic variables from ( , ) to ( , ), which can be made explicit 1 for this kind of EoS, we obtain the monophasic stiffened gas law: If we denote then we can express ℎ and as functions of ( , ) Making a change of thermodynamic variables from ( , ) to ( , ) we could also write and is given as a function of ℎ, by Notice that, with NASG EoS, the specific heat at constant pressure is always constant while the coefficient depends on . We also note that, with the stiffened gas EoS (i.e. when = 0), the enthalpy depends only on the temperature, but not on the pressure.
Let us use the index =l for the liquid phase or =g for the vapour phase. If all parameters involved in pure phase equations of state are given (i.e. , , , , , , ), then , ,˜and are deduced. Let us denote by sat ( ) the solution of the equationl( , ) =g( , ) (the so-called temperature at saturation). We can then define ℎ sat ( ) def = ℎ ( sat ( ), ) and sat ( ) def = ( sat ( ), ).
Notice that, at a fixed pressure * , all these quantities are constant and satisfy the equalities

Iso-equation of state for a mixture
We consider each phase ( =g orl) as a compressible fluid characterised by its thermodynamic properties, i.e. each fluid is governed by its own given EoS (see Sect. 3.1). The two-fluid mixture is constructed according to isobar and isothermal assumptions: when fluids coexist (i.e. when 0 < < 1), they have the same pressure and the same temperature 2 , so that we consider the volume , the internal energy and the entropy as functions of and .
The mixture specific volume and the mixture internal energy are defined by where is the mass fraction. Let ℎ be the enthalpy of the phase and ℎ the enthalpy of the mixture. When the pressure is the same in both fluids, recalling that the internal energy is linked to the enthalpy by the relation We now assume that each fluid is described by its own Noble-Abel stiffened gas EoS. Using (3.2), we can write Therefore the temperature verifies we can hence obtain the following expression for the specific volume (ℎ, , ) From now on, let us drop the dependency upon ≡ * . Hence we express, using (3.3), Hence the mixture can be considered as a generalised Noble-Abel stiffened gas in the sense that, at constant pressure, coefficients , and˜only depend on (as in [5,6] for a generalised stiffened gas).

Relaxation term
In the 3-Lmnc model [8], the two fluids are the liquid and the vapour phases of the same component. The mixture is supposed to be at thermodynamic equilibrium and the mass fraction is computed to take into account the phase transition (i.e. an instantaneous mass transfer from one phase to the other). The rate of mass transfer from the liquid phase to the vapour phase is due to their difference of Gibbs potentials.
Here the mixture is supposed to be only at isothermal and isobar equilibrium and mass transfer can be modelled by introducing a relaxation source term ℛ allowing the exchange of mass between the two phases with a given characteristic time .
We choose to model the mass transfer as in [4][5][6][31][32][33] by setting where the coefficient represents the relaxation time and the mass fraction (ℎ) is computed to ensure the saturation of the mixture (equality of the pressure, temperature and Gibbs potential of each phase, see Sect. 5.2 for a complete analysis), that is Three regimes can be considered: -the instantaneous relaxation regime: it corresponds to → 0, and we recover formally the equation of state at saturation and the 3-Lmnc model (which corresponds to the equality of chemical potentials); -the infinite relaxation regime: when → ∞, only the convective part is involved (the thermodynamic is too slow to affect the hydrodynamic motion, no mass transfer between phases); -in between: finite values of > 0 lead to an actual relaxation system (non-instantaneous relaxation).
In this paper we shall study two regimes: -the non-instantaneous relaxation regime ( > 0) in Section 4, -the instantaneous relaxation regime ( → 0) in Section 5.

Properties of the 4-lmnc
Let us first state some properties of model (2.3). The well-posedness of the system is not tackled in this work, and we assume the existence of a solution smooth enough for the following computations to make sense (in particular, the velocity is assumed to belong to 1,∞ (Ω) in order to apply the method of characteristics).
A second result that can be proved about this model is the positivity of the relaxation term ℛ which describes the effect of the relaxation process: there is a discrepancy between the phenomena "ℎ becomes greater than ℎ sat l " and " becomes positive", due to the non-instantaneous relaxation.
Concerning the long-time behaviour, we can obtain an explicit form for the steady-state solution for ℎ and the flow rate (defined in (2.4c)), as it is stated in the following straightforward proposition. This steady-state solution will then be used in the next paragraph, where we derive a well-balanced scheme for the model. . Let us assume that ℎ , and the flow rate defined by (2.4c) are independent from , and that the density function Φ depends only on . Then, for any EoS, model (2.3) admits the following steady state: Proof. The steady-state system in the conservative form is the following From the first equation we deduce that is constant in space, so that ( )( ) = for all ∈ [0; ]. The second equation becomes ℎ = Φ with ℎ( = 0) = ℎ so that ℎ(

A well-balanced numerical scheme
In previous papers [8,19], numerical schemes were designed for the 1D 3-Lmnc model based on the method of characteristics. It was proved to be (temperature) positivity-preserving and second-order in space and time.
In the present paper, there is an additional equation, namely a transport equation for the mass fraction . The previous algorithm could have been used coupled to the method of characteristics applied to the additional equation. However another strategy is proposed in this paper in order to mimic the theoretical property stated in the previous proposition: a well-balanced strategy to recover the asymptotic states (Prop. 4.3). The idea is to ensure the stability of some numerical steady states that are consistent with the continuous steady state in the same spirit as what it is done in the framework of hyperbolic equations [28,30,52]. Such a well-balanced approach is well-suited since, in a nuclear core setting, simulations are often devoted to steady-state or near steady-state flows.
The well-balanced property is based on the mass conservation property for the steady-state flow: where sat is defined by (3.11) and where we used the relation (3.12) for and ′ .
Our well-balanced scheme strongly relies on the discrete analogue of this computation, which is stated in the following Lemma.
we have that, for any , where we dropped the time superscripts for simplicity, since this property is true for any .
Proof. Let us expand the right hand side of the claimed equality, using the definitions of ,˜and : Since˜= (˜g −˜ℓ) +˜l and = (g − ℓ ) +l this implies so that, in the previous equation, the terms involving˜,˜l,˜g vanish, and we obtain From the definition of , let us compute which finishes the proof.
Let us now introduce the numerical scheme (4.6d) Observe that since equation (4.5b) can be stiff for small values of due to the right hand side (4.6c), the source term is discretised implicitly. In fact, this does not induce longer computational costs, since the equation is linear in and can be solved explicitly. The CFL condition is thus only related to the transport equation (4.5a) for ℎ.  Proof. Assume that ( , ℎ , ) = (︀ +1 , ℎ +1 , +1 )︀ . Dropping the time indices, Scheme (4.5) reduces to We first prove (4.7), which is the discrete analogue of the mass conservation property (4.4). This is equivalent to proving ( − −1 ) = ( − −1 ) .
Using successively the equations on , ℎ and of the scheme (4.9), we have We then prove (4.8) by writing ℎ − ℎ −1 ∆ = Φ · Remark 4.6. Let us observe that (4.9b) implies thanks to (4.7), which is consistent with the ODE in Proposition 4.5.
Notice that a similar strategy can be adapted to the 3-Lmnc model for which the steady state is completely known (see Appendix A) without any ODE to solve.

Semi-analytical steady state solution
To assess the well-balanced property of the present scheme, we have to compute a steady state solution. We apply Proposition 4.3 to compute the asymptotic solution: -the enthalpy ℎ ∞ ( ) is given by (4.3); -we compute the mass fraction ∞ by solving the Cauchy problem with an explicit Runge-Kutta scheme at order 6 (with 7 stages) on a very fine grid using 51201 points; -we can then compute the velocity ∞ using the EoS:

Constant Φ test case
In the first test, we investigate the ability of our model to deal with two-phase flows with non-instantaneous phase transition ( = 1.0 × 10 −1 ).
The power density is set constant in space and time and equal to Φ = 170 × 10 6 W · m −3 . The boundary and initial conditions are ℎ 0 ( ) = ℎ ( ) = 0.9ℎ sat l , 0 ( ) = ( ) = 0.4 m · s −1 , so that = 306.2 m · kg · s −1 · m −3 and 0 ( ) = ( ) = 0 . With these parameters, the domain [0; ] with = 4.2 m is initially filled with liquid. With these constant values, we can apply the algorithm presented at Section 4.4 to compute a semi-analytic asymptotic solution. Figure 1 displays numerical results at instants = 0 s (blue) and = 6.57 s (red) and the steady solution (magenta dotted) for the enthalpy ℎ, the mass fraction and the velocity . At this last time the solution has already reached the asymptotic regime. The (green dotted) line represents the solution when the mixture is at instantaneous equilibrium (i.e. = (ℎ)). The computation is performed on a grid with 101 points and the CFL constant is equal to 0.99. At final time we have which means that the scheme is well-balanced even with few points. Moreover, we observe where the enthalpy is equal to ℎ sat l and ℎ sat g respectively. We remark that, for > 2.566 m, the fluid is still a mixture (i.e. 0 < < 1) even if ℎ( ) > ℎ sat g , which enlightens the influence of the time delaying the mass transfert between the phases.
Notice that, although the steady enthalpy is the same for both 4-Lmnc and 3-Lmnc models, (which is reached within the same time), the mass fraction evolution is delayed with respect to the enthalpy (with a time scale of order ). This can induce some difficulties to approach numerically the position of appearance of the pure phase in the 4-Lmnc model, and rounding errors can lead to a large error on the position of the interface. Nevertheless, in the context of the 4-Lmnc model, it is not crucial to identify if the fluid is a pure phase ( = 0 or 1) or if there is a small fraction of the other phase, since the definition of the iso-Tp mixture also describes pure phases continuously (cf. Rem. 3.1). This is a remarkable difference with respect to the closure law of the 3-Lmnc model, and represents a major benefit of the 4-Lmnc model since it is less sensitive to small errors in determining the parameters.

Sinusoidal Φ test case
In the second test, we consider a space-dependent constant in time power density function equal to The boundary and initial conditions, the domain and are the same as for the test of Section 4.5.2. As previously, Figure 2 displays numerical results at instants = 0 s (blue) and = 4.09 s (red) and the steady solution (magenta dotted) for the enthalpy ℎ, the mass fraction and the velocity . At this last time the solution has already reached the asymptotic regime. The (green dotted) line represents the solution when the mixture is at instantaneous equilibrium (i.e. = (ℎ)). The computation is performed on a grid with

The instantaneous relaxation regime
As we already stated, the 4-Lmnc model describes a two-phase flow under the assumption of instantaneous mechanical and thermal equilibrium (but the two phases will in general not be at chemical equilibrium). In this section we wish to derive the 3-Lmnc model [8], where the phase change is instantaneous, as the instantaneous relaxation limit of the above 4-Lmnc model. The relaxation term deals with phase change and forces the phases towards chemical equilibrium.
We first recall the equations governing the 3-Lmnc model and its closure law (isothermal, isobar and isochemical potential equation of state). We then study the relaxation limit of the 4-Lmnc model towards the 3-Lmnc model, and we finally introduce an asymptotic-preserving numerical scheme for the 4-Lmnc model.

Systems of equations
In a non conservative formulation, the 3-lmnc and the 4-lmnc models can be written as follows.

Iso-equation of state for a mixture
We recall that = * , thus we drop the dependence on in the description of the equation of state. We proved in [8,19] that the EoS at saturation reads as follows, given sat , ℎ sat , sat , as described in Section 3.1.
-The specific volume in the 3-Lmnc model is given by -The mass fraction in the 3-Lmnc model is given by (4.2). -Temperature in the 3-Lmnc model is expressed by -The coefficient is not defined in the iso-Tpg equation of state, since which is equal to zero in the mixture at saturation.
To study the relaxation of the 4-Lmnc towards the 3-Lmnc system, we first compare the equations of state of the iso-and the iso-mixtures.
then, by definition (3.11) of sat ( ) Finally, for the temperature, we can compute In pure phases, all quantities ( ,˜, , and ) clearly coincide. Moreover, although phase change is described by = 0 (resp. 1) in the iso-Tp EoS whereas it is described by ℎ = ℎ sat l (resp. ℎ sat g ) in the iso-Tpg EoS, the two descriptions coincide when = (ℎ), and phase change thus occurs at the same point.
Then we consider the following modified system: Let us expand any variable of Model (5.3) as = (0) + (1) + ( 2 ). At order −1 , we get from the form Moreover, at order 0 , we have from (5.3c) and from (5.3b) bounded which leads to the second equation of (5.2) for ℎ (0) up to a correction of order . Further, we compute from the definition of S The assumption on the EoS leads to which is exactly the first equation of (5.2) for (0) up to a correction of order . Last, letting tend to zero concludes the proof.

An asymptotic-preserving numerical scheme for the 4-lmnc model
In this section, an asymptotic-preserving strategy is proposed to be able to simulate the 4-lmnc model in stiff regimes (small values of ) with a reasonable computational cost. The scheme relies on a time splitting between the relaxation step (stiff, thus implicit) and the transport one (explicit). It is inspired by the standard asymptotic-preserving approach introduced by Jin [24,34,35].
Let us consider the following semi-discrete (in time) scheme is an approximation of and (·) is any discretisation of (·). As for the previous scheme, since the second equation of (5.6a) can be stiff for small values of , the source term is discretised implicitly. Again, this does not induce longer computational costs, since the equation is linear in and can be solved explicitly. The CFL condition is thus only related to the transport equation on the enthalpy at speed (two last equations of (5.6a)).
Proposition 5.5. Under the previous assumptions, Scheme (5.6) is weakly asymptotic-preserving in the mixture, in the sense that its numerical solution is an ( )-approximation of the numerical solution to the scheme in the mixture.
Proof. Let us denote where we used a Taylor expansion of (which is a C 1 function with respect to ℎ and ), and the constants 1 and 2 are given by Then, due to the first equation of (5.6a). Hence Moreover, we have from the first equation of (5.6a) Hence, using the second equation of (5.6a) . Therefore thanks to equation (5.6b) and the third equation of (5.6a) becomes together with previous equalities Applying a Taylor expansion to the function with respect to , we obtain Hence Recalling that (ℎ) = (ℎ − )/ and ( +1 ) = Φ/ , we finally obtain where we used the first equation of (5.6a). Then, the second part of the splitting is computed as Likewise, in the mixture, we have that ) due to the properties of the operator .

Asymptotic-preserving behaviour of the scheme
In this section, we illustrate numerically the good behaviour of the numerical scheme (5.6). Let us set the discretisation parameter ∆ = 2.1 × 10 −2 . The CFL condition is equal to 0.99. Since the velocity remains less than 10 in this test case, the time step ∆ is always greater than 2.1 × 10 −3 , independently of the value of .
We compare the solutions computed with the 4-lmnc and the 3-Lmnc models for different values of the relaxation parameter , from 10 −2 to 10 −7 . Figure 4 displays the 2 norm of the relative error between these solutions in semi-log-y scale for ∈ [0 s; 2.5 s] for the enthalpy ℎ (left), the mass fraction (center) and the velocity (right). We see that the scheme has actually the behaviour of a relaxed asymptotic-preserving scheme, meaning that even if the initial conditions are not the same between the two models, after some time the solution of the 4-lmnc model converges to the solution of the 3-lmnc model up to an error of order . Figure 5 displays the error calculated at time = 2.5 s for the enthalpy ℎ, the mass fraction and the velocity as a function to of in loglog-scale. We see that the scheme is at order 1 in for all variables.

Spatial coupling of the two regimes
In this section, we investigate the spatial coupling of non-instantaneous and the instantaneous regimes, and compare qualitatively the results with a classic compressible test case (coupling spatially HRM and HEM Figure 5. Convergence in of the asymptotic-preserving scheme. models). To this end, we chose a similar physical setting as in [4]. More precisely, we set = 80 m, which represents a part of the coolant circuit containing the core on the left half of the domain. The power density is thus located in the region 18 m ≤ ≤ 23 m, and is set in this region to the constant value Φ = 75 × 10 6 W · m −3 . The computation is done on a coarse grid, and the discretisation parameter is ∆ = 0.4. The CFL condition is equal to 0.99, and in this test case, the time step ∆ is always greater than 0.03.
Concerning the thermodynamical parameters, it is well-known that the values computed in [39] or [40] are not adapted to PWR simulations. In [19], we showed the drawbacks related to these values for the pressure of interest * = 155 bar (compared to experimental values [41]), and we introduced a new strategy to compute an incomplete EoS which is exact at saturation and rich enough to close the 3-lmnc model. Here another complete set of parameters is computed by fitting to experimental data around a reference point adapted to PWR simulations. First, we take ℎ sat , sat and sat from NIST experimental values at = * : sat = 617.939 K and sat l = 1.68243 × 10 −3 m 3 · kg −1 sat g = 9.81065 × 10 −3 m 3 · kg −1 ℎ sat l = 1.629 × 10 6 J · K −1 ℎ sat g = 2.596 × 10 6 J · K −1 .
Then, we choose to set = 0 and , to the values of [39]: ,l = 4.268 × 10 3 J · K −1 · kg −1 ,g = 1.4874 × 10 3 J · K −1 · kg −1 . On Figure 6, we plot the liquid mass fraction 1 − at time = 3 s, as well as the one for = 1 in the whole domain and the one for = 10 −10 in the whole domain. The asymptotic-preserving scheme handles perfectly the coupling, without having to determine artificial boundary conditions on the coupling interface. Moreover, we can observe the influence of the mass transfer relaxation time when = 1. Of course, comparisons with compressible test cases are limited to low Mach number regimes, where no shock waves appear.

Conclusion
In this paper, we derived the low Mach number approximation of a 4-equation HRM system, with equality of velocities, pressures and temperatures between phases. This model takes into account a relaxation term allowing the exchange of mass between phases through a given characteristic time. In order to describe properly the thermodynamics of the fluid, we derived an equation of state for the model, and compared it with the one at saturation (used in the 3-lmnc model). One particularly important feature of this equation of state is the fact that it is not defined piecewise, avoiding discontinuities which may lead to large numerical errors due the approximation of the physical parameters.
For the 4-lmnc model, we were able to prove some results on the analytical steady-state solution, as well as positivity properties on the unknowns and the relaxation term. A well-balanced scheme has been designed, which preserves the properties of a steady-state solution for any value of the relaxation parameter (thus being possibly far from the explicit steady-state solution of the 3-lmnc model). This allows in particular to be extremely accurate on the flux conservation, as well as on the slope of the enthalpy. A similar scheme is also given in Appendix for the 3-lmnc model. This work could be a preliminary study to more realistic 2D/3D computations, taking into account the coupling with a simplified neutronics model, as it in [17,18].
Finally, we studied the formal convergence of the 4-lmnc model towards the 3-lmnc one, when the relaxation parameter tends to zero. An asymptotic-preserving scheme has also been derived, for which we proved and observed numerically that it provides, after some time iterations, a numerical solution close to the one of the 3-lmnc model up to an error of order . This scheme can be particularly useful for the coupling of spatial regions where the unrelaxed 4-equation model is needed with spatial regions where the 3-equation HEM model may be used (instantaneous relaxation of the chemical potentials).
A natural extension of this paper will be to investigate the low Mach number models that can be obtained from the hierarchy of models mentioned in the introduction, with possible disequilibria in pressures, temperatures, chemical potentials and velocity. This shall be the topic of further works.
-if ℎ −1 ≤ ℎ sat g ≤ ℎ , by using the definition of * g , we have -otherwise, ℎ and ℎ −1 are in the same phase, thus (ℎ ) = (ℎ −1 ), so that We can conclude that, in any case, equation (A.1b) is equivalent to We can thus prove the two well-balanced properties. We first compute which completes the proof.