David Le Touzé

Zhe Li
email: zhe.li@ec-nantes.fr

Guillaume Oger

A finite volume WENO scheme for immiscible inviscid two-phase flows

Keywords: Two-phase finite volume method, WENO reconstruction, approximate Riemann flux solver, free-surface flows, oscillation-free

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Multiphase flows are widely present in nature as well as in industrial applications, such as ocean waves, sloshing/impacting problems in a liquid tank and atomization of a fuel jet, etc. This type of fluid flows is usually immiscible and of high density ratios between the two phases, which makes it challenging for computational fluid dynamics (CFD) methods to handle high density ratios and to maintain sharp phase-interfaces in a relatively efficient way.

Given these physical characteristics of two-phase flows, it is of primary importance to choose an appropriate mathematical model to describe the necessary physics, which is expected to be mathematically well-posed and easy to solve using a numerical scheme. Insightful reviews on the mathematical models for two-phase flows can be found in [START_REF] Drew | Mathematical modelling of two-phase flow[END_REF][START_REF] Ransom | Hyperbolic two-pressure models for two-phase flow[END_REF][START_REF] Stewart | Two-Phase Flow: Models and Methods[END_REF]. In the present work, we adopt a weakly-compressible two-phase model based on a hyperbolic system of three equations, which are the mass conservation equation for both phases and the momentum conservation equation for the two-phase mixture, closed with a barotropic equation of state (EOS).

Notice that, although without energy equation, this reduced three-equation system is sufficient for representing the physics in the targeted weakly-compressible two-phase flows, in which the physical compressibility is not of interest. Thus, highly-compressible multiphase flows, such as strong shock waves, are not considered in the scope of the present work. Initially tested by Chanteperdrix [START_REF] Chanteperdrix | Modélisation et simulation numérique d'écoulements diphasiques à interface libre. Application à l'étude des mouvements de liquides dans les résservoirs de véhicules spatiaux[END_REF], this model has been recently applied by Grenier et al. [START_REF] Grenier | An accurate low-Mach scheme for a compressible two-fluid model applied to free-surface flows[END_REF] for simulating two-phase free-surface flows with an implicit low-Mach scheme. In this two-phase model, it is assumed that both phases have the same velocity and the pressure equilibrium is instantaneously accomplished, i.e. at any given location the pressure of the two phases are always the same. As stated by Drew [START_REF] Drew | Mathematical modelling of two-phase flow[END_REF], this instantaneous pressure equilibrium is an adequate assumption for the fluid flows where the speed of sound in each phase is large compared with velocities of interest, which is the case in the present work for weakly-compressible flows. More importantly, it can be proven that the adopted 'single-pressure & singlevelocity' model is a well-posed hyperbolic system [START_REF] Stewart | Two-Phase Flow: Models and Methods[END_REF], hence it has less mathematical difficulties to solve using a numerical scheme. Recently, this assumption has also been applied in a general formulation for various types of multiphase flows [START_REF] Saurel | A general formulation for cavitating, boiling and evaporating flows[END_REF]. Here it is worth noting that surface tension and viscosity effects are not considered in the present paper and will be studied in the future work. The reason for this choice is that without physical surface tension and viscous terms in the governing equations, it is straightforward to assess the properties such as the stability and the numerical dissipation of the proposed scheme.

In the present work, we consider to apply a high-order (≥ 3) scheme in order to have a good compromise between the CPU time and precision [START_REF] Bigay | A weakly-compressible cartesian grid approach for hydrodynamic flows[END_REF]. To do so, we propose a finite volume (FV) formulation based on a 5 thorder weighted essentially non-oscillatory (WENO) scheme for simulating weakly-compressible two-phase flows.

Initially proposed by Liu et al. [START_REF] Liu | Weighted essentially non-oscillatory schemes[END_REF], WENO schemes [START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF] have been widely adopted for capturing discontinuities with appropriate stencils and for attaining desired high order accuracies in smooth regions. Promising results have been obtained in several recent works [START_REF] Coralic | Finite-volume WENO scheme for viscous compressible multicomponent flows[END_REF][START_REF] Dumbser | High order space-time adaptive ADER-WENO finite volume schemes for non-conservative hyperbolic systems[END_REF][START_REF] Han | Adaptive multi-resolution method for compressible multi-phase flows with sharp interface model and pyramid data structure[END_REF][START_REF] Johnsen | Implementation of WENO schemes in compressible multicomponent flow problems[END_REF][START_REF] Shu | High order WENO and DG methods for time-dependent convection-dominated PDEs: A brief survey of several recent developments[END_REF][START_REF] Rodriguez | A high-order accurate five-equations compressible multiphase approach for viscoelastic fluids and solids with relaxation and elasticity[END_REF][START_REF] Shukla | An interface capturing method for the simulation of multi-phase compressible flows[END_REF][START_REF] Wang | An incremental-stencil WENO reconstruction for simulation of compressible two-phase flows[END_REF], in which WENO schemes were used for handling high-Mach compressible multiphase flows in the presence of the interaction between phase interfaces and shock waves. One of the difficulties in these types of numerical simulations is that pressure oscillations may occur near the phase interface, which are not only non-physical (purely numerical error) but also may cause numerical instabilities. Indeed, as stated by Abgrall et al. [START_REF] Abgrall | How to prevent pressure oscillations in multicomponent flow calculations: a quasi conservative approach[END_REF][START_REF] Abgrall | Computations of compressible multifluids[END_REF], it is generally a major difficulty for a multiphase numerical scheme to avoid this type of pressure oscillations, even in the case of the advection of an isolated material interface with uniform pressure and velocity fields, which is named as 'the Abgrall test-case' hereafter in the present paper.

In 2006, Johnsen and Colonius [START_REF] Johnsen | Implementation of WENO schemes in compressible multicomponent flow problems[END_REF] have proposed a FV-WENO scheme which can avoid pressure oscillations near the material interface in the Abgrall test-case, by reconstructing the primitive variables instead of the conservative variables in the WENO reconstruction step. The desired formal convergence order of WENO scheme (3 rd -and 5 th -orders) have been obtained in a 1D test-case. In 2014, Coralic and Colonius [START_REF] Coralic | Finite-volume WENO scheme for viscous compressible multicomponent flows[END_REF] have shown that, despite of the advantage of oscillation-free feature, if one directly reconstructs the primitive variables [START_REF] Johnsen | Implementation of WENO schemes in compressible multicomponent flow problems[END_REF], the convergence order of the WENO scheme is limited to 2 nd -order in a 2D vortex advection test-case. In order to avoid pressure oscillations near material interfaces and to retain the convergence order of the WENO scheme simultaneously, Coralic and Colonius [START_REF] Coralic | Finite-volume WENO scheme for viscous compressible multicomponent flows[END_REF] proposed a reconstruction procedure, summarized as: (i) first-time WENO reconstruction of the conservative variables to obtain the primitive variables at the Gauss quadrature points within each FV-cell; (ii) Gauss quadrature to have the cell-average values of the primitive variables; (iii) second-time WENO reconstruction of the primitive variables to get the fluid state at each cell-interface. Notice that, in the dimension-by-dimension fashion [START_REF] Titarev | Finite-volume WENO scheme for three-dimensional conservation laws[END_REF], each WENO reconstruction requires twice 'sweeps' in 2D cases in each direction in order to get the interpolated values at cell-interfaces. Thanks to this reconstruction procedure, they succeeded in obtaining a 5 th -order convergence in a 2D vortex advection test-case.

Here we propose a FV-WENO scheme which can avoid pressure oscillations and retain the convergence order, using only once the WENO reconstruction. First, we choose to reconstruct the conservative variables in order to retain the convergence order of the WENO scheme. Second, we summarize two necessary conditions for ensuring the oscillation-free feature of a numerical scheme in the Abgrall test-case

• Condition-I: Given uniform (constant in space) velocity and pressure fields, the velocity and pressure reconstructed by means of the WENO scheme at cell-interfaces should remain constant;

• Condition-II: With the numerical flux computed by an approximate Riemann flux solver, the velocity and pressure fields at next instant (updated using an time integration scheme) should remain constant in space.

One may notice that it is easy to satisfy Condition-I when using the primitive variables (pressure and velocity) in the reconstruction step. However, this choice will degrade the convergence order of WENO schemes down to 2 nd -order in multi-dimensional cases, as shown in [START_REF] Coralic | Finite-volume WENO scheme for viscous compressible multicomponent flows[END_REF]. In the present work, we propose a modified WENO scheme satisfying Condition-I by means of: (i) the use of a common smoothness indicator for the conservative variables ρ1 and ρ2 ; (ii) the reconstruction of the mixture density ρ for deriving the velocity at cell-interfaces. The main advantages of the proposed scheme are its relatively simple algorithmics and small computational efforts it requires. Nevertheless, the proposed numerical framework cannot be directly applied to simulate violent and challenging cases involving significant volume changes, which is hence out of the scope of the present paper. In addition, inspired of the works of Tokareva and Toro [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF][START_REF] Tokareva | A flux splitting method for the Baer-Nunziato equations of compressible two-phase flow[END_REF], we utilize the generalized Riemann invariants (GRI) [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF] to provide an analysis on the mathematical property of the adopted three-equation hyperbolic system and then propose a linearized two-phase Riemann flux solver based on the GRI. Furthermore, the classical 4 th -order Runge-Kutta scheme is applied to integrate in time the semi-discrete system of equations in an explicit way. It is shown that the proposed linearized Riemann flux solver and the 4 th -order Runge-Kutta time integrator satisfy Condition-II, i.e. the constant velocity and pressure fields are retained at next instant. As a consequence, the proposed FV-WENO scheme satisfies the above two conditions, which means that the oscillation-free feature is ensured. It is here worth noting that, as did in [START_REF] Coralic | Finite-volume WENO scheme for viscous compressible multicomponent flows[END_REF][START_REF] Titarev | Finite-volume WENO scheme for three-dimensional conservation laws[END_REF], two Gauss quadrature points are used for the surface integral of flux at cell-interfaces. Despite the lower theoretical accuracy order that could be expected, we have observed a 5th-order convergence of the proposed scheme in the advection of a multiphase vortex test-case, which is consistent with the results presented in [START_REF] Coralic | Finite-volume WENO scheme for viscous compressible multicomponent flows[END_REF][START_REF] Titarev | Finite-volume WENO scheme for three-dimensional conservation laws[END_REF].

The rest of the paper is organized as follows. Section 2 presents the adopted three-equation mathematical model for describing weakly-compressible two-phase flows with a barotropic closure law. The finite volume formulation is briefly presented in Section 3. In Section 4, the reason of pressure oscillations in the present numerical framework is firstly given, followed by the proposed modifications to the WENO reconstruction scheme. Based on the GRI, a linearized Riemann flux solver is presented in Section 5. A short discussion on the stability issues and some limiters is given in Section 6. The proposed FV-WENO scheme is validated using a series of numerical test-cases in the presence of two-phase free-surface flows, which are shown in Section 7.

Finally, the conclusions are drawn in Section 8.

Mathematical models

A three-equation model for weakly-compressible two-phase flows

In the present work, we adopt a three-equation model [START_REF] Chanteperdrix | Modélisation et simulation numérique d'écoulements diphasiques à interface libre. Application à l'étude des mouvements de liquides dans les résservoirs de véhicules spatiaux[END_REF][START_REF] Grenier | An accurate low-Mach scheme for a compressible two-fluid model applied to free-surface flows[END_REF] for describing weakly-compressible inviscid two-phase fluid flows. The system of governing equations consists of the mass conservation equation for the two phases and the momentum conservation equation for the mixture, given as

             ∂ ρ1 ∂t + ∇ • (ρ 1 u) = 0, ∂ ρ2 ∂t + ∇ • (ρ 2 u) = 0, ∂ (ρu) ∂t + ∇ • (ρu ⊗ u + pI) = ρg, (1)
where ρk (x, t) = α k (x, t)ρ k (x, t) with α k (x, t) and ρ k (x, t) being respectively the volume fraction and density of the k th phase (k = 1, 2) at the location x and time t. By definition, one has

k α k (x, t) = α 1 (x, t)+α 2 (x, t) = 1.
In addition, u(x, t) denotes the velocity of the mixture, which is shared by both of the two phases, as shown in Equation [START_REF] Abgrall | How to prevent pressure oscillations in multicomponent flow calculations: a quasi conservative approach[END_REF]. The density and pressure of the mixture are defined as ρ(x, t) = k α k (x, t)ρ k (x, t) and p(x, t) = k α k (x, t)p k (x, t), respectively. p k (x, t) denotes the individual pressure of the k th phase, which is related to ρ k (x, t) by means of the chosen equation of state (EOS) discussed subsequently. I is the identity matrix and g is the body force per unit mass.

Moreover, the instantaneous pressure equilibrium assumption is adopted, which means that the individual pressures of the two phases are always the same, i.e. p(x, t) = p 1 (x, t) = p 2 (x, t), ∀x ∈ Ω f , ∀t ∈ [0, +∞). As commented by Drew [START_REF] Drew | Mathematical modelling of two-phase flow[END_REF], this assumption is adequate for the cases where the speed of sound is sufficiently large compared to the velocity of interest, or in other words, the M a number is sufficiently small, which is the case for weakly-compressible fluid flows. As a result, the present system of governing equations can be categorized as a 'single-pressure & single-velocity' two-phase model [START_REF] Stewart | Two-Phase Flow: Models and Methods[END_REF], which is hyperbolic with the adopted EOS, as shown subsequently.

Barotropic closure law

The individual pressure p k is related to each density ρ k by means of an EOS. As did in [START_REF] Grenier | An accurate low-Mach scheme for a compressible two-fluid model applied to free-surface flows[END_REF], we adopt a barotropic EOS in order to close the system of equations [START_REF] Abgrall | How to prevent pressure oscillations in multicomponent flow calculations: a quasi conservative approach[END_REF], which is given as

p k = p 0 + c 2 k (ρ k -ρ k,0 ), (2)
where p 0 , ρ k,0 and c k are three constants denoting the reference pressure, the reference density and the speed of sound for the k th phase, respectively, with k = 1, 2. This barotropic EOS can be adopted as no thermal effects are considered in the present work, since there is no energy or temperature equation in the system of governing equations (1).

More importantly, with this barotropic closure law, the adopted three-equation two-phase model is a hyperbolic system possessing real eigenvalues, which are related to the speed of propagation of information (the speed of sound) in the fluid domain.

Speed of sound in the two-phase mixture

In the present weakly-compressible multiphase model, the speed of sound in the two-phase mixture, noted as c(x, t), can be expressed as

c = 1 ρ α 1 ρ 1 c 2 1 + α 2 ρ 2 c 2 2 , (3)
of which the derivation is provided in Appendix A. As an example, Figure 1 shows the evolution of the speed of sound in an air-water mixture, as the volume fraction of water α 1 varies from 0 to 1. One may observe that when the volume fraction α 1 (x, t) is close to 0.5, the local speed of sound in the mixture can be much smaller than any of the two individual speeds of sound c 1 and c 2 . Similar features of the speed of sound in the mixture for other multiphase models have been reported in [START_REF] Le Martelot | Towards the direct numerical simulation of nucleate boiling flows[END_REF][START_REF] Murrone | A five equation reduced model for compressible two phase flow problems[END_REF][START_REF] Saurel | A general formulation for cavitating, boiling and evaporating flows[END_REF][START_REF] Stewart | Two-Phase Flow: Models and Methods[END_REF]. This phenomenon has also been observed experimentally by Costigan and Whalley [START_REF] Costigan | Measurements of the speed of sound in air-water flows[END_REF] in their measurements of speed of sound in air-water flows.

Equilibrium volume fraction

As shown in Equation ( 1), there is no advection equation for the volume fraction field α k (x, t) in the system of governing equations. This means that no extra equation is needed to update the volume fraction at each new instant. As performed in [START_REF] Grenier | An accurate low-Mach scheme for a compressible two-fluid model applied to free-surface flows[END_REF], given ρ1 and ρ2 , we can determine the volume fraction α 1 so that the pressure equilibrium condition is ensured, by solving the following equation

p 0 + c 2 1 ρ1 α 1 -ρ 1,0 = p 0 + c 2 2 ρ2 1 -α 1 -ρ 2,0 . (4)
Once α 1 is obtained, the volume fraction of the second phase is computed by α 2 = 1 -α 1 . It is then easy to prove that with positive ρ1 and ρ2 , there exists a unique valid equilibrium volume fraction α 1 by solving Equation (4). As did in [START_REF] Chanteperdrix | Modélisation et simulation numérique d'écoulements diphasiques à interface libre. Application à l'étude des mouvements de liquides dans les résservoirs de véhicules spatiaux[END_REF], let us rewrite Equation (4) as

ρ2 c 2 2 γ 2 -(q -q)γ -ρ1 c 2 1 = 0, (5)
where γ = α 1 /(1 -α 1 ), q = ρ 2,0 c 2 2 -ρ 1,0 c 2 1 and q = ρ2 c 2 2 -ρ1 c 2 1 . Given ρ1 > 0 and ρ2 > 0, this quadratic equation [START_REF] Bigay | A weakly-compressible cartesian grid approach for hydrodynamic flows[END_REF] possesses two real solutions

         γ + = q -q + (q -q) 2 + 4ρ 2 c 2 2 ρ1 c 2 1 2ρ 2 c 2 2 > 0, γ -= q -q -(q -q) 2 + 4ρ 2 c 2 2 ρ1 c 2 1 2ρ 2 c 2 2 < 0, (6)
of which γ + is the unique valid solution, since α 1 must satisfy the condition 0 < α 1 < 1. Finally, the unique equilibrium volume fraction is obtained as α 1 = γ + /(1 + γ + ).

Finite volume formulation

In the present work, the finite volume method is adopted to solve the governing equations (1

where U, K, n and S denote respectively the vector of conservative variables, the flux matrix, the outward normal vector on the surface Γ e of Ω e and the vector of source term, given as

U =         ρ1 ρ2 ρu ρv         , K = F G =         ρ1 u ρ1 v ρ2 u ρ2 v ρuu + p ρuv ρvu ρvv + p         , n =   n x n y   and S =         0 0 ρg x ρg y         . ( 8
)
In addition, F and G denote the flux vectors in the x and y directions respectively. Equation ( 8) is expressed here in 2D but its extension to 3D is straightforward.

In a fixed Cartesian mesh, Equation ( 7) can be written for the mesh cell over

[x i-1/2 , x i+1/2 ]×[y j-1/2 , y j+1/2 ]
as follows

dU i,j dt = - 1 ∆x∆y    y j+1/2 y j-1/2 F(x i+1/2 , y, t) dy - y j+1/2 y j-1/2 F(x i-1/2 , y, t) dy    - 1 ∆x∆y    x i+1/2 x i-1/2 G(x, y j+1/2 , t) dx - x i+1/2 x i-1/2 G(x, y j-1/2 , t) dx    + S i,j , (9)
where U i,j = U i,j (t) and S i,j = S i,j (t) denote respectively the cell-average values of the conservative variables and the source term, defined as

                   U i,j (t) = 1 ∆x∆y y j+1/2 y j-1/2 x i+1/2 x i-1/2 U(x, y, t) dx dy, S i,j (t) = 1 ∆x∆y y j+1/2 y j-1/2 x i+1/2 x i-1/2
S(x, y, t) dx dy, [START_REF] Ducrozet | CN-Stream: Open-source library for nonlinear regular waves using stream function theory[END_REF] with ∆x = x i+1/2 -x i-1/2 and ∆y = y j+1/2 -y j-1/2 being the cell spacing in x-and y-directions, respectively. Before carrying out the time integration of Equation [START_REF] Drew | Mathematical modelling of two-phase flow[END_REF], one needs to evaluate the surface integrals on the right hand side. As suggested in [START_REF] Coralic | Finite-volume WENO scheme for viscous compressible multicomponent flows[END_REF][START_REF] Titarev | Finite-volume WENO scheme for three-dimensional conservation laws[END_REF], a two-point Gauss quadrature rule is utilized to compute the surface integral of flux at each cell boundary, for example

y j+1/2 y j-1/2 F(x i+1/2 , y, t) dy σ w σ F(x i+1/2 , y σ , t), (11)
where w σ = ∆y/2 denotes the quadrature weight for the σ th -Gauss point at y σ = y j ± ∆y/(2 √ 3) and the flux function is evaluated with an approximate Riemann flux solver, i.e. F(x i+1/2 , y σ , t)

F(U L σ , U R σ ) with U L σ and U R
σ being the reconstructed states by WENO scheme at the left and right sides of the cell-interface.

In the present work, we apply the classical 4 th -order Runge-Kutta time integrator to Equation ( 9) in order to update the cell-average state to the next instant, i.e. from U

n i,j = U i,j (t n ) to U n+1 i,j = U i,j (t n+1
). Because of the explicit feature of the scheme, the time-step ∆t is chosen as [START_REF] Bigay | A weakly-compressible cartesian grid approach for hydrodynamic flows[END_REF] ∆t = CFL min e (∆x, ∆y)

max e (|u| + c) , (12)
where the CFL number is set equal to 0.7 in all the numerical test-cases presented in the following. As an explicit weakly-compressible scheme [START_REF] Vittoz | Comparisons of weakly-compressible and truly incompressible approaches for viscous flow into a high-order Cartesian-grid finite volume framework[END_REF], the present method has a certain number of advantages, such as there is no need to solve a Poisson equation, it is easy to implement and it usually maintains good scaling properties, although the time-step is restricted because of the explicit feature.

WENO reconstruction scheme

As mentioned previously, the WENO reconstruction scheme is used to interpolate the variables at each cellinterface with high-order accuracy. In this section, the WENO reconstruction step of the proposed multiphase scheme, which satisfies Condition-I given in Introduction, will be presented with details.

Reconstructed variables

At the beginning of each time-step, the only available data are the newly updated cell-average values of the conservative variables which are obtained by means of the time integration. Hence, it is straightforward to reconstruct the conservative variables at each cell-interface using the cell-average states within the stencils.

For example, in 1D cases with uniform mesh spacing, let φ(x) be any one of the conservative variables ρ1 , ρ2 and ρu, the 5 th -order accurate reconstructed values of φ(x) on the left side of the cell-interface at x L i+1/2 can be expressed as

[φ] L i+1/2 = ω 0 6 2φ i-2 -7φ i-1 + 11φ i + ω 1 6 -φ i-1 + 5φ i + 2φ i+1 + ω 2 6 2φ i + 5φ i+1 -φ i+2 , (13)
where ω l denotes the weighting coefficient for the l th -stencil (l = 0, 1, 2) and ω l = d l when there is no WENO smoothness effect with d l being the ideal weighting coefficient which are equal to d 0 = 1/10, d 1 = 3/5 and

d 2 = 3/10 for [φ] L i+1/2
. Moreover, we note hereafter [•] × as the WENO reconstruction operator that gives the interpolated value of the variable '•' on the '×' side of the cell-interface at the position ' '. Now, we shall demonstrate that, when there is no WENO smoothness effect, i.e. ω l = d l , the reconstructed state at the cell-interface satisfies Condition-I.

To do so, let us take the left side of the cell-interface at x L i+1/2 as an example. Using the formula given in Equation ( 13), the reconstructed variables ρ1 = α 1 ρ 1 and ρ2 = α 2 ρ 2 can be obtained and rewritten as

         [ρ 1 ] L i+1/2 = m ζ m (α 1 ρ 1 ) m , [ρ 2 ] L i+1/2 = m ζ m (α 2 ρ 2 ) m , (14)
where ζ m denotes the coefficient for the m th -cell-average value (α k ρ k ) m within the three stencils with

m = i -2, i -1, i, i + 1, i + 2. With the condition ω l = d l , one can have ζ i-2 = d 0 /3 = 1/30, ζ i-1 = -7d 0 /6 -d 1 /6 = -13/60, ζ i = 11d 0 /6 + 5d 1 /6 + d 2 /3 = 47/60, ζ i+1 = d 1 /3 + 5d 2 /6 = 9/20 and ζ i+2 = -d 2 /6 = -1/20. Note
here that, since there is no smoothness effect, ζ m is entirely determined by the mesh spacing and is the same

for obtaining [ρ 1 ] L i+1/2 and [ρ 2 ] L i+1/2 . Given a uniform pressure field, e.g. p(x, t) = p 0 , ∀x ∈ Ω f at the current instant t, one has ρ 1 = ρ 1,0 and ρ 2 = ρ 2,0 , hence          [ρ 1 ] L i+1/2 = m ζ m (α 1 ) m ρ 1,0 , [ρ 2 ] L i+1/2 = m ζ m (α 2 ) m ρ 2,0 . (15)
These reconstructed values

[ρ 1 ] L i+1/2 and [ρ 2 ]
L i+1/2 are then used to get the equilibrium volume fraction α L 1,i+1/2 at the cell-interface by solving Equation (4) as follows

c 2 1 [ρ 1 ] L i+1/2 α L 1,i+1/2 -ρ 1,0 = c 2 2 [ρ 2 ] L i+1/2 1 -α L 1,i+1/2 -ρ 2,0 , (16)
which, with the help of Equation ( 15) and (α 2 ) m = 1 -(α 1 ) m , can be written as

c 2 1    ρ 1,0 m ζ m (α 1 ) m α L 1,i+1/2 -ρ 1,0    = c 2 2    ρ 2,0 -ρ 2,0 m ζ m (α 1 ) m 1 -α L 1,i+1/2 -ρ 2,0    , (17)
of which one obvious solution is

α L 1,i+1/2 = m ζ m (α 1 ) m . ( 18
)
Since it is previously shown that there is a unique valid solution, Equation ( 18) is just the sought solution for the equilibrium volume fraction on the left side of the cell-interface at x L i+1/2 . It can then be easily verified that the pressure at the cell-interface is equal to the constant value, i.e. p L i+1/2 = p 0 .

Similarly, in order to evaluate the velocity u L i+1/2 at x L i+1/2 , one needs to firstly get the reconstructed conservative variable ρu as

[ρu] L i+1/2 = m ζ m (ρu) m , (19)
which gives the velocity at the cell-interface as

u L i+1/2 = [ρu] L i+1/2 ρ L i+1/2 = [ρu] L i+1/2 [ρ 1 ] L i+1/2 + [ρ 2 ] L i+1/2 = m ζ m (ρu) m m ζ m (α 1 ρ 1 ) m + m ζ m (α 2 ρ 2 ) m = m ζ m (ρu) m m ζ m (ρ) m . ( 20
)
Given a uniform velocity field, e.g. u(x, t) = u 0 , ∀x ∈ Ω f , the velocity at the cell-interface is then evaluated as

u L i+1/2 = u 0 m ζ m (ρ) m m ζ m (ρ) m = u 0 . ( 21
)
As a consequence, in the absence of smoothness effect, given uniform pressure and velocity fields, the original WENO reconstruction scheme satisfies Condition-I, which means that it can ensure that the reconstructed pressure and velocity at each cell-interface have constant values. Similarly, one can obtain the same conclusion

for p R i-1/2 and u R i-1/2 .

Smoothness effect

The idea of WENO reconstruction scheme is to assign a weighting coefficient for each stencil, based on its smoothness, such that the reconstruction essentially utilizes the information from the relatively smooth stencils.

Using a well defined smoothness indicator, a WENO scheme based on three stencils can retain a 5 th -order accuracy in smooth areas of the domain, while being able to provide a lower order reconstructed value at the locations where there exists a strong discontinuity. In the present work, we adopt the widely used WENO smoothness indicators given by Jiang and Shu [START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF], defined as

β l = x i+1/2 x i-1/2 ∆x ∂P l ∂x 2 + ∆x 3 ∂ 2 P l ∂x 2 2 dx, (22)
where β l and P l denote the smoothness measurement and the local 2 nd -order polynomial for the l th -stencil. For l = 0, 1 and 2, one has

             β 0 = 13 12 φ i-2 -2φ i-1 + φ i 2 + 1 4 φ i-2 -4φ i-1 + 3φ i 2 , β 1 = 13 12 φ i-1 -2φ i + φ i+1 2 + 1 4 φ i-1 -φ i+1 2 , β 2 = 13 12 φ i -2φ i+1 + φ i+2 2 + 1 4 3φ i -4φ i+1 + φ i+2 2 , (23)
with which the weight ω l is then obtained applying

ω l = ξ l ξ 0 + ξ 1 + ξ 2 , with ξ l = d l (β l + ε) 2 , ( 24
)
where ε is a small enough value in case of β l = 0, which is set as ε = 10 -6 as suggested by Jiang and Shu [START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF]. It is worth noting here that the smoothness indicators β l in Equation ( 23) are not dimensionless. Hence using the same value of ε for all the variables, i.e. ρ1 , ρ2 and ρu, may generate some problems for choosing the appropriate stencils [START_REF] Fedkiw | Simplified discretization of systems of hyperbolic conservation laws containing advection equations[END_REF]. In the present work, we propose to renormalize β l with a characteristic quantity φ 2 0 for each variable, in such a way that the weight is now computed by

ω l = ξl ξ0 + ξ1 + ξ2 , with ξl = d l βl + ε 2 , and βl = β l φ 2 0 , (25)
where φ 0 = ρ 1,0 , ρ 2,0 and ρ 1,0 u 0 for the reconstructed variables ρ1 , ρ2 and ρu, respectively.

Using the weighting coefficient given in Equation ( 25) for the WENO reconstruction scheme [START_REF] Embid | Mathematical analysis of a two-phase continuum mixture theory[END_REF], one can avoid utilizing the information from the stencils containing strong discontinuities. However, as shown in the following, such WENO scheme cannot ensure Condition-I, when considering the smoothness effect.

As mentioned previously, each reconstructed variable possesses its own WENO smoothness indicator β l , e.g. let us assign β1,l for ρ1 and β2,l for ρ2 . Because β1,l and β2,l are generally not identical, although they are indeed related via the volume fraction field, the weight are generally not the same, i.e. ω 1,l = ω 2,l . Given a uniform pressure field (p = p 0 ), the reconstructed ρ1 and ρ2 at the cell-interface can then be expressed as

         [ρ 1 ] L i+1/2 = m ζ 1,m (α 1 ) m ρ 1,0 , [ρ 2 ] L i+1/2 = m ζ 2,m (α 2 ) m ρ 2,0 , (26)
where one should note that ζ 1,m = ζ 2,m in the general case.

To obtain the volume fraction

α L 1,i+1/2 , [ρ 1 ] L i+1/2 and [ρ 2 ] L i+1/2
given by Equation ( 26) are then substituted in Equation ( 16), which gives

c 2 1    ρ 1,0 m ζ m,1 (α 1 ) m α L 1,i+1/2 -ρ 1,0    = c 2 2    ρ 2,0 -ρ 2,0 m ζ m,2 (α 1 ) m 1 -α L 1,i+1/2 -ρ 2,0    , (27)
of which the solution

α L 1,i+1/2 is generally not equal to m ζ 1,m (α 1 ) m or m ζ 2,m (α 1 ) m , since ζ 1,m = ζ 2,m .
Hence the pressure at the cell-interface x L i+1/2 is not equal to the constant value, i.e. p L i+1/2 = p 0 . As a consequence, a pressure oscillation may occur from this cell-interface.

Besides, given a uniform velocity field (u = u 0 ), the interface velocity u L i+1/2 is obtained with Equation [START_REF] Lax | Hyperbolic systems of conservation laws II[END_REF] as

u L i+1/2 = u 0 m ζ ρ,m (ρ) m m ζ 1,m (α 1 ρ 1 ) m + m ζ 2,m (α 1 ρ 1 ) m , (28)
where ζ ρ,m denotes the coefficient related to the smoothness of the mixture density field ρ(x, t). Obviously, one may observe that u L i+1/2 = u 0 , because the coefficients ζ ρ,m , ζ 1,m and ζ 2,m are generally not the same.

WENO smoothness indicators in the proposed method

With the previous analysis on the smoothness effects, in the present work, we propose a modified WENO reconstruction scheme, which is based on (i) Common smoothness indicator for ρ1 and ρ2 . As shown previously, the pressure oscillation is basically due to the fact the β1,l is generally not equal to β2,l . Hence, we propose a common dimensionless smoothness indicator βc,l , computed as βc,l = β1,l β2,l ,

which is a reasonable choice, because βc,l embodies the smoothness in both ρ1 and ρ2 fields which are correlated through the volume fraction α 1 , especially when the two densities ρ

u L i+1/2 = [ρu] L i+1/2 [ρ] L i+1/2 = u 0 m ζ ρ,m (ρ) m m ζ ρ,m (ρ) m = u 0 , (30)
which means that the velocity at the cell-interface has the constant value u 0 . Notice that two candidates for the value of the mixture density ρ L i+1/2 at the cell-interface are available, which are

[ρ] L i+1/2 and [ρ 1 ] L i+1/2 + [ρ 2 ] L i+1/2 .
In smooth regions, the difference between these two values is of order O(∆x 5 ). In order to ensure a strict oscillation-free property, it is suggested to use [ρ] L i+1/2 exclusively for evaluating the velocity at the cell-interface.

Once again, the same conclusion can be obtained for p R i+1/2 and u R i+1/2 . Consequently, with the above two modifications (i) and (ii), the proposed WENO scheme can satisfy Condition-I. The influence on the convergence order is discussed with details subsequently.

Convergence order

One of the advantages of the WENO scheme is that the numerical solver can choose the relatively smooth stencils, based on the smoothness indicators, by means of a unified formula for the whole computational domain.

Hence, the smoothness effect should not degrade the high convergence order in smooth areas. As shown in the following, the proposed WENO scheme can retain the original convergence order in smooth areas.

As demonstrated by Jiang and Shu [START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF], the 5 th -order convergence of WENO scheme can be retained, if

ω l = d l + O ∆x 2 , ( 31
)
where ω l is given in Equation ( 25) and d l is the optimal weight. To satisfy Equation [START_REF] Shu | High order WENO and DG methods for time-dependent convection-dominated PDEs: A brief survey of several recent developments[END_REF], it suffices to have

βl = D 1 + O ∆x 2 , ( 32
)
where βl denotes the smoothness indicator given in Equation ( 25) and D is a non-zero quantity independent of l.

Equation ( 32) can be easily verified with a Taylor expansion of Equation ( 23). However, Taylor expansion cannot be directly carried out, since φ i denotes the cell-average value of the function φ, but not the value of φ at x i . Nevertheless, one can assume there exists a function ϕ(x), such that

ϕ i = ϕ(x i ) = 1 ∆x xi+1/2∆x xi-1/2∆x φ(x) dx = φ i , (33)
with which the Taylor expansion for β0 , β1 and β2 can be carried out at x i as

                 β0 = 1 φ 2 0 ϕ 2 i ∆x 2 + 1 φ 2 0 13 12 ϕ 2 i - 2 3 ϕ i ϕ i ∆x 4 - 1 φ 2 0 13 6 ϕ i ϕ i - 1 2 ϕ i ϕ i ∆x 5 + O ∆x 6 , β1 = 1 φ 2 0 ϕ 2 i ∆x 2 + 1 φ 2 0 13 12 ϕ 2 i + 1 3 ϕ i ϕ i ∆x 4 + O ∆x 6 , β2 = 1 φ 2 0 ϕ 2 i ∆x 2 + 1 φ 2 0 13 12 ϕ 2 i - 2 3 ϕ i ϕ i ∆x 4 + 1 φ 2 0 13 6 ϕ i ϕ i - 1 2 ϕ i ϕ i ∆x 5 + O ∆x 6 . (34)
Obviously, dividing the smoothness indicator β l by φ 2 0 still allows one to express βl in the form as shown in Equation [START_REF] Shukla | An interface capturing method for the simulation of multi-phase compressible flows[END_REF], which means this renormalization does not degrade the convergence order of the WENO scheme. Now it is necessary to assess whether the common smoothness indicator βc,l can retain the convergence order.

To do so, let us express the smoothness indicators β1,l and β2,l for the variables ρ1 and ρ2 , respectively, as

   β1,l = D 1 1 + O ∆x 2 , β2,l = D 2 1 + O ∆x 2 , ( 35
)
where D 1 and D 2 are two quantities independent of l. Hence, the common smoothness indicator βc,l can then be expressed as βc,l = β1,l β2,l = D 1 D 2 (1 + O (∆x 2 )), [START_REF] Tokareva | A flux splitting method for the Baer-Nunziato equations of compressible two-phase flow[END_REF] which, by means of a Taylor expansion analysis, can be rewritten as

βc,l = D c 1 + O ∆x 2 , with D c = D 1 D 2 , (37)
which means that the common dimensionless smoothness indicator βc,l proposed in Equation ( 29) satisfies the sufficient condition [START_REF] Shukla | An interface capturing method for the simulation of multi-phase compressible flows[END_REF].

Hence, in the proposed WENO scheme, all the variables can be reconstructed at the cell-interface with 5 th -order accuracy, i.e.

[φ]

L i+1/2 = φ(x i+1/2 ) + O(∆x 5
), where φ = ρ1 , ρ2 , ρ, ρu, ρv. Furthermore, the proposed scheme can avoid oscillations when given uniform pressure and velocity fields in an efficient way, compared with the strategy proposed in [START_REF] Coralic | Finite-volume WENO scheme for viscous compressible multicomponent flows[END_REF] requiring twice the WENO reconstructions.

Approximate Riemann flux solver

Given the reconstructed variables U L σ and U R σ at both sides of each cell-interface, a numerical flux solver F(U L σ , U R σ ) is then utilized to evaluate the flux F at the cell-interface, as shown in Equation [START_REF] Dumbser | FORCE schemes on unstructured meshes II: Non-conservative hyperbolic systems[END_REF]. In the present work, we propose a linearized Riemann flux solver, based on the generalized Riemann invariants.

Eigenstructure of the hyperbolic system

To estimate the flux at cell-interface in the x-direction, one usually analyzes the quasi-linear x-split governing equations expressed as

∂W ∂t + A ∂W ∂x = 0, (38)
where W = [ρ 1 , ρ2 , u, v] . In addition, the matrix A is given as

A =          u 0 ρ1 0 0 u ρ2 0 c2 1 ρ c2 2 ρ u 0 0 0 0 u          , with        c1 = ∂p ∂ ρ1 , c2 = ∂p ∂ ρ2 , (39)
so that the pressure gradient has been expressed as

∂p ∂x = ∂p ∂ ρ1 ∂ ρ1 ∂x + ∂p ∂ ρ2 ∂ ρ2 ∂x = c2 1 ∂ ρ1 ∂x + c2 2 ∂ ρ2 ∂x . (40)
The matrix A possesses the following eigenvalues

           λ 1 = u -c, λ 2 = λ 3 = u, λ 4 = u + c, (41)
where c denotes the speed of sound in the two-phase mixture, computed as

c = ρ1 c2 1 + ρ2 c2 2 ρ , (42)
which is equivalent to Equation (3), as shown in Appendix A. As the individual density ρ k and the volume fraction α k are positive in Equation ( 3), c has a real value, so do the eigenvalues. Consequently, the present three-equation two-phase model is a hyperbolic system, of which the corresponding right eigenvectors are given as

r 1 =         ρ1 ρ2 -c 0         , r 2 =         c2 2 -c 2 1 0 0         , r 3 =         0 0 0 1         and r 4 =         ρ1 ρ2 c 0         . ( 43
)

Generalized Riemann invariants

Based on the simple wave theory [START_REF] Lax | Hyperbolic systems of conservation laws II[END_REF], the generalized Riemann invariants (GRI) condition [START_REF] Embid | Mathematical analysis of a two-phase continuum mixture theory[END_REF][START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF] states the fact that an invariant relation holds across each eigenvalue-related wave in the x-t plane, as shown in Figure 2. For the present hyperbolic system, λ 1 and λ 4 are related respectively to two genuinely non-linear waves, and λ 2 and λ 3 are both related to a contact wave which is linearly degenerated. In the recent works of Tokareva and Toro [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF][START_REF] Tokareva | A flux splitting method for the Baer-Nunziato equations of compressible two-phase flow[END_REF], the GRI condition has been successfully applied to solve the Baer-Nunziato equations [START_REF] Baer | A two-phase mixture theory for the Deflagration-to-Detonation Transition (ddt) in reactive granular materials[END_REF] for compressible flows. Across the λ a -related wave, the variation dW a is parallel to the corresponding right eigenvector r a in the phase space, which gives

• across the λ 1 -related wave:

dρ 1 ρ1 = dρ 2 ρ2 = du -c = dv 0 , hence            ρ1 du + c dρ 1 = 0, ρ2 du + c dρ 2 = 0, dv = 0, (44)
• across the λ 4 -related wave:

dρ 1 ρ1 = dρ 2 ρ2 = du c = dv 0 , hence            ρ1 du -c dρ 1 = 0, ρ2 du -c dρ 2 = 0, dv = 0, (45)
• across the λ 2 -and λ 3 -related wave (contact wave):

dρ 1 c2 2 = dρ 2 -c 2 1 = du 0 = dv 0 and dρ 1 0 = dρ 2 0 = du 0 = dv 1 , hence    c2 1 dρ 1 + c2 2 dρ 2 = 0, du = 0. ( 46
)
By analyzing the GRI conditions ( 44), ( 45) and ( 46), we summarize several properties of the present threeequation two-phase hyperbolic system, given as follows (i) since dv = 0 across the λ 1 -and λ 4 -related waves, there is no variation of v across these two waves, i.e.

v * ,L = v L and v * ,R = v R ;

(ii) the velocity component u is constant in the star region, because across the contact wave, one has du = 0 which indicates that u * ,L = u * ,R = u * ;

(iii) the pressure p is constant in the star region, because across the contact wave, the variation of the pressure is equal to zero, since

dp = ∂p ∂ ρ1 dρ 1 + ∂p ∂ ρ2 dρ 2 = c2 1 dρ 1 + c2 2 dρ 2 = 0, (47)
which gives p * ,L = p * ,R = p * ;

(iv) with the adopted EOS (2) and the property (iii), the density for each phase remains constant across the contact wave, i.e. ρ * ,L k = ρ * ,R k = ρ * k , with k = 1, 2.

Linearized two-phase Riemann flux solver

Based on the GRI conditions, we propose a linearized Riemann flux solver, inspired of the work in [START_REF] Tokareva | A flux splitting method for the Baer-Nunziato equations of compressible two-phase flow[END_REF].

The idea is to linearize the GRI conditions across the λ 1 -and λ 4 -waves, which are both related to the speed of sound c, in order to obtain an approximate state in the star region of the two-phase Riemann problem at each cell-interface.

Given Equation ( 44), the variation of pressure across the λ 1 -wave can be expressed as

dp = c2 1 dρ 1 + c2 2 dρ 2 , = c2 1 - ρ1 c du + c2 2 - ρ2 c du , = - ρ1 c2 1 + ρ2 c2 2 c du, = -ρc du. (48)
Similarly, based on Equation (45), across the λ 4 -wave, one has dp = ρc du. Now, we adopt a linearization step by taking ρc ρ L c L across the λ 1 -wave and ρc ρ R c R across the λ 4 -wave, which gives

   p L -p * + ρ L c L (u L -û * ) = 0, p R -p * -ρ R c R (u R -û * ) = 0, (49)
where p * and û * denote respectively the approximate values of p * and u * , which can then be computed as

       p * = ρ L c L ρ R c R (u L -u R ) + (ρ R c R p L + ρ L c L p R ) ρ L c L + ρ R c R , û * = (p L -p R ) + (ρ L c L u L + ρ R c R u R ) ρ L c L + ρ R c R . ( 50
)
As shown in the property (iv) given previously, we can get the density for each phase in the star region with the EOS (2), i.e. ρ * k = ρ k,0 + (p * -p 0 )/c 2 k , with k = 1, 2. Besides, from the GRI conditions ( 44) and (45), one can obtain    ρ du + c dρ = 0, across the λ 1 -wave, ρ du -c dρ = 0, across the λ 4 -wave, (51) which can be linearized, giving

   ρ L (û * -u L ) + c L (ρ * ,L -ρ L ) = 0, ρ R (û * -u R ) -c R (ρ * ,R -ρ R ) = 0, (52)
where ρ * ,L and ρ * ,R denote the approximate density of the mixture at the left and right star region, respectively, which are given as

       ρ * ,L = ρ L - ρ L c L (û * -u L ), ρ * ,R = ρ R + ρ R c R (û * -u R ). (53)
Now, given ρ * k , ρ * ,L and ρ * ,R , we can obtain the volume fraction in the star region by means of the definition of the density of mixture

   ρ * ,L = α * ,L 1 ρ * 1 + (1 -α * ,L 1 )ρ * 2 , ρ * ,R = α * ,R 1 ρ * 1 + (1 -α * ,R 1 )ρ * 2 .
(54)

Finally, we have the approximate states of variables within the whole star region, i.e. Û * ,L and Û * ,R , which allow us to obtain the numerical flux Fσ at the σ th -Gauss point as

Fσ =    F( Û * ,L σ ), if û * ≥ 0, F( Û * ,R σ ), if û * < 0, (55)
where F(U) is the flux function given in Equation [START_REF] Costigan | Measurements of the speed of sound in air-water flows[END_REF]. It is here worth noting that Equation ( 55) is sufficient for weakly-compressible two-phase flows, because the speed of sound is much greater than the fluid velocity. In case of local hypersonic flow, we need to include another two situations

Fσ =    F(U L ), if S L > 0, F(U R ), if S R < 0, (56)
where S L and S R denote the estimation of the speed of wave related to the eigenvalues λ 1 = u-c and λ 4 = u+c, as shown in Figure 2. In the literature [START_REF] Batten | On the choice of wavespeeds for the HLLC Riemann solver[END_REF][START_REF] Toro | Restoration of the contact surface in the HLL-Riemann solver[END_REF], several ways have been proposed for estimating the speeds of waves, e.g.

S L = min(u L -c L , u R -c R ) and S R = max(u L + c L , u R + c R ).
Now, we shall show that the proposed linearized two-phase Riemann flux solver satisfies Condition-II given in Introduction, i.e. by means of the present numerical flux and the adopted RK4 time integrator, the velocity and pressure of the updated state at the next instant still remain constant.

To do so, let us recall that the proposed WENO scheme satisfies Condition-I, hence one has p L = p R = p 0 and u L = u R = u 0 at all the cell-interfaces. By means of Equation (50), one can obtain that p * = p 0 and û * = u 0 , which gives ρ * ,L = ρ L and ρ * ,R = ρ R with Equation (53). In addition, since p * = p 0 , one has ρ * 1 = ρ 1,0 and ρ * 2 = ρ 2,0 with the EOS. Consequently, by solving Equation (54), one can have α * ,L

1 = α L 1 and α * ,R 1 = α R 1 .
Hence, the numerical flux for the i th -cell at x i-1/2 and x i+1/2 are computed as (e.g.

u 0 > 0) Fi-1/2 = F( Û * ,L i-1/2 ) =      α L 1,i-1/2 ρ 1,0 u 0 α L 2,i-1/2 ρ 2,0 u 0 ρ L i-1/2 u 2 0 + p 0      and Fi+1/2 = F( Û * ,L i+1/2 ) =      α L 1,i+1/2 ρ 1,0 u 0 α L 2,i+1/2 ρ 2,0 u 0 ρ L i+1/2 u 2 0 + p 0      , ( 57
)
with which the first intermediate state U

(1) i

for the RK4 scheme can be obtained as

U (1) i =      (α 1 ρ 1 ) (1) i (α 2 ρ 2 ) (1) i (ρu) (1) i      =      (α 1 ρ 1 ) n i (α 2 ρ 2 ) n i (ρu) n i      - ∆t 2∆x           α L 1,i+1/2 ρ 1,0 u 0 α L 2,i+1/2 ρ 2,0 u 0 ρ L i+1/2 u 2 0 + p 0      -      α L 1,i-1/2 ρ 1,0 u 0 α L 2,i-1/2 ρ 2,0 u 0 ρ L i-1/2 u 2 0 + p 0           =      (α 1 ) (1) i ρ 1,0 (α 2 ) (1) i ρ 2,0 (ρ) (1) i u 0      , ( 58
) with              (α 1 ) (1) i = (α 1 ) n i - ∆t 2∆x α L 1,i+1/2 -α L 1,i-1/2 u 0 , (α 2 ) (1) i = (α 2 ) n i - ∆t 2∆x α L 2,i+1/2 -α L 2,i-1/2 u 0 , (ρ)
(1) i have constant pressure and velocity fields. Hence the proposed linearized Riemann flux solver and the adopted time integration scheme satisfy Condition-II.

= (ρ) n i - ∆t 2∆x ρ L i+1/2 -ρ L i-1/2 u 0 . (

Stability and limiters

As presented in Section 2, in order to prove the uniqueness of the solution of the volume fraction, one needs to verify that ρ1 and ρ2 are both positive. Indeed, maintaining the positivity of ρ1 and ρ2 are of great importance in the present scheme, not only because of their physical meaning, but also due to the fact that negative values of ρ1 or ρ2 could not give real eigenvalues, which sometimes leads to ill-posed problems with unphysical instabilities [START_REF] Ransom | Hyperbolic two-pressure models for two-phase flow[END_REF]. However, it is well known that the conventional WENO scheme cannot guarantee that the reconstructed values of ρ1 or ρ2 are positive at cell-interfaces. In the literature, there exist several positivity-preserving schemes and limiters [START_REF] Perthame | On positivity preserving finite volume schemes for Euler equations[END_REF][START_REF] Zhang | On maximum-principle-satisfying high order schemes for scalar conservation laws[END_REF][START_REF] Zhang | Positivity-preserving high order finite difference WENO schemes for compressible Euler equations[END_REF].

In the present work, for the sake of efficiency, we choose to apply a simple limiter for the reconstructed ρ1 and ρ2 at each Gauss point. For both of the two phases

(k = 1, 2), if the reconstructed value [ρ k ] L σ or [ρ k ] R
σ is less than a critical value 0.9ε α ρ k,0 , where ε α = 10 -4 , then it is set to be the cell-average value (ρ k ) i,j . Although this operation might theoretically degrade the convergence order, it is observed that in practice this limiter is activated only in areas where there exist strong discontinuities of the volume fraction field α 1 . In the two subsequent test-cases, it is shown that this limiter does not affect 5 th -order convergence of the present WENO scheme.

Besides, we observe that numerical instabilities may occur in relatively violent two-phase flows, in the presence of solid wall boundaries, such as the dam-break flow. This is because the cells near solid walls have only one available stencil for the WENO reconstruction. Using the stencil selection algorithm based on the comparison of smoothness, one always takes this stencil for reconstruction, even though it contains huge discontinuities, which is the case for the dam-break flow, especially when the water front impacts on the solid wall at the corner of the tank. Hence, for this type of two-phase flows, we adopt a limiter similar as the one used in [START_REF] Titarev | Finite-volume WENO scheme for three-dimensional conservation laws[END_REF], which sets the reconstructed value [φ] L σ to be the cell-average value (φ

) i,j , if [ρ k ] L σ -(ρ k ) i,j > 0.9(ρ k ) i,j .

Numerical results

In the following test-cases, all the two-phase flows have a density ratio of 1000, i.e. ρ 1,0 /ρ 2,0 = 1000, k = 1 for water with ρ 1,0 = 1000 kg/m 3 and k = 2 for air with ρ 2,0 = 1 kg/m 3 . In addition, as the physical compressibility is not of primary importance in all following test-cases, an artificial speed of sound in each fluid c k is chosen according to the weakly-compressible hypothesis, by comparing with the characteristic fluid velocity in each test-case, for the sake of numerical efficiency. As presented in Introduction, viscous effects are not considered in the governing equations, hence in the following numerical examples, either there is no viscous term in the analytical solution or the viscosity effects are not dominant in the considered physical problems. Accordingly, the free-slip boundary condition is imposed on solid walls, if there are. This is realized by using a mirrored state, e.g. U R σ , when computing the numerical flux on the solid wall, which possesses the same densities and tangential component but opposite normal velocity component as the ones of the reconstructed state U L σ .

Interface advection

1D advection

In order to verify numerically the oscillation-free feature of the proposed scheme in the Abgrall test-case, we first consider an 1D interface advection, of which the initial condition is given as

α 1 (x, 0) =      1 -ε α , if L 2 -r 0 < x < L 2 + r 0 , ε α , otherwise, (60)
where ε α = 10 -4 , L = 1 m denotes the length of the 1D computational domain with x ∈ [0, L] and r 0 = L/10

is the radius of the 1D liquid droplet. In addition, ρ 1 (x, 0) = ρ 1,0 = 1000 kg/m 3 and ρ 2 (x, 0) = ρ 2,0 = 1 kg/m 3 , so that p(x, 0) = p 0 = 0 Pa. The speed of sound for both fluids are set as c 1 = c 2 = 70 m/s and the whole fluid domain is initialized with a uniform velocity field, i.e. u(x, 0) = u 0 = 1 m/s. Periodic boundary conditions are imposed at the two extremities, at x = 0 and x = L.

Figure 3 shows the density fields of the mixture at the initial (t = 0 s) and the final (t = 1 s) instants of the numerical simulation, obtained with the proposed scheme using 64 uniform cells for spacial discretization.

The 1D droplet of water moves from left to right and gets back to its initial position after 1 s. It can be clearly observed that the phase-interface becomes smeared or diffused as time evolves. Figure 4 shows the pressure and velocity oscillations in the numerical results obtained with the present and original WENO schemes. We observe that the present WENO scheme can ensure zero oscillation in the pressure and velocity fields. As for the computational efficiency, we observe that the present WENO scheme is approximately 1.65 times more costly than the original WENO scheme. The extra cost is majorly due to the reconstruction of the mixture's density and the computation of the common WENO smoothness for ρ1 and ρ2 . As mentioned in Introduction, the objective of reconstructing conservative variables is to attain high convergence order, i.e. 5 th -order in the present work. In order to carry out a convergence study, we initialize the shape of the 1D water droplet with a smoother function

α 1 (x, 0) = (1 -2ε α ) exp - |x -L/2| 2 r 2 0 + ε α , (61)
which is shown in Figure 5-(a). Notice that all the other variables remain the same as in the previous sharp interface advection test-case. The convergence study is carried out using five mesh resolutions with the number of cells N = 64, 128, 256, 512 and 1024. The L 1 , L 2 and L ∞ errors of the mixture's density at t = 1 s are computed as

                         L 1 -Error = 1 N i ρ num i -ρ ref i ρ 1,0 , L 2 -Error = 1 N i ρ num i -ρ ref i 2 ρ 2 1,0 , L ∞ -Error = max i ρ num i -ρ ref i ρ 1,0 , (62)
where ρ num i and ρ ref i denote respectively the numerical and reference cell-average values of the mixture's density of the i th -cell. Table 1 shows the errors and the convergence orders evaluated by computing Order = ln Error q Error q+1 ln(

where Error q denotes the error corresponding to the q th -mesh resolution. As shown in Table 1 and Figure 5-(b), the present numerical scheme possesses the expected 5 th -order convergence in this 1D advection test-case. Next, we assess the oscillation-free feature of the proposed numerical scheme in a 2D interface advection test-case. The initial volume fraction field α 1 within the (i, j) th -cell is imposed as where x i,j and y i,j denote the coordinates in x-and y-directions of the barycenter of each cell, and x 0 and y 0 the coordinates of the center of the computational domain [-L, L] × [-L, L] with L = 1 m. As in the 1D test-case, periodic conditions are imposed at the four boundaries. The initial uniform velocity and pressure fields are set as u(x, y, 0) = 2 m/s, v(x, y, 0) = 2 m/s and p(x, y, 0) = 0 Pa. The speed of sound is chosen as

α 1 =    1 -ε α , if (x i,j -x 0 ) 2 + (y i,j -y 0 ) 2 < r 2 0 , ε α , otherwise, (64)
c 1 = c 2 = 30 m/s.
Figure 6 shows the initial and final density fields of the mixture in the simulation using a mesh of size 64∆x × 64∆y. From the pressure field at the final instant as shown in Figure 7, one can observe that a strict zero pressure oscillation is ensured by the proposed scheme in this 2D interface advection test-case, as in the previous 1D advection.

Multiphase vortex advection

As presented in [START_REF] Coralic | Finite-volume WENO scheme for viscous compressible multicomponent flows[END_REF], a numerical scheme possessing a 5 th -order convergence property in 1D test-cases might fail to converge with a 5 th -order in multi-dimensional cases, such as the one in [START_REF] Johnsen | Implementation of WENO schemes in compressible multicomponent flow problems[END_REF]. Hence, it is necessary to assess the convergence order of the present scheme in 2D where the gradients of all variables have non-zero values. For this purpose, we carry out a convergence study with a 2D multiphase vortex advection test-case, of which the analytical solution is provided in Appendix B.

The initial volume fraction and pressure fields at t = 0 s within the periodic computational domain [-L, L]× [-L, L] are imposed as

                               α 1 (r, θ, 0) = (1 -2ε α ) exp - r 2 r 2 0 + ε α , p(r, θ, 0) = P 1 -exp - r 2 r 2 0 , u(r, θ, 0) = U -sin θ r ρ ∂p ∂r , v(r, θ, 0) = V + cos θ r ρ ∂p ∂r , (65)
where r and θ denote the cylindrical coordinates with x = r cos θ and y = r sin θ, ε α = 10 -4 , r 0 = 0.2L with L = 1 m, and P = 0.1ρ 2,0 c 2 2 with c 1 = c 2 = 30 m/s. The density of both phases ρ 1 and ρ 2 can be obtained with the EOS. Finally, U = V = 2 m/s denote the two components of the advection velocity. The initial volume fraction field and the relative velocity vector with respect to the moving frame are shown in Figure 8-(a). Table 2 shows the errors and convergence orders of different norms computed using the formulas (62) and (63). In addition, based on these values of errors, the convergence slopes are also given in Figure 8-(b). As expected, a 5 th -order convergence is obtained in this 2D multiphase vortex advection test-case using the proposed scheme. Here it is worth noting that, in this convergence study, all physical fields including the volume fraction α 1 should be sufficiently smooth. In the presence of discontinuities, the convergence order will be limited to 1 st -order at most.

Linear sloshing

Now, we consider a 2D linear sloshing test-case [START_REF] Chanteperdrix | Modélisation et simulation numérique d'écoulements diphasiques à interface libre. Application à l'étude des mouvements de liquides dans les résservoirs de véhicules spatiaux[END_REF][START_REF] Grenier | An accurate low-Mach scheme for a compressible two-fluid model applied to free-surface flows[END_REF], of which the configuration is shown in Figure 9.

Initially, air and water are both at rest under the gravity effect in a closed tank which is composed of four free-slip solid walls. As t > 0, a constant horizontal acceleration a is applied to the whole fluid domain in the positive x-direction, which is equivalent to the case where the tank suddenly moves in the negative x-direction with an acceleration -ae x . The geometrical and physical parameters are given in Table 3.

where

         κ n = nπ L , ω n = gκ n (ρ 1,0 -ρ 2,0 ) ρ 1,0 coth(κ n h 1 ) + ρ 2,0 coth(κ n h 2 ) . (67)
Figure 10 shows the time histories of the air-water interface elevation on the left and right walls of the tank, which are obtained using three uniform meshes of size 64∆x × 144∆y, 128∆x × 288∆y and 256∆x × 576∆y. In addition, the shapes of the air-water interface, plotted by means of the contour line of α 1 = 0.5 at four instants t = 1 s, 2 s, 3 s and 4 s, are provided in Figure 11. A good agreement can be observed between the present numerical result and the analytical solution, which shows the capacity of the proposed two-phase FV-WENO scheme to simulate accurately free-surface flows of high density ratio in such type of linear sloshing cases with small wave amplitude.

Wave propagation

The objective of the following test-case is to assess the capability of the present weakly-compressible twophase finite volume scheme for handling long-period wave propagations near the free-surface (air-water interface), which is a challenging task for most either incompressible or compressible two-phase CFD methods due to the numerical diffusion effects introduced by the numerical schemes.

The configuration of this test-case is shown in Figure 12 and the parameters are given in Table 4. We the left and right hand sides of the computational domain. Free-slip boundary conditions are imposed on the upper and lower solid walls of the channel. Based on the parameters given in Table 4, the analytical solution can be obtained using the non-linear potential flow theory [START_REF] Ducrozet | CN-Stream: Open-source library for nonlinear regular waves using stream function theory[END_REF][START_REF] Rienecker | A Fourier approximation method for steady water waves[END_REF], which is given as

             η(x, t) = a 0 2 + +∞ n=1 a n cos (κ n (x -c p t)) , Φ(x, y, t) = b 0 (x -c p t) + +∞ n=1 b n cosh (κ n (y + h)) cosh(κ n h) sin (κ n (x -c p t)) , (68)
where η and Φ denote respectively the free-surface elevation and the potential function, the wave number κ n = 2πn/λ and the phase velocity c p = λ/T with T 0.7018 s being the wave propagation period obtained from the non-linear potential flow theory. Once the coefficients a n and b n are determined, the velocity and pressure fields can then be obtained as

               u(x, y, t) = ∂Φ(x, y, t) ∂x , v(x, y, t) = ∂Φ(x, y, t) ∂y , p(x, y, t) ρ 1,0 = R -gy - 1 2 c 2 p - ∂Φ(x, y, t) ∂t - u 2 (x, y, t) + v 2 (x, y, t) 2 , ( 69
)
where R is the Bernoulli constant. (ρ 1,0 /ρ 2,0 = 1000) between water and air, the numerical solution in the water is expected not to be significantly affected by this choice of initial condition.

The computational domain is discretized with a uniform mesh of size 1024∆x×512∆y. This spatial resolution leads to λ = 102.4∆x and H = 24.576∆y with ∆x/∆y = 3.3675. The total physical simulation time is set as 40T , which is sufficiently long to assess the numerical diffusion. In the present work, this test-case has been carried out using the proposed WENO scheme and the 2 nd -order accurate MUSCL [START_REF] Van Leer | Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method[END_REF] scheme with the 'minmod' slope-limiter. Both of the two numerical results are compared with the analytical solution obtained

with the non-linear potential flow theory.

Figure 13 shows the shapes of the free-surface at each 10 periods (10T ), which are drawn using the contour line of α 1 = 0.5. In this comparison, one may observe that the present WENO scheme almost retains the wave amplitude even after 40 periods (40T ), although the shape of the free-surface is slightly shifted, comparing to the analytical solution. Whereas in the result obtained with the MUSCL scheme, the wave amplitude is significantly damped and a large phase error can also be observed.

As for the computational cost, it is found that the numerical computation using the present WENO scheme is approximately 9.3 times slower than the one with the MUSCL scheme. In order to have a meaningful and fair comparison between the present WENO and MUSCL schemes, we carry out another simulation using the MUSCL scheme with a refined mesh (Mesh2: 2048∆x × 1024∆y), which costs almost the same as the present WENO scheme with the previous coarser mesh (Mesh1: 1024∆x × 512∆y). More precisely, in these conditions the WENO scheme is just 1.079 times slower than the MUSCL one. Figure 14 gives the numerical result of the MUSCL scheme using the refined mesh, which is compared with the previous ones. One may observe that the MUSCL scheme with the refined mesh still introduces more numerical diffusion than the proposed WENO scheme. This difference is clearly shown in Figure 15, where a fast Fourier transform with a moving window of size T has been carried out on the free-surface elevation η m (t) = η(L/2, t) in the middle of the domain. By comparing the first two harmonic amplitudes, it can be easily observed that the present WENO scheme is less diffusive and provides more accurate result.

In this specific test-case involving long-time wave propagation, the present 5 th -order WENO scheme provides promising results in good agreement with the analytical solution and appears to be superior to the adopted 2 nd -order MUSCL scheme. Nevertheless, it is worth noting that MUSCL schemes are often more robust in the presence of strong discontinuities such as shock waves.

Dam-break flow

At last, we tried to apply the proposed weakly-compressible two-phase finite volume scheme to simulate a 2D dam-break flow involving moderately more violent impacts and higher dynamics than the previous test-cases. Figure 16 shows the initial configuration of this test-case and the computational parameters are given in Table 5. Initially, a column of water at rest under the gravity is located at the right lower corner of a closed tank which is composed of four free-slip walls. The fluid domain is discretized with a uniform mesh of size 256∆x × 128∆y.

Water Air

Figure 17 shows the volume fraction field α 1 of water at three instants, which are compared with the experimental measurements obtained by Lobovsky et al. [START_REF] Lobovský | Experimental investigation of dynamic pressure loads during dam break[END_REF]. As time evolves, the water column falls down under the gravity effect and then impacts on the opposite solid wall. In addition, the water front position and the pressure signal at the left lower corner (Sensor 1 in [START_REF] Lobovský | Experimental investigation of dynamic pressure loads during dam break[END_REF]) are shown in Figure 18. From these comparisons between the present numerical results and the experimental measurements, one may observe that the proposed scheme provides reasonably accurate results in good agreement with the experimental ones. Nevertheless, some pressure oscillations can be found, as shown in Figure 18-(b), for the pressure sensor at the left lower corner of the tank.

This pressure oscillation is due to the fact that only one stencil (even not smooth) is available for the WENO reconstruction at this corner where huge gradients can be present, especially when the water front impacts on the vertical solid wall. Improving the numerical stability of the present scheme near solid walls is of great importance and will be considered in future works.

Conclusions

In the present paper, a finite volume WENO scheme is proposed for simulating weakly-compressible twophase flows. Based on a three-equation hyperbolic system of equations, the semi-discrete equations in the present method are explicitly integrated in time using the Runge-Kutta scheme. A modified WENO reconstruction scheme is proposed, which consists in using a common smoothness indicator for ρ1 and ρ2 and reconstructing the mixture's density ρ for obtaining the velocity at cell-interfaces. In addition, based on the linearization of the generalized Riemann invariants, an approximate two-phase Riemann flux solver is proposed to compute the numerical flux at cell-interfaces. It is then shown that the proposed method can retain the formal convergence order of WENO scheme in multi-dimensional cases and avoid strictly pressure oscillations in the Abgrall testcase. In the numerical examples of free-surface flows such as the linear sloshing and the wave propagation test-cases, it is observed that the proposed scheme introduces a very small numerical diffusion and predicts correctly the behaviors of free-surface waves even during a long-time simulation. In the dam-break flow testcase, the present FV-WENO scheme provides a numerical result in good agreement with the experimental data. Some pressure oscillations occur at the corner of the tank, due to the fact that only one stencil is available at this place for the WENO reconstruction. In the future work, viscous and surface tension terms will be considered for more complex two-phase flows.

Figure 1 :

 1 Figure 1: Speed of sound in an air-water mixture with ρ 1 = 1000 kg/m 3 , ρ 2 = 1 kg/m 3 , c 1 = 1500 m/s and c 2 = 350 m/s.

1 and ρ 2

 2 vary little under the weakly-compressible hypothesis. Using this common smoothness indicator βc,l , the coefficients ζ 1,m and ζ 2,m in Equation (27) become identical, i.e. ζ 1,m = ζ 2,m , which gives a unique solution of α L 1,i+1/2 that verifies p L i+1/2 = p 0 . (ii) Reconstruction of the mixture density ρ. In order to avoid the velocity oscillation, we propose to reconstruct the mixture density at the cell-interface x L i+1/2 , instead of using the sum of [ρ 1 ] L i+1/2 and [ρ 2 ] L i+1/2 , which is then used to get the velocity u L i+1/2 at the cell-interface as

Figure 2 :

 2 Figure 2: Wave structure of the present hyperbolic system.

Figure 3 :Figure 4 :

 34 Figure 3: 1D interface advection test-case: density of the mixture at t = 0 s and t = 1 s.

Figure 5 :

 5 Figure 5: Convergence study in the 1D advection test-case: (a) initial volume fraction field; (b) convergence order.

Figure 6 :

 6 Figure 6: Density fields of the mixture in the 2D interface advection test-case: (a) t = 0 s; (b) t = 1 s.

Figure 7 :

 7 Figure 7: Pressure field at the final instant (t = 1 s) in the 2D interface advection test-case.

Figure 8 :

 8 Figure 8: 2D multiphase vortex advection test-case: (a) initial configuration in the moving frame; (b) convergence order.

Figure 9 :

 9 Figure 9: Initial configuration of the linear sloshing test-case.

Figure 10 :

 10 Figure 10: Time history of the interface elevations η L = η(0, t) and η R = η(L, t) on the left and right walls in the linear sloshing test-case with three mesh resolutions: 64∆x × 144∆y, 128∆x × 288∆y and 256∆x × 576∆y.

Figure 11 :

 11 Figure 11: Shape of the air-water interface in the linear sloshing test-case using the mesh resolution 256∆x×576∆y at four instants: (a) t = 1 s; (b) t = 2 s; (c) t = 3 s; (d) t = 4 s.

Figure 12 :

 12 Figure 12: Configuration of the wave propagation test-case.

Figure 13 :Figure 14 :Figure 15 : 1 with A ana 1 2 2 with A ana 2 3

 1314151122 Figure 13: Shape of the free-surface in the wave propagation test-case using the proposed WENO scheme and the MUSCL scheme with the same mesh of resolution 1024∆x × 512∆y at four instants: (a) t = 10T ; (b) t = 20T ; (c) t = 30T ; (d) t = 40T .

Figure 16 : 2

 162 Figure 16: Initial configuration of the dam-break flow test-case.

Figure 17 :

 17 Figure 17: Comparison of flow patterns between the numerical results (left) obtained with the proposed WENO scheme and the experimental measurements (right) carried out by Lobovsky et al. [23] (reprinted with permission from Elsevier) at three dimensionless instants t * = t g/H: (a) t * 1.58; (b) t * 3.27; (c) t * 6.67.

Figure 18 :

 18 Figure 18: Water front position and pressure signal in the dam-break flow test-case: (a) time history of the water front position; (b) pressure signal at the left lower corner of the tank.

  59)Equation(58) shows that, after one intermediate time integration, α 1 , α 2 and ρ vary in time, while ρ 1 , ρ 2 (hence p) and u remain constant for all the cells. It can then be easily proven that all the intermediate states as well as the state at the next time-step U

	n+1
	i

Table 1 :

 1 Numerical errors in the 1D advection test-case.

	64	8.7542e-4	-	1.5857e-3	-	4.7826e-3	-
	128	3.7302e-5	4.5527	6.8293e-5	4.5372	2.1892e-4	4.4493
	256	1.2148e-6	4.9405	2.2702e-6	4.9108	8.7785e-6	4.6403
	512	3.7243e-8	5.0276	7.1105e-8	4.9967	2.9355e-7	4.9023
	1024	1.0935e-9	5.0899	2.1816e-9	5.0265	8.8813e-9	5.0467
	7.1.2. 2D advection

Number of cells L

1 -Error L 1 -Order L 2 -Error L 2 -Order L ∞ -Error L ∞ -Order

Table 2 :

 2 Numerical errors in the 2D multiphase vortex advection test-case.Number of cells L1 -Error L 1 -Order L 2 -Error L 2 -Order L ∞ -Error L ∞ -Order

	64	3.9112e-4	-	1.2847e-3	-	8.8276e-3	-
	128	1.9138e-5	4.3531	6.0041e-5	4.4193	4.5100e-4	4.2908
	256	6.2019e-7	4.9476	1.9616e-6	4.9358	1.5642e-5	4.8496
	512	1.7145e-8	5.1769	5.6494e-8	5.1178	4.6175e-7	5.0822

Table 3 :

 3 Geometrical and physical parameters of the linear sloshing test-case.

		Parameter				Value
		Tank width (L)			1 m
		Initial water height (h 1 )		1 m
		Initial air height (h 2 )		1.25 m
		Gravity (g)				9.81 m/s 2
		Horizontal acceleration (a)	0.01g
		Speed of sound in water (c 1 )	100 m/s
		Speed of sound in air (c 2 )	100 m/s
	Based on the linearized potential flow theory, the analytical solution of the interface elevation η(x, t) is given
	as [6, 16]
	η(x, t) =	a g	x -	L 2	+	+∞ n=0	4 2n+1 Lκ 2	cos(ω 2n+1 t) cos(κ 2n+1 x) ,

  The initial state in the water is imposed by means of Equation (69) with t = 0. However, since there is no analytical solution available in the air, in the present work, we choose to adopt a simple strategy in which both the pressure and velocity fields are initially equal to zero in the air. Due to the high density ratio

		0.015						0.015
				Present WENO				Present WENO
		0.01			Analytical		0.01			Analytical
		0.005						0.005
	y/L	0					y/L	0
		-0.005						-0.005
		-0.01						-0.01
				t = 1 s						t = 2 s
		-0.015						-0.015
		0	0.2	0.4	0.6	0.8	1	0	0.2	0.4	0.6	0.8	1
				x/L						x/L
				(a)						(b)
		0.015						0.015
				Present WENO				Present WENO
		0.01			Analytical		0.01			Analytical
		0.005						0.005
	y/L	0					y/L	0
		-0.005						-0.005
		-0.01						-0.01
				t = 3 s						t = 4 s
		-0.015						-0.015
		0	0.2	0.4	0.6	0.8	1	0	0.2	0.4	0.6	0.8	1
				x/L						x/L
				(c)

Table 4 :

 4 Geometrical and physical parameters of the wave propagation test-case.

	Parameter	Value
	Length of the computational domain (L)	8.082 m
	Water depth (h)	0.6 m
	Height of the channel	2h
	Wave length (λ)	0.1L
	Wave height (H)	0.0576 m
	Gravity (g)	9.81 m/s 2
	Speed of sound in water (c 1 )	30 m/s
	Speed of sound in air (c 2 )	30 m/s

Acknowledgement

This work was performed by using HPC resources of the Centrale Nantes Supercomputing Centre on the cluster Liger.

Appendix A: Speed of sound in the two-phase mixture

The speed of sound of the mixture is given in Equation [START_REF] Wang | High-order CFD methods: current status and perspective[END_REF], or as

in which c2 1 and c2 2 are defined in Equation [START_REF] Van Leer | Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method[END_REF]. Using the instantaneous pressure equilibrium condition p = p 1 = p 2 and the barotropic EOS (2), one can have

Meanwhile, the equality p 1 = p 2 gives

where α 1 (or α 2 = 1 -α 1 ) can be considered as a function of ρ1 and ρ2 . Now, taking the partial derivative of Equation (A.3) with respect to ρ1 yields

which gives

Similarly, following the same procedure for the second phase (k = 2), one can obtain

Finally, substituting Equation (A.6) and Equation (A.7) into Equation (A.1) gives

which can then be written as

from which one can easily find the speed of sound of the mixture Equation (3).

Appendix B: Analytical solution of the multiphase vortex advection test-case

Similar to the strategies proposed in [START_REF] Dumbser | FORCE schemes on unstructured meshes II: Non-conservative hyperbolic systems[END_REF][START_REF] Franquet | RungeKutta discontinuous Galerkin method for the approximation of Baer and Nunziato type multiphase models[END_REF][START_REF] Wang | High-order CFD methods: current status and perspective[END_REF], in the present work we derive firstly the steady analytical solution of the governing equations (1) with a cylindrical coordinate system, in which one has x = r cos θ, y = r sin θ and u = u θ (r, θ, t)e θ + u r (r, θ, t)e r with u θ and u r being the velocity components in tangential and radial directions, respectively. By definition, one has u r (r, θ, t) = 0 and ∂u θ /∂θ = 0. It can then be easily verified that the two mass conservation equations are automatically satisfied. The momentum conservation equation for the two-phase mixture then becomes

(B.1)

where P , r 0 and α are chosen as P = 0.1ρ 2,0 c 2 2 , r 0 = 0.2L and α = 10 -4 with c 2 = 30 m/s and L = 1 m. With these two conditions, one can obtain the velocity field u θ by solving Equation (B.1). The velocity components u and v in x-and y-directions of the Cartisian coordinate system is expressed as

v(r, θ, t) = + cos θ r ρ ∂p ∂r ,

(B.3) which will be used as the initial velocity fields.

Finally, in order to have an unsteady test-case, we add two constants U = 2 m/s and V = 2 m/s to the previous velocity fields, so that this multiphase vortex is advected in the diagonal direction of the computational domain with a constant speed.


Nana Kofi

Sacha Hodencq

Manar Amayri

Stéphane Ploix

Frédéric Wurtz

Optimal Sizing of Stationary Battery Storage Taking into Account Indirect Flexibility in Tertiary Buildings: Use case of an Electric Vehicle Community

This article proposes a method for taking into account the indirect flexibility (charging of Electric Vehicles) available in buildings for the sizing of stationary battery storage systems (direct flexibility). A Linear Programming approach was applied to data sourced from the Predis-MHI platform (a living lab) such that the day-to-day charging of EVs, as well as the scheduling of the charging and discharging of the proposed battery, were optimized whilst simultaneously dimensioning the battery capacity.

Our results indicate that based on the percentage increase in self-consumption with reference to the base case, it is possible to reduce the battery capacity required by up to 100% as compared to a methodology that doesn't account for the indirect flexibility. Whilst relevant, the sizing approach proposed in this article assumes optimal human behavior, which is usually difficult to achieve. Our proposed approach can be adapted and used for the dimensioning of direct flexibilities for both residential and commercial/public buildings.

INTRODUCTION

The topic of climate change is a cause for immediate action and in recent times has gained a justifiably large response from the press, researchers, and the general populace. The resulting consequence has been a global shift towards sustainable practices in all sectors of the economy. Transportation and electricity generation have been identified as major contributors of Green House Gases (GHGs). Statistical data indicates that transport accounted for 25.8%, energy industries (electricity generation, petroleum refinery, heat production, etc.) 24.1.0%, and the building sector 14.2% of the EU's GHG emissions in 2019 [START_REF] Olivier | Infographics: greenhouse gas emissions in the European Union[END_REF]. Given these figures, these two sectors have seen an enormous shift and uptake of sustainable policies, innovations, and interventions (renewable energy and Emobility).

Additionally, Energy Communities (a group of individuals with a common goal, who actively participate in the energy sector and organize collective and citizen-driven energy actions [START_REF] Caramizaru | Energy communities : an overview of energy and social innovation[END_REF]) have emerged and are a key component in achieving the EU's energy transition goals. These energy communities bring to the table energy services, however, they play an ancillary role of increasing energy awareness and potentially streamlining the behavior of their members (in terms of energy use) in a desirable and much-needed direction [START_REF] Bauwens | What Roles for Energy Cooperatives in the Diffusion of Distributed Generation Technologies?[END_REF].

As renewable energy and e-mobility (electric vehicles, electric bicycles, etc.) adoption increases, they pose a threat to the reliability and robustness of existing grid infrastructure. The Smart Grid concept has been proposed, as a means to mitigate if not alleviate these potential challenges [START_REF] Ipakchi | Grid of the future[END_REF], which will emerge. Tertiary sector buildings offer a perfect paradigm, in that they are normally used during the day (when solar energy is available), and its occupants who own personal transportation units would have them on or close to the site. Thus, it is possible to use e-mobility solutions (Electric Vehicles) as an effective means of indirect flexibility. In reality, this is difficult to achieve since in many energy systems (buildings included) humans are expected to interact with the said system as such [START_REF] Wurtz | Smart buildings' integrated in 'smart grids': A key challenge for the energy transition by using physical models and optimization with a 'humanin-the-loop' approach[END_REF], humans are thus an (often overlooked [START_REF] Janda | Buildings don't use energy: People do[END_REF]) important parameter to consider for all phases of a building's life cycle (conception, design, and operation).

In this article, we propose a Linear Programming (LP) approach for sizing stationary battery storage systems (direct flexibilities) taking into account an optimal behavior from EV users (indirect flexibility) in a tertiary sector building set-up (the Predis-MHI platform of G2ELab).

DIRECT AND INDIRECT FLEXIBILITIES IN LITERATURE

Flexibility is a key concept for achieving energy transition goals. This has become particularly feasible and necessary owing to the increased adoption and integration of Distributed Energy Resources (DERs). The French Agency for Environment and Energy Management (Agence de l'Environnement et de la Maîtrise de l'Energie -ADEME) has proposed a plan for 100% Renewable Energy penetration in France by 2050 and estimates that buildings would be capable of offering flexibility services of approximately up to 18% of the peak demand (100 GW) [START_REF] Dubilly | A 100% renewable electricity mix? Analyses and optimisations Testing the boundaries of renewable energy-based electricity development in metropolitan France by[END_REF]. At the scale of the building, and given the increasing growth of prosumers (DERs) in the energy sector, flexibility is key in maximizing Self-Consumption (SC).

In this paper, we propose to classify Flexibility into two distinct groups, direct and indirect flexibility. Direct flexibility can be defined as devices within the global energy system whose actions can be fully automated/controlled using software or a control signal and can significantly affect the overall performance of the global energy system (e.g. heat pumps, batteries, etc.).

Paper n° 1275

CIRED 2022 Workshop 2/5

Conversely, indirect flexibility is defined as actors or subsystems in the global system that cannot be controlled directly with a control signal but whose actions directly influence the performance of the energy system, mainly referring to the human in the loop [START_REF] Twum-Duah | Assessing the Impact of Direct and Indirect Flexibilities on the Self-Consumption of A Tertiary Building: The Case of Predis-MHI, A Smart Tertiary Building[END_REF].

Hodencq et al. [START_REF] Hondencq | Including greenhouse gas emissions and behavioural responses for PV self-sufficient optimal design[END_REF] proposed an approach for sizing PV systems that takes into account the Global Warming Potential (GWP) of the components and the operational carbon impact, positive feedback (reduced demand with time) and the degradation of the individual components of the PV system with self-sufficiency as a constraint (i.e. the energy system must achieve the defined SC). The results showed that small PV systems are effective at reducing the carbon intensity of a building's energy system whilst systems designed to achieve higher self-sufficiency tend to increase the carbon footprint of the said building. The optimal level of self-production to reap maximum benefit from the addition of a solar PV system is, however, dependent on many factors, which include the energy mix of the grid electricity, and the demand profile.

Current approaches for sizing storage systems both for residential and public buildings do not take into account indirect flexibility. Tostado-Veliz et al [START_REF] Tostado-Véliz | A novel methodology for optimal sizing photovoltaic-battery systems in smart homes considering grid outages and demand response[END_REF] proposed a MILP strategy for sizing solar PV systems taking into account grid failures and demand response (using direct flexibilities) with the objective function being to reduce the system cost. Previous works have focused on sizing taking into account financial costs [START_REF] Mulleriyawage | Optimally sizing of battery energy storage capacity by operational optimization of residential PV-Battery systems: An Australian household case study[END_REF], and environmental impacts [START_REF] Mehrtash | An Enhanced Optimal PV and Battery Sizing Model for Zero Energy Buildings Considering Environmental Impacts[END_REF], [START_REF] Hondencq | Including greenhouse gas emissions and behavioural responses for PV self-sufficient optimal design[END_REF], however, similar to other buildingrelated models the human in the loop has been ignored in these solutions. Taking into account the indirect flexibility is bound to yield a smaller system size and financial commitment reducing the building's carbon impact.

USE CASE: THE PREDIS-MHI PLATFORM

Figure 1 Predis-MHI electricity system (represented using OMEGAlpes graphical representation [START_REF] Hodencq | OMEGAlpes, an open-source optimisation model generation tool to support energy stakeholders at district scale[END_REF])

This article considers the Predis-MHI platform (a living lab), which is a platform of the Grenoble Electrical Engineering Laboratory (G2ELab) dedicated to the study of smart, highly energy-efficient buildings and occupant behavior. The platform is a 600m² multi-floor space composed of lecture rooms, a demonstration room, a student's foyer, and offices housed in the Grenoble Energy-Education and Research (GreEn-ER) building [START_REF] Delinchant | GreEn-ER Living Lab -A Green Building with Energy Aware Occupants[END_REF]. The energy system of the platform is composed of a 22 kWp solar PV system, 50kWh lithium-ion battery, lighting, Heating, Ventilation and Air-Conditioning (HVAC), electrical outlets, and 4 EVs charging stations (2 -14 kW chargers and 2 -44 kW added in August 2021).

Figure 1 shows the energy system of the platform.

Data Management

To allow for further improvements of the proposed model, this article follows the Open and Reproducible Use Cases for Energy (ORUCE) principles [START_REF] Hodencq | Open and Reproducible Use Cases for Energy ( ORUCE ) methodology in systems design and operation : a dwelling photovoltaic selfconsumption example[END_REF]. Thus, the data, associated code and notebooks are available online and can be modified and applied to different use cases and scenarios.

Data Analysis

Figure 2 Annual hourly charger utilization vs PV production

Figure 2 shows the charging patterns of the occupants of the building compared to the average annual solar PV production. It can be deduced from Figure 2 that the behavior of EV owners concerning the charging of their EVs is less than ideal with 9 and 10 AM being the favored time (occupants are plugging and charging their vehicles when they arrive). Ultimately, Figure 2 shows that there exists some degree of indirect flexibility that can be used to improve the SC rate.

LP FORMULATION

The optimal sizing of direct flexibilities implies that both the long-term and short-term constraints of the system need to be accounted for. The addition of indirect flexibilities further implies that mid-term constraints have to be introduced into the optimization. In the case of the proposed optimization, the mid-term constraint was to ensure that the energy consumed within a defined period (1 day for the study) is exactly met by the charging schedule proposed by the optimizer whilst respecting the operating hours of the building (8:00 -20:00). For this LP optimization, the objective function was defined as:

𝑜𝑏𝑗𝑒𝑐𝑡𝑖𝑣𝑒 = 𝑚𝑖𝑛[𝐶𝑎𝑝𝑎𝑐𝑖𝑡𝑦 𝑏𝑎𝑡𝑡𝑒𝑟𝑦 ] 1
To ensure that the battery state of charge (SOC) stays within defined operational bounds (20% -100% considered for the study):

SOC Min × 𝐶𝑎𝑝 𝑏𝑎𝑡 ≤ E 𝑏𝑎𝑡 (pd, t) ≤ 𝑆𝑂𝐶 𝑀𝑎𝑥 × 𝐶𝑎𝑝 𝑏𝑎𝑡 2
Where SOC Min × 𝐶𝑎𝑝 𝑏𝑎𝑡 and 𝑆𝑂𝐶 𝑀𝑎𝑥 × 𝐶𝑎𝑝 𝑏𝑎𝑡 refers to the minimum and maximum allowed battery energy respectively and E 𝑏𝑎𝑡 (pd, t) is the electric charge in the battery for the period pd (day) at time step t, expressed as:

E 𝑏𝑎𝑡 (𝑝𝑑, 𝑡) = E 𝑏𝑎𝑡 (𝑝𝑑, 𝑡 -1) × [1 -∆ 𝑠𝑑 ] + [𝑃 𝑏𝑎𝑡 (𝑝𝑑, 𝑡) × 𝜂 𝑏𝑎𝑡 × 𝑡𝑠] 3
Where 𝑃 𝑏𝑎𝑡 (𝑝𝑑, 𝑡) is the charge/discharge power of the battery for period p at time t, ∆ 𝒔𝒅 is the self-discharge rate of the battery and 𝜂 𝑏𝑎𝑡 is the battery efficiency. Additionally, to ensure continuity between the various periods, the battery's SOC at the end of each period was defined as:

E bat (p, T + 1) = E bat (p, T) × [1 -∆ sd ] + [P bat (p, T) × η bat × ts] 4
Where 𝑇 is the last time step in the set {0,1,2 … 𝑇}, and the final battery energy is constrained as defined in equation ( 2). Thus, the starting battery energy for the various periods was defined as:

E 𝑏𝑎𝑡 (pd, t) = { 𝐸 𝑏𝑎𝑡 (𝑝𝑑, 𝑇), 𝑖𝑓 𝑝𝑑 = 0 𝑎𝑛𝑑 𝑡 = 0 E 𝑏𝑎𝑡 (𝑝𝑑 -1, 𝑇 + 1), 𝑖𝑓 𝑝𝑑 > 0 𝑎𝑛𝑑 𝑡 = 0 5

Where E 𝑏𝑎𝑡 (pd, t) is the final energy at the end of the simulation, i.e. the optimizer determines the ideal SOC to initialize the battery. This condition is especially important as the SC increases and the required battery capacity increases in response (at higher SC, the initial battery SOC impacts the dimensioning of the battery). Further, to ensure that the battery charging and discharging power stayed within defined safe operating bounds:

𝑅 𝑑𝑖𝑠𝑐ℎ𝑎𝑟𝑔𝑒 × 𝐶𝑎𝑝 𝑏𝑎𝑡 𝑡𝑠 ≤ P 𝑏𝑎𝑡 (pd, t) ≤ 𝑅 𝑐ℎ𝑎𝑟𝑔𝑒 × 𝐶𝑎𝑝 𝑏𝑎𝑡 𝑡𝑠 6

Where 𝑅 𝑑𝑖𝑠𝑐ℎ𝑎𝑟𝑔𝑒 and 𝑅 𝑐ℎ𝑎𝑟𝑔𝑒 are factors that represent the proportion of the nominal capacity that can be charged or discharged from the battery at any given time step. The charging of EVs was defined mathematically as

𝐸 𝑐ℎ𝑎𝑟𝑔𝑒𝑟𝑋 𝑜𝑟𝑖𝑔𝑖𝑛𝑎𝑙 (𝑝𝑑) = ∑ 𝑃 𝑐ℎ𝑎𝑟𝑔𝑒𝑟𝑋 (𝑝𝑑, 𝑡) × 𝑡𝑠 𝑇 𝑡 7
Where 𝐸 𝑐ℎ𝑎𝑟𝑔𝑒𝑟𝑋 𝑜𝑟𝑖𝑔𝑖𝑛𝑎𝑙 (𝑝𝑑) is the energy originally consumed by charger X (X a set of charger IDs [START_REF] Olivier | Infographics: greenhouse gas emissions in the European Union[END_REF][START_REF] Caramizaru | Energy communities : an overview of energy and social innovation[END_REF][START_REF] Bauwens | What Roles for Energy Cooperatives in the Diffusion of Distributed Generation Technologies?[END_REF][START_REF] Ipakchi | Grid of the future[END_REF]) during period pd for the optimization horizon and 𝑃 𝑐ℎ𝑎𝑟𝑔𝑒𝑟𝑋 (𝑝𝑑, 𝑡) is the power of charger X for period pd at time t. Additionally, to respect the time constraints imposed by the operating hours of the building:

∑ 𝑃 𝑐ℎ𝑎𝑟𝑔𝑒𝑟 𝑥 (𝑝, 𝑡) 𝑋 𝑥 × 𝑡𝑠 = 0, 𝑖𝑓 𝑡 < 𝑇 𝑜𝑝𝑒𝑛 𝑜𝑟 𝑡 > 𝑇 𝑐𝑙𝑜𝑠𝑒 8
Where 𝑇 𝑜𝑝𝑒𝑛 and 𝑇 𝑐𝑙𝑜𝑠𝑒 refer to the opening and closing times of the GreEN-ER building respectively. Equation 9below prevents mitigates the use of the grid for charging the battery.

𝑃 𝑖𝑚𝑝𝑜𝑟𝑡 (𝑝𝑑, 𝑡) ≤ 𝑃 𝑙𝑜𝑎𝑑 (𝑝𝑑, 𝑡) + ∑ 𝑃 𝑐ℎ𝑎𝑟𝑔𝑒𝑟 𝑥 (𝑝𝑑, 𝑡)

𝑋 𝑥 9
Where 𝑃 𝑖𝑚𝑝𝑜𝑟𝑡 (𝑝, 𝑡) and 𝑃 𝑙𝑜𝑎𝑑 (𝑝𝑑, 𝑡) are the grid import power and demand power of the platform for period pd at time t respectively. The grid import and export are further constrained by ensuring that the SC for the evaluation period is equal to the defined rate, and is expressed as [START_REF] Luz | Coordinating shiftable loads for collective photovoltaic self-consumption: A multiagent approach[END_REF]:

∑ ∑ 𝑃 𝑙𝑜𝑎𝑑 (𝑝𝑑, 𝑡) -𝑃 𝑖𝑚𝑝𝑜𝑟𝑡 𝑇 𝑡 (𝑝𝑑, 𝑡) 𝑃𝐷 𝑝𝑑 ∑ ∑ 𝑃 𝑝𝑣 (𝑝, 𝑡) 𝑇 𝑡 𝑃𝐷 𝑝𝑑 = 𝑅 𝑠𝑐 10
Where 𝑅 𝑠𝑐 , is an input to the optimization, and acts as a constraint for the LP problem (i.e. it is the desired SC rate for which the direct flexibility should be dimensioned).In addition, the power balance is ensured by: 𝑃 𝑃𝑉 (𝑝𝑑, 𝑡) -𝑃 𝑏𝑎𝑡 (𝑝𝑑, 𝑡) + 𝑃 𝑔𝑟𝑖𝑑 (𝑝𝑑, 𝑡) + 𝑃 𝑡𝑜𝑡 (𝑝𝑑, 𝑡) = 0 11

Given:

𝑃 𝑡𝑜𝑡 (𝑝𝑑, 𝑡) = [𝑃 𝑙𝑜𝑎𝑑 (𝑝𝑑, 𝑡) + ∑ 𝑃 𝑐ℎ𝑎𝑟𝑔𝑒𝑟 𝑥 (𝑝𝑑, 𝑡) 𝑋 𝑥 ] 12
𝑃 𝑔𝑟𝑖𝑑 = 𝑃 𝑖𝑚𝑝𝑜𝑟𝑡 (𝑝𝑑, 𝑡) + 𝑃 𝑒𝑥𝑝𝑜𝑟𝑡 (𝑝𝑑, 𝑡) 13

RESULTS

Using the optimization problem as defined above (solved using Gurobi [START_REF] Optimization | Gurobi Optimizer Reference Manual[END_REF]), a parametric study was conducted to determine the effect of indirect flexibility (the human actors, EV users in this case that do plug-in their vehicles

Paper n° 1275

CIRED 2022 Workshop 4/5

at the right moment) on the sizing of battery storage. Figure 4 shows the optimal battery size for a defined annual SC rate with and without the indirect flexibility and the percentage reduction in the required storage capacity with the direct flexibility only approach as reference.

From figure 4, it can be observed that an SC rate of 54% can be attained with no need for battery storage, however, in the case of the direct flexibility only scenario, a 6.00kWh battery is required to achieve the same 54% SC rate. Thus, it can be deduced that there exists a potential of increasing the SC rate by approximately 8% (with reference to the natural SC) using only the indirect flexibility.

Figure 4 Comparison of battery sizing considering only direct flexibility and both direct and indirect flexibility

Figure 4 also shows the potential savings in terms of both financial and environmental impact (relative to the direct only approach), with the reduction in required battery capacity to achieve an SC rate of 70% being approximately 30.20%.

Figure 5 Comparison of charging frequency for direct only and direct & indirect strategies (SC =70%)

For the platform, it was found that the optimal SC with respect to the results in figure 4 (i.e. the SC for which the gain of using indirect flexibility is maximized) was approximately 70%. Considering this SC rate, figure 5 shows a comparison of the hourly charging frequency for the reference case, the direct only, and direct and indirect flexibilities using the battery capacities determined for an SC rate of 70%. From figure 5, by using the indirect flexibility, we were able not only to achieve a lower rate of SC but also to maximize the use of the indirect flexibility (i.e. more charging occurring at peak sun hours).

To assess the potential effort required to move from the sub-optimal (reference case) to the optimal case (direct only and direct and indirect approaches), an evaluation indicator is proposed such that:

𝑒𝑓𝑓𝑜𝑟𝑡 = ∑ |𝑓𝑟𝑒𝑞 𝑜𝑝𝑡𝑖𝑚𝑎𝑙 (𝑡) -𝑓𝑟𝑒𝑞 𝑟𝑒𝑓𝑒𝑟𝑒𝑛𝑐𝑒 (𝑡)| 𝑡 ∑ 𝑓𝑟𝑒𝑞 𝑟𝑒𝑓𝑒𝑟𝑒𝑛𝑐𝑒 (𝑡) 𝑡 14
Where 𝑓𝑟𝑒𝑞 𝑜𝑝𝑡𝑖𝑚𝑎𝑙 (𝑡) and 𝑓𝑟𝑒𝑞 𝑟𝑒𝑓𝑒𝑟𝑒𝑛𝑐𝑒 (𝑡) are the hourly charging frequency (i.e. number of times a vehicle was charging for said hour during the evaluation period) of charging in the optimal and reference cases respectively.

Where a value of zero implies zero to minimal effort would be required, whilst a value of one would imply the greatest level of effort would be required to reach the optimal case. Thus, the effort required to reach the direct only and direct and indirect optimals were calculated to be 0.51 and 0.52 respectively, indicating that a lower level of effort is required from EV users to reach the direct only scenario optimal relative to the effort level required to attain the direct & indirect scenario optimals.

Figure 6 Comparison of monthly SC rates for different battery capacities

Additionally, the growth of the battery capacity with increasing SC was observed to be linear in both cases until approximately 76%, beyond which the battery capacity required grows exponentially and consequently diminishes the gain of using the indirect flexibility. We evaluate the impact of the battery on the monthly SC rate of the platform. From Figure 6, which shows the monthly SC rates for different battery capacities (using both the direct and indirect flexibilities), it can be seen that the impact of the significantly larger battery is greatest during the spring and autumn seasons when PV production is relatively high (more than in winter but lower than in summer). For the summer periods, however, the significantly larger battery capacity has little to no impact on the SC rates, indicating that, the production is in excess of the load and the batteries are underutilized for these months.

CONCLUSION

The indirect flexibility potential of a building (and in extension its occupants' behaviors) is an important parameter that should be taken into account for the dimensioning of energy storage systems (both for tertiary sector and residential buildings). Whilst it is difficult to ensure optimal behavior from occupants, methods and tools exist and are constantly under development to influence occupant behaviors towards an optimal point. The results of this study indicate that indirect flexibility alone can improve a building's SC, however, by taking into account the indirect flexibility available, it is possible to achieve the same SC rates as a sizing, which only takes into account the direct flexibility using a significantly smaller battery capacity. Indirect flexibilities are a key aspect of the energy system of a building and are indispensable for achieving the energy transition and climate mitigation goals.

Our work did not take into consideration the long-term effects of behavioral change of occupants (i.e. reduced demand due to increased energy awareness) or the growing adoption of EVs for private transportation needs. Additionally, we assumed optimal behavior from the EV users in this sizing approach, which in reality might be difficult to achieve. Further work is required on mechanisms and tools to influence human actors towards a globally optimum behavior. Energy communities, nudging and the use of feedback indicators present a suitable use case for such a study.

SUPPLEMENTARY MATERIALS

A Gitlab repository including the code and data sources is available: https://gricad-gitlab.univ-grenoblealpes.fr/NanaKofi/battery_sizing

Figure 3

 3 Figure 3 Annual hourly demand and PV productionAdditionally, comparing the annual average demand of the Predis-MHI platform (without EV charging) to the average annual PV production, see Figure3, we show that there exists a significant potential to increase the SC of the platform using the indirect flexibility (EV charging). The natural SC (SC of a building without any technical

AKNOWLEDGEMENT

This work was funded and supported by ADEME RETHINE, CARNOT Interebat and ORCEE and the interdisciplinary project IDEX Eco-Sesa (see: https://ecosesa.univ-grenoble-alpes.fr), itself supported by ANR project ANR-15-IDEX-02.