Nore

J Varela

S Brun

B Dubrulle

C Nore

Role of boundary conditions in helicoidal flow collimation: Consequences for the von Kármán sodium dynamo experiment

Keywords: number(s): 47.20.Ky, 47.65.-d, 47.27.Cn, 52.72.+v

come

Role of boundary conditions in helicoidal flow collimation: Consequences for the von Kármán sodium dynamo experiment

I. INTRODUCTION

The dynamo effect is the main driver of magnetic field in celestial bodies through the transformation of mechanical energy in magnetic energy [1][START_REF] Moffatt | Magnetic Field Generation in Electrically Conducting Fluid[END_REF][START_REF] Brun | [END_REF]. Given their huge dimensions and the extreme conditions prevailing in their environment, natural dynamos are characterized by control parameters that are presently out-of-reach of analytical or numerical grasp. To make progress, several liquid-metal dynamo experiments have been built [4][5][6], allowing observation and study of several regimes of astrophysical or geophysical interest with interesting open questions regarding the dynamo mechanisms. When the level of fluctuations is low, like in the Riga or Karlsruhe dynamo, the dynamo mechanism appears well described by a "laminar" mechanism, by which the mean velocity flow stretchs and folds the magnetic field lines, resulting in magnetic field growth [START_REF] Moffatt | Magnetic Field Generation in Electrically Conducting Fluid[END_REF]. The observed magnetic field is then similar to that predicted by kinematic computations. When the level of fluctuations is large, like in the VKS dynamo, the observed dynamo magnetic field is dipolar aligned with the symmetry axis of the setup, in contradiction with the "laminar" mechanism that predicts an equatorial dipole [7]. It has therefore been suggested that the dominant mechanism in VKS is, like in many celestial bodies, "turbulent" [8][9][10], with a magnetic field generation occurring either via the interaction of the differential rotation and the nonaxisymmetric velocity perturbations [8], or via self-interaction of the helical perturbations [11]. However, many features of the VKS dynamo are left unexplained. The VKS experiment generates turbulence through the counter-rotation of two impellers (i.e., disks fitted with blades). Both the geometry and the composition of the impellers influence the dynamo: (i) when the blades curvature is decreased, the dynamo threshold increases [12]; (ii) when the blades and the disks are made of soft-iron, with large magnetic permeability, the dynamo threshold is at least two to three times smaller than when the blades and/or the disks are made of copper, with large electrical conductivity, or with stainless steel [13]; (iii) when the ferromagnetic impellers are magnetized, the dynamo threshold decreases [12]. This shows that a better understanding of the VKS dynamo requires investigation of the complex interaction of the large-scale flow, generated by the impeller, with magnetic field and current at the impeller. Previous investigations of the role of impeller material on the dynamo mechanism in VKS were restricted to the role of the boundary conditions onto the magnetic field growth with prescribed velocity [10,14,15].

In this paper, we study the influence of a prescribed magnetic field (presumably issued from the disk magnetization) on the velocity fluctuations. We use magnetohydrodynamic (MHD) numerical simulations in a simplified geometry mimicking the flow structure in the vicinity of the impeller. We show that the interplay between the magnetic and velocity fields results in a strong enhancement of the helicity only when ferromagnetic boundary conditions are used.

II. NUMERICAL MODEL

We use the PLUTO code with a viscous and resistive MHD single fluid model in 3D Cartesian coordinates [16]. The simulation domain is plotted in Fig. 1(a) and focused on the near impeller region and the flow in between two blades, with X, Y , and Z directions corresponding to local azimuthal (toroidal), radial, and vertical (poloidal) directions. For simplicity, we consider straight blades instead of curved blades and walls without thickness. The gray surfaces on Fig. 1(a) represent the blades (at X = 0 and X = 2), the impeller disk (at Z = 0) and the cylinder outer wall (at Y = 4). The influence of the blade's geometry is taken into account via the velocity boundary condition, through , the ratio of the poloidal to toroidal mean velocity that varies from 0.9 to 0.4 as the blade's curvature is changed from 72 o to -72 o [17]. In the blades and the impeller disk we impose either perfect ferromagnetic ( B × n = 0) or perfect conductor ( B • n = 0) boundary conditions ( n is the surface unitary vector), null velocity, and constant slope for the density (ρ) and pressure (p). For the wall at Y = 4 and at the other boundaries, the magnetic field is fixed to 10 -3 T and oriented in the azimuthal X direction, mimicking an azimuthal disk magnetization as observed in the VKS experiment [18]. The value of 10 -3 T has been chosen to match the order of magnitude of the remnant magnetic field observed in the impeller, after a dynamo has been switched off. At the top boundary, the velocity is fixed to V = (10, 0, -10 ) m/s, mimicking the impinging velocity field due to Ekman pumping towards the impeller. The density is fixed to 931 kg/m 3 in the left wall outside the blade (X = 0) and has a constant slope in the rest. The pressure is calculated as p = ρv 2 s /γ with γ = 5/3 the specific heat ratio and v s = 250 m/s the sound speed. The v s value is one order smaller than the real sound speed in liquid sodium to keep a time step large enough for the simulation to remain tractable. The consequence is a small enhancement of the compressible properties of the flow (subsonic low Mach number flow or pseudoincompressibility regime); nevertheless the impact on the simulations is small and the incompressible nature of the liquid sodium flow is preserved. The numbers of grid points are typically 128 in the X and Z directions and 256 in the Y direction for the simulations with kinetic Reynolds number Re n = ρV L/η = 200, with L = 1 m and η the dynamic viscosity. For the simulations with Re n = 1000 we double the resolution in each dimension. We chose to illustrate our results with the Re n = 200 cases, which are more didactic to explain the physical mechanisms at play than the Re n = 1000 cases, since in the latter the system enters a chaotic regime due to the precession of the whirl vortex, leading to a more complex regime to analyze. However, we demonstrate in the text and Table I that the trends in both configurations are the same. The robustness of the model was further tested for different values of magnetic diffusivity (λ) as showed in Fig. 1(c) for a range of simulations with magnetic Reynolds number Re m = V L/λ taking values between 10 to 10 4 . The magnetic energy decreases when the magnetic diffusivity increases; in the simulation with Re m = 10 no enhancement of the magnetic field is observed, while for a simulation with Re m = 100 the magnetic energy is around six times smaller than the Re m = 10 4 case. The effective magnetic Reynolds number of the numerical magnetic diffusion due to the model resolution corresponds to Re m ≈ 6 × 10 3 in the simulation with Re n = 200. The simulations we show in this paper are performed using the numerical magnetic diffusion of the code for simplicity. We perform hydrodynamic (HD) simulation with no magnetic field and MHD simulations, with perfect ferromagnetic and perfect conductor boundary conditions at the blades and impeller disk.

III. VORTEX LOCATION AND AMPLITUDE

We have run several HD simulations with increasing Reynolds number Re n and ranging from 0.6 to 1. In all cases, we observe that the impinging of the velocity at the impeller generates a radial helicoidal vortex, with a kinetic helicity increasing from 50 to 330 ms -2 when Re n is varied from 10 to 1000, meaning that the vortex becomes more and more focused. A visualization of the resulting helicoidal vortex is displayed in Fig. 2 This intensification of the vortex and its displacement closer to the wall can potentially make it more sensitive to any current circulating within the blades. These results are compatible with HD simulations performed by other authors [19]. To study this effect, we thus introduce a background magnetic field to investigate this interplay between the vortex flow and the magnetic field.

IV. INFLUENCE OF BOUNDARY CONDITIONS

We now fix Re n = 200, = 0.8 (corresponding to TM 73 impeller rotating in the unscooping direction in the VKS experiment) and apply a large-scale magnetic field in the azimuthal (X) direction at the boundaries, with ferromagnetic (denoted as Bvecn) or perfect conductor (denoted as Bdotn) boundary conditions at the impeller. In both cases, there is a generation of a magnetic field in the radial (Y ) direction. The magnetic field lines are collimated by the vortex leading to a local enhancement of the magnetic field. In the (Bvecn) case [Fig. 2(a)] the magnetic field is almost five times larger than in the (Bdotn) case [Fig. 2(b)] pointing out the essential role played by the magnetic field boundary conditions in making the interaction between the flow and the magnetic field more efficient. An explanation of this can be found by investigating the geometry of the current streamlines (Fig. 3). The currents in the ferromagnetic walls are parallel to the surface (gray arrows) and are shorted out with the current lines inside the fluid, as can be observed for the current lines plotted nearby the upstream blade and the impeller base [black lines, Fig. 3(a)], not connected with the wall and creating a current whirl. This short-circuit avoids the transfer of magnetic energy from the fluid to the wall, resulting in a tighter collimation and enhanced magnetic field amplitude overall. For the conducting walls, the scenario is the opposite [Fig. 3(b)]: the currents are perpendicular to the blades and the impeller base, connected with the current lines inside the fluid. Therefore, the fluid accumulates a smaller amount of magnetic energy resulting in a weaker magnetic field collimation and enhancement. To quantify this effect, we define the magnetic energy as a volume average in a region nearby the whirl [ME] = MEd x/ d x, provided in Fig. 4(d). The amount of magnetic energy in the ferromagnetic wall case is seven times the one computed in the conducting wall case. In the sequel, we keep this notation for volume averaged quantities.

V. LORENTZ FORCE FEEDBACK ON VORTEX FLOW

The magnetic field reaches locally values of 10 -2 T (100 gauss) but there is no discernible influence on the mean-flow helicity, as can be seen in Fig. 4 smaller than the kinetic helicity. As a result, the total helicity

(H e T = [ B • J /ρ -v • ω]
) is dominated by the kinetic term, so that it is not influenced by the magnetic field, whatever the boundary conditions [Fig. 4(e)]. In contrast, the helicity of the fluctuations, defined as He

f = [ B • J /ρ -v • ω ],
where the denotes the fluctuating part with respect to the time average, is sensitive to the boundary conditions [Figs. 4(c) and 4(f)]. Splitting the fluctuating helicity into a current J H and a kinetic KH parts, we see in Fig. 4(f) that, in the Bdotn case, J H KH , while in the Bvecn case, J H dominates, showing that the magnetic field is strong enough to perturb the velocity fluctuations.

To further quantify the variations of the fluctuating helicity for our five different setups (HD reference run, ferromagnetic, and conductor cases for both Re n number values), we have computed the nine components of the helicity tensor h ij (data not shown), defined as h ij = ikn u k ∂ j u n . Due to the moderate Reynolds number in the Re n = 200 simulations, the fluctuations are not very strong (of the order of 1% of the mean velocity), so that the representative amplitudes of these quantities (maximum of the order of 10 -3 m/s 2 ) are several orders smaller than in the Re n = 1000 simulations. In order to analyze what a hypothetical α 2 dynamo loop based on the components of the helicity tensor (e.g., B y → h yy ,h zy B x → h xx ,h zx B y ) would yield, one can define a gain factor G using the ratio between the time averaged components of the dynamo loop for the ferromagnetic and conductor cases (e. largest dynamo loop enhancement. In Table I, we show first for the Re n = 200 case that the gain factor is for most components greater than 1. This result confirms that the simulations with ferromagnetic boundary conditions are more efficient to collimate the flow and to induce the magnetic field, a trend that is further reinforced for the simulations with Re n = 1000 where all components of the gain factor are now greater than 1.

Figures 5(a) and 5(b) show that regions of alignment or antialignment of the magnetic and the velocity fields are associated with maxima of the magnetic energy. To clarify the magnetic energy transfer from the fluid to the boundary we calculate the Poynting vector on the upstream blade, defined as P y = 1 μ 0 ( E × B), with E the electric field. The Poynting vector in the Bvecn case is parallel to the surface so there is not a net flux of magnetic energy leaving the system, while in the Bdotn simulation the Poynting vector is not parallel to the surface so there is a net flux of magnetic energy. As a consequence, the magnetic energy in the Bvecn case remains in the fluid while in the Bdotn simulation it can be transferred to the wall.

VI. DISCUSSION

We have shown that the impinging of the large-scale flow upon the impeller generates a coherent radial vortex between the blades. This flow collimates any existing magnetic field that in return modifies locally the velocity fluctuations and the helicity tensor. From Table 1, we see that the products of the dynamo loop h yy h zx , h yy h zy , and h zy h zx increase when changing the wall properties from conducting to ferromagnetic material in the Re n = 200 simulations, a trend that is even more robust in the nonstationary Re n = 1000 cases where all the products are larger for a ferromagnetic material. Within the first-order smoothing approximation, the α tensor is defined as α = τ h [11], where τ is the correlation time. This effect therefore leads to a net increase of the efficiency of the α 2 dynamo mechanism, e.g., B y → h yy ,h zy B x → h xx ,h zx B y in the case of ferromagnetic boundary conditions as soon as the disks are magnetized. This effect is stronger when the vortex is closer to the wall, which is piloted by the parameter . This effect may therefore explain why changing the blade's orientation and the material of the impeller has an impact on the dynamo threshold and validates the use of larger values of α when ferromagnetic blades are used and the disks are magnetized [10,14]. It has often been said that the VKS dynamo is just a rotating magnet because of the essential role of the soft iron. If this were the case, the dynamo properties would not be changed by the flow geometry, as it is observed. In the present work, we show that the soft iron may actually enhance the dynamo cycle by retroaction of the magnetic field onto the velocity fluctuations responsible for the α effect, effect mediated by an helicoidal radial vortex that collimates the magnetic field. Even though the present scenario was elaborated at moderate values of the kinetic Reynolds number, we believe that it remains valid for the conditions realized in the experiments. The next step is to further increase the degree of turbulence and to push toward even more realistic values, such as, for instance, lower magnetic Prandtl number. From an experimental point of view, it would be interesting to seek alternative ways to reinforce the α effect between the blades just by mechanics, by enhancing flow collimation rather than through the use of material favoring strong localized magnetic fields.

FIG. 1 .

 1 FIG. 1. (Color online) (a) Simulation domain for a portion of the flow in between two blades: X, Y , and Z directions correspond to local azimuthal, radial, and vertical direction with X ∈ [0,2], Y ∈ [0,4], Z ∈ [0,2]. (b) Whirl location versus value: X coordinate (dotted line), Z coordinate (dashed line), and module (solid line). (c) Magnetic field for simulations with different magnetic diffusion: Re m = 10 4 (dash-dot line), Re m = 10 3 (dotted line), Re m = 10 2 (dashed line), and Re m = 10 (solid line)

  FIG. 2. (Color online) Magnetic and velocity field stream lines at Re n = 200 and = 0.8 in the perfect ferromagnetic (a) and conductor (b) simulations in the all domain. Isocontour of the magnetic field (0.005 T in the ferromagnetic and 0.001 T in the conductor cases). Magnetic field module in the Y = 2 plane.

FIG. 3 .

 3 FIG. 3. (Color online) Angle between the current and the unitary surface vector in the Bvecn configuration (a) and the Bdotn configuration (b) in a zoomed box. The arrows indicate the orientation of the current on the blades and black lines are current stream lines nearby the upstream blade.

  (a): the kinetic helicity time evolutions ([ v • ω], with ω the vorticity) overlap for the Bvecn, Bdotn, and HD simulations. Moreover, the current helicity ([ B • J /ρ], with J the current density) is five times larger in the Bvecn simulation [Fig. 4(b)], but still two orders of magnitude

FIG. 4 .

 4 FIG. 4. (a) Kinetic helicity. Ferromagnetic case (solid line), conductor case (dashed line), and hydro case (dotted line) (b) Current helicity. (c) Magnetic energy. (d) Total helicity. (e) Helicity of fluctuations. (f) Helicity of fluctuations in Bvecn case: kinetic helicity (solid line), current helicity (dashed line); and Bdotn case: kinetic helicity (dotted line), current helicity (dash-dot-dot line). All values are averaged in a volume localized around the whirl.

  g., G = ( |h ij ||h km | ) ferro /( |h ij ||h km | ) cond ). This allows us to assess what boundary condition configurations lead to the

FIG. 5 .

 5 FIG. 5. (Color online) Angle between the velocity and the magnetic fields for (a) the Bvecn and (b) the Bdotn simulation (plane Y = 1) at t = 1 s. Properties of the Poynting vector in the upstream blade: (c) Bvecn, (d) Bdotn cases. The arrows direction (resp. color) codes its orientation (respectively, its angle with the unitary vector normal to the surface). The color in the blade codes its modulus.

TABLE I .

 I Gain factor (G) between the ferromagnetic and conductor cases for the simulations with Re n = 200 and Re n = 1000. Re n = 200 h yy h zx = 1.90 h yy h xx = 0.01 h zy h xx = 0.01 h zy h zx = 53.50 Re n = 1000 h yy h zx = 2.48 h yy h xx = 2.33 h zy h xx = 1.10h zy h zx = 1.17

1539-3755/2015/92(6)/063015(5) 063015-1 ©2015 American Physical Society

ACKNOWLEDGMENTS

We have received funding by the Labex PALM/P2IO/LaSIPS (VKStars Grant No. 2013-02711), INSU/PNST and ERC PoC Grant No. 640997 Solar Predict. We thank M. Cemeljic, W. Herreman, and the VKS team for fruitful discussions.


H Lochon

F Daude

P Galon

J.-M Hérard

HLLC-type Riemann solver with approximated two-phase contact for the computation of the Baer-Nunziato two-fluid model

The computation of compressible two-phase flows with the Baer-Nunziato model is addressed. Only the convective part of the model that exhibits non-conservative products is considered and the source terms of the model that represent the exchange between phases are neglected. Based on the solver proposed by Tokareva & Toro [42], a new HLLC-type Riemann solver is built. The key idea of this new solver lies in an approximation of the two-phase contact discontinuity of the model. Thus the Riemann invariants of the wave are approximated in the "subsonic" case. A major consequence of this approximation is that the resulting solver can deal with any Equation of State. It also allows to bypass the resolution of a non-linear equation based on those Riemann invariants. We assess the solver and compare it with others on 1D Riemann problems including grid convergence and efficiency studies. The ability of the proposed solver to deal with complex Equations Of State is also investigated. Finally, the different solvers have been compared on challenging 2D test-cases due to the presence of both material interfaces and shock waves: a shock-bubble interaction and underwater explosions. When compared with others, the present solver appears to be accurate, efficient and robust.

Introduction

Different strategies are used to model compressible two-phase flows. Unlike the homogeneous approaches, the two-fluid approaches do not consider a mixture of the two phases with supposed equilibrium between phases. The two-fluid models are built using an averaging procedure of the local conservation laws [START_REF] Ishii | Thermo-fluid dynamic theory of two-phase flows[END_REF][START_REF] Drew | Theory of Multicomponent Fluids[END_REF]. If no phasic equilibrium is assumed, two-pressure two-fluid models are obtained. As shown by Ransom and Hicks [START_REF] Ransom | Hyperbolic two-pressure models for two-phase flow[END_REF], who proposed one of the first two-pressure two-fluid models, this class of models exhibits pleasant mathematical properties, particularly hyperbolic systems. More recently, Baer and Nunziato [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] proposed another model of this class to study granular explosive materials. Then, the model has been used in a wider range of applications and many different closure laws have been proposed to that purpose [START_REF] Glimm | Renormalization group solution of two-phase flow equations for Rayleigh-Taylor mixing[END_REF][START_REF] Saurel | A Multiphase Godunov Method for Compressible Multifluid and Multiphase Flows[END_REF][START_REF] Gavrilyuk | Mathematical and Numerical Modeling of Two-Phase Compressible Flows with Micro-Inertia[END_REF][START_REF] Coquel | Closure laws for a two-fluid two-pressure model[END_REF][START_REF] Gallouët | Numerical Modeling Of Two-Phase Flows Using The Two-Fluid Two-Pressure Approach[END_REF]. The mathematical properties of the Baer-Nunziato model without any source term were first studied by Embid and Baer [START_REF] Embid | Mathematical analysis of a two-phase continuum mixture theory[END_REF]. The main difficulty of the convective part of the model lies in the fact that it cannot be written under a conservative form. However, a detailed analysis of the system shows that the non-conservative products are associated with a linearly degenerate characteristic field due to the particular choice of closure laws of the Baer-Nunziato model. Thus the non-conservative products are not active across shock waves resulting in unique jump conditions. Many Finite Volume schemes have been proposed to tackle the numerical resolution of the Baer-Nunziato model, or similar forms, with different treatments of the non-conservative terms. Saurel and Abgrall [START_REF] Saurel | A Multiphase Godunov Method for Compressible Multifluid and Multiphase Flows[END_REF] proposed to use a "free streaming" condition to discretize the non-conservative terms together with Rusanov and HLL schemes. Complete Riemann solvers that take into account all the characteristic waves of the model have then been proposed. Schwendeman et al. [START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF] and later Deledicque and Papalexandris [START_REF] Deledicque | An exact Riemann solver for compressible two-phase flow models containing non-conservative products[END_REF] built Godunov-type schemes to that purpose. For the treatment of the non-conservative terms, Schwendeman et al. [START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF] proposed the "thin layer approximation": the associated linearly degenerate wave is approximated by a thin layer across which the solution is assumed to be smooth. This approximation has also been used in the HLLC-type solver proposed by Tokareva and Toro [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF]. Simpler but incomplete Riemann solvers have been proposed based on HLLtype solver [START_REF] Ambroso | A Godunov-type method for the seven-equation model of compressible two-phase flow[END_REF] and HLLC-type solvers [START_REF] Liang | Solving seven-equation model for compressible two-phase flow using multiple GPUs[END_REF][START_REF] Furfaro | A simple HLLC-type Riemann solver for compressible non-equilibrium two-phase flows[END_REF]. The so-called VFRoe-ncv scheme used by Gallouët et al. [START_REF] Gallouët | Numerical Modeling Of Two-Phase Flows Using The Two-Fluid Two-Pressure Approach[END_REF] is another attempt to tackle the non-conservative framework. Extension of the Osher scheme to non-conservative systems has also been applied to the Baer-Nunziato model [START_REF] Dumbser | A Simple Extension of the Osher Riemann Solver to Non-conservative Hyperbolic Systems[END_REF]. Finally, relaxation approaches have been applied to both the isentropic and the complete Baer-Nunziato model [START_REF] Ambroso | Relaxation and numerical approximation of a two-fluid two-pressure diphasic model[END_REF][START_REF] Coquel | A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model[END_REF][START_REF] Coquel | A Positive and Entropy-Satisfying Finite Volume Scheme for the Baer-Nunziato Model[END_REF].

In this paper, we propose a new approximation of the linearly degenerate wave associated with the non-conservative terms. This approximation is based on a physically-relevant simplification of the Riemann invariants of the corresponding wave in the "subsonic" case. We propose to use it in a HLLCtype Riemann solver based on the work of Tokareva and Toro [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF]. It results in a complete Riemann solver that can deal with any Equation Of State in a straightforward way thanks to the proposed approximation. Another important feature of the present solver is that no iterative method to solve a non-linear equation is required unlike the original HLLC-type solver. Indeed, the intermediate states of the Riemann problem are approximated as a solution of a linear system. Regarding the integration of the non-conservative terms, our new approximation is in agreement with the "thin layer approximation", therefore we treat them as in [START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF][START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF]. We compare this new solver with the Rusanov scheme and the HLLC-type solver [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF] on both 1D and 2D test-cases. A detailed analysis of the accuracy of the solver is carried out, based on several grid convergence studies on 1D Riemann problems. We also pay a particular attention to the case of vanishing phases as it is involved in a lot of practical configurations and may raise robustness problems. The proposed scheme provides similar results as the original HLLC-type solver in most numerical test-cases. Moreover it appears to be very robust and computationally efficient.

Thus the paper is organized as follows. Section 2 introduces the governing equations of the Baer-Nunziato model. We also recall the main mathematical properties of the model in order to describe the general solution of the Riemann problem, particularly the two subsonic configurations. Section 3 is devoted to the numerical scheme. We first present the approximation of the two-phase contact which is the keystone of the scheme. Then we detailed the HLLC-type Riemann solver based on this approximation and we provide the complete numerical flux for the associated Finite Volume scheme. In section 4, we assess the proposed numerical scheme. Seven 1D Riemann problems have been carefully chosen to compare it with the Rusanov scheme and the HLLC-type scheme proposed by Tokareva and Toro [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF]. Grid convergence studies are carried out on half of the test-cases in order to provide a deeper analysis. Besides, we also assess the ability of the solver to deal with complex Equations of State. Finally, three 2D test-cases are studied: a shock-bubble interaction and two underwater explosions. They both include material interfaces and shock waves (with very high pressure differences in the case of underwater explosions) so they are particularly suited to investigate both the accuracy and the robustness of the schemes.

2 The Baer-Nunziato model

Governing equations

The governing equations of the Baer-Nunziato model [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] without source terms are the following (when restricting to the one-dimensional case):

       ∂ t (α 1 ) + u I ∂ x (α 1 ) = 0 ∂ t (α k ρ k ) + ∂ x (α k ρ k u k ) = 0 ∂ t (α k ρ k u k ) + ∂ x α k ρ k u 2 k + α k p k + (-1) k p I ∂ x (α 1 ) = 0 ∂ t (α k ρ k e k ) + ∂ x (α k ρ k e k u k + α k p k u k ) + (-1) k p I u I ∂ x (α 1 ) = 0 (1)
The subscript k is the index of the phase (k = 1, 2). Thus, α k which is the volume fraction of phase k, lies in ]0, 1[ and complies with the following relation:

α 1 + α 2 = 1
ρ k , u k and p k respectively refer to the mean density, the mean velocity and the mean pressure of phase k. The specific total energy e k is given by:

e k = ε k + 1 2 u 2 k
where ε k is the specific internal energy of phase k, linked to the density and pressure by an Equation Of State (EOS):

ε k = ε k (ρ k , p k ) (2)
Thanks to the EOS, the celerity c k of acoustics waves in phase k is defined using the following:

ρ k c 2 k = (∂ p k ε k ) -1 p k ρ k -ρ k (∂ ρ k ε k ) (3)
The specific entropy s k of phase k complies with:

c 2 k (∂ p k s k ) + (∂ ρ k s k ) = 0 (4)
and the temperature T k is defined as:

1 T k = (∂ p k ε k ) -1 (∂ p k s k ) (5)
The closure laws for the interfacial velocity and pressure (u I , p I ) used by Baer and Nunziato [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] are the following:

(u I , p I ) = (u 2 , p 1 ) (6)
Many other closure laws have been proposed in the literature since then. For example, the approach proposed by Coquel et al. [START_REF] Coquel | Closure laws for a two-fluid two-pressure model[END_REF] and used in [START_REF] Gallouët | Numerical Modeling Of Two-Phase Flows Using The Two-Fluid Two-Pressure Approach[END_REF] is based on two requirements: the enforcement of an entropy inequality and a linearly degenerate field associated with λ = u I . The following closure laws comply with these two conditions:

u I = au 1 + (1 -a) u 2 , a = χα 1 ρ 1 χα 1 ρ 1 + (1 -χ)α 2 ρ 2 p I = bp 1 + (1 -b) p 2 , b = (1 -a)T 2 aT 1 + (1 -a)T 2 with χ ∈ [0, 1]
. We can notice that the case when χ = 0 (or χ = 1) corresponds to the closure laws of Baer and Nunziato [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] and the case when χ = 1 2 to the closure laws proposed in [START_REF] Coquel | Closure laws for a two-fluid two-pressure model[END_REF][START_REF] Gallouët | Numerical Modeling Of Two-Phase Flows Using The Two-Fluid Two-Pressure Approach[END_REF]. Those three choices (χ = 0, 1 2 , 1) have been compared for steam-water transient flows in [START_REF] Lochon | Comparison of two-fluid models on steam-water transients[END_REF]. Other approaches can also be found in the literature [START_REF] Glimm | Renormalization group solution of two-phase flow equations for Rayleigh-Taylor mixing[END_REF][START_REF] Jin | Entropy of averaging for compressible two-pressure two-phase flow models[END_REF][START_REF] Saurel | A Multiphase Godunov Method for Compressible Multifluid and Multiphase Flows[END_REF][START_REF] Chinnayya | Modelling detonation waves in heterogeneous energetic materials[END_REF] and a wider class of models is introduced in appendix A of [START_REF] Lochon | Comparison of two-fluid models on steam-water transients[END_REF]. In the sequel, we use the closure laws proposed by Baer and Nunziato, detailed in equation [START_REF] Batten | On the Choice of Wavespeeds for the HLLC Riemann Solver[END_REF].

System (1) may be recast in the following vectorial form:

∂ t (U) + ∂ x (F(U)) + H(U)∂ x (α 1 ) = 0 (7)
with

U =     α 1 α k ρ k α k ρ k u k α k ρ k e k     , F(U) =     0 α k ρ k u k α k ρ k u 2 k + α k p k α k ρ k e k u k + α k p k u k     and H(U) =     u I 0 (-1) k p I (-1) k p I u I    

Properties and Riemann Problem

The homogeneous part of the Baer-Nunziato model has been first studied by Embid and Baer [START_REF] Embid | Mathematical analysis of a two-phase continuum mixture theory[END_REF].

We first recall its main mathematical properties and then tackle the Riemann problem:

Property 1 (Hyperbolicity and structure of waves). System (1) is hyperbolic, but not strictly hyperbolic. Indeed, it admits the following seven eigenvalues:

λ 1 = u 1 -c 1 , λ 2 = u 1 , λ 3 = u 1 + c 1 , λ 4 = u 2 -c 2 , λ 5,6 = u 2 , λ 7 = u 2 + c 2 (8)
and associated righteigenvectors span the whole space

R 7 if |u 1 -u 2 | = c 1 .
Fields associated with eigenvalues λ 2,5,6 are linearly degenerate (LD) whereas fields associated with eigenvalues λ 1,4,3,7 are genuinely non linear (GNL).

Property 2 (Riemann invariants through LD waves).

• The Riemann invariants of the LD 2-wave are {α 1 , u 1 , p 1 , ρ 2 , u 2 , p 2 }. Thus the volume fraction is constant across the wave and the 2-wave has the classical structure of a monophasic contact discontinuity with respect to phase 1.

• The 5 Riemann invariants of the LD 5,6-wave are the following:

I 1 = u 2 I 2 = α 1 ρ 1 (u 1 -u 2 ) I 3 = α 1 p 1 + α 2 p 2 + α 1 ρ 1 (u 1 -u 2 ) 2 (9)
I 4 = ε 1 + p 1 /ρ 1 + (u 1 -u 2 ) 2 /2 I 5 = s 1
The 5,6-wave is the only one involving a jump of volume fraction therefore it will be referred as two-phase contact in the sequel. Although the structure of this wave is complex, some schemes [START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF][START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF] directly use its 5 Riemann invariants when solving the Riemann problem.

Property 3 (Jump conditions through GNL waves). Unique jump conditions hold within each isolated GNL field:

             [α k ] R L = 0 -σ [ρ k ] R L + [ρ k u k ] R L = 0 -σ [ρ k u k ] R L + ρ k u 2 k + p k R L = 0 -σ [ρ k e k ] R L + [ρ k e k u k + p k u k ] R L = 0 (10) where [ϕ] R L = ϕ R -ϕ L
is the jump between the Left and Right states on each side of a shock wave traveling at speed σ.

The fact that the void fraction remains constant across the GNL waves is closely linked to the LD structure of the two-phase contact and leads to classical monophasic jump conditions through GNL waves.

Property 4 (Entropy of the system). Assume that (u I , p I ) comply with the following relation:

k 1 T k (-1) k (p k -p I ) (u k -u I ) = 0 (11)
Then smooth solutions of system (1) are such that:

∂ t k α k ρ k s k + ∂ x k α k ρ k s k u k = 0 (12)
Proof. Smooth solutions of system (1) comply with:

∂ t (α k ρ k s k ) + ∂ x (α k ρ k s k u k ) = 1 T k (p k -p I ) (u k -u I ) ∂ x (α k )
So the equation of the total entropy [START_REF] Crouzet | Validation of a two-fluid model on unsteady liquid-vapor water flows[END_REF] is obtained thanks to the equation ( 11) by summing over the phases.

Thus the enforcement of the entropy inequality directly provides the closure law for the interfacial pressure p I thanks to equation [START_REF] Crouzet | Approximate solutions of the Baer-Nunziato Model[END_REF] once the interfacial velocity u I has been set. If we set u I = u 2 for example, we obtain p I = p 1 which is the closure law proposed by Baer and Nunziato given in equation [START_REF] Batten | On the Choice of Wavespeeds for the HLLC Riemann Solver[END_REF]. Moreover, the same strategy has been used to provide closure laws with three, or more, phases [START_REF] Hérard | A three-phase flow model[END_REF][START_REF] Müller | Closure conditions for non-equilibrium multi-component models[END_REF].

Thanks to all those theoretical properties, we address the Riemann problem for system [START_REF] Chinnayya | Modelling detonation waves in heterogeneous energetic materials[END_REF]:

       ∂ t (U) + ∂ x (F(U)) + H(U)∂ x (α 1 ) = 0 U(t = 0, x) = U L , if x < 0 U R , if x > 0 ( 13
)
As pointed out in [START_REF] Andrianov | The Riemann problem for the Baer-Nunziato two-phase flow model[END_REF] (or in [START_REF] Gallouët | Numerical Modeling Of Two-Phase Flows Using The Two-Fluid Two-Pressure Approach[END_REF] in the case of other closure laws for interfacial quantities), the existence and uniqueness of the solution of such a problem is not guaranteed. Despite this, we can exhibit the general structure of a solution: it consists of several intermediate states connected by the 6 waves of the model. The previous properties provide the connection across each isolated wave, even shock waves thanks to the unique jump conditions across GNL waves detailed in property 3. An important difficulty of this Riemann problem lies in the fact that the wave-ordering is not known, different wave configurations are possible. Two cases are usually considered depending on the relative velocity (u 1 -u 2 ): the case when (u 1 -u 2 ) 2 < c 2 1 holds is called the "subsonic" one and the opposite case, when (u 1 -u 2 ) 2 > c 2 1 , is called "supersonic". The subsonic case is considered to be more physicallyrelevant and we only consider the subsonic case in the sequel like many other authors in the literature [START_REF] Coquel | A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model[END_REF][START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF][START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF]. Readers interested in the supersonic case are referred to [START_REF] Daude | On the computation of the Baer-Nunziato model using ALE formulation with HLL-and HLLC-type solvers towards fluid-structure interactions[END_REF][START_REF] Deledicque | An exact Riemann solver for compressible two-phase flow models containing non-conservative products[END_REF][START_REF] Andrianov | The Riemann problem for the Baer-Nunziato two-phase flow model[END_REF].

In the subsonic case, two wave configurations, presented in figure 1, are possible. We denote by (S k ) L (respectively (S k ) R ) the GNL waves associated with λ = u k -c k (respectively u k + c k ). In the same way, we denote (S k ) M the LD waves of the model. Thus, configuration (a) of figure 1 is the case (S 2 ) M ≤ (S 1 ) M and configuration (b) the opposite case (S 1 ) M < (S 2 ) M . For both configurations, the previous properties provide us the relations between the different intermediate states. Besides, as pointed out previously, the GNL waves have a similar structure as in the monophasic case therefore they don't interact with each other and we don't need to consider separate cases when (S 1 ) K < (S 2 ) K or (S 1 ) K > (S 2 ) K for K = L, R. The notation W k used on figure 1 is defined by:

W k =   ρ k u k p k   (14)
Moreover, the volume fraction α k only changes across the two-phase contact (S 2 ) M as pointed out on the figure. We underline that the LD structure of this wave is crucial in order to obtain unique jump conditions. Indeed, this LD structure implies that the non-conservative products are not active in the other waves, especially the GNL waves (equation [START_REF] Coquel | A Positive and Entropy-Satisfying Finite Volume Scheme for the Baer-Nunziato Model[END_REF]). Thus, those non-conservative products are not "true" non-conservative products (i.e. active in GNL fields) which would have required particular definitions [START_REF] Maso | Definition and weak stability of nonconservative products[END_REF]. As a consequence, shock solutions of the Riemann problems are defined in a nonambiguous way. This is of major importance from a computational point of view: without such properties, numerical grid converged shock solutions may differ due to the inner numerical diffusion of the scheme. To the authors' knowledge, this point has received very little attention in the literature on two-fluid models but striking numerical examples are available in [START_REF] Guillemaud | Modélisation et simulation numérique des écoulements diphasiques par une approche bifluide à deux pressions[END_REF] (pages 121-142).

HLLC-type Riemann solver with approximated two-phase contact

In this section, we introduce a HLLC-type Riemann solver to approximate the solutions of the Baer-Nunziato model based on the scheme proposed by Tokareva and Toro [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF]. The key idea of the present scheme is an approximation of the two-phase contact whereas it is treated using the thin layer approximation proposed by Schwendeman et al. [START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF] in the original HLLC-type solver [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF]. Thus, we first introduce the approximated two-phase contact then detail the intermediate states of the solution of the Riemann problem obtained thanks to the approximation and finally, we define the numerical flux of the solver based on the estimated states.

Approximated two-phase contact

The Riemann invariants of the two-phase contact are recalled in equation [START_REF] Coquel | A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model[END_REF]. We introduce the coefficient γ k defined by the following:

γ k = ρ k c 2 k p k (15)
In the particular case of the ideal gas equation of state (EOS), γ k is equal to the ratio of specific heats γ k . Thanks to this notation, some Riemann invariants can be recast in the following form:

I 3 = α 1 ρ 1 1 γ 1 c 2 1 + (u 1 -u 2 ) 2 + 1 γ 2 α 2 ρ 2 c 2 2 I 4 = ε 1 + 1 γ 1 c 2 1 + 1 2 (u 1 -u 2 ) 2
In the subsonic case, we already assume that (u 1 -u 2 ) 2 < c 2 1 . In order to approximate the two-phase contact, we propose a stronger hypothesis:

(u 1 -u 2 ) 2 1 γ1 c 2
1 . This hypothesis is physically relevant since u 2 1 c 2 1 in a lot of practical configurations in the subsonic case. Moreover, the difference (u 1 -u 2 ) 2 should be even smaller due to the velocity relaxation included in a complete model with convective terms and source terms such as drag effects [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF]. Finally, 1 γ1 c 2 1 can be considered to be in the same order of magnitude as c 2 1 , particularly for most usual gases.

Thanks to the proposed hypothesis, the Riemann invariants of the wave are approximated in the following way:

I 1 = u 2 I 2 = α 1 ρ 1 (u 1 -u 2 ) I 3 ≈ α 1 p 1 + α 2 p 2 I 4 ≈ ε 1 + p 1 /ρ 1 I 5 = s 1
Therefore, the specific internal enthalpy ε 1 + p 1 /ρ 1 of phase 1 seems to be a relevant approximation of I 4 . However, the specific entropy s 1 (p 1 , ρ 1 ) of phase 1 is the Riemann invariant I 5 of the wave. Preserving both is equivalent to the preservation of the primitive variables (ρ 1 , p 1 ). Thus, we define the approximated two-phase contact with the following new set of Riemann invariants:

I 1 = u 2 I 2 = α 1 (u 1 -u 2 ) I 3 = α 1 p 1 + α 2 p 2 (16)
I 4 = p 1 I 5 = ρ 1
If we now focus on the consequences of this approximation in the structure of the solution of the 1D-Riemann problem (figure 1), we can notice that the approximation of the state (W 1 ) * M is directly obtained thanks to the states (W 1 ) * K , K = L, R. For example, if we consider configuration (a), we have:

(W 1 ) * M =   (ρ 1 ) * L (u 1 ) * R (p 1 ) * R  
whereas, in configuration (b), we get:

(W 1 ) * M =   (ρ 1 ) * R (u 1 ) * L (p 1 ) * L  

HLLC-type intermediate states

This section is devoted to the estimations of the intermediate states of the solution of the Riemann problem (figure 1) using the approximation of the two-phase contact. In a HLLC-type approach, the intermediate states are obtained by using the averaged Rankine-Hugoniot relation across the GNL waves (S k ) K (k = 1, 2 and K = L, R). The speeds of those waves need to be estimated a priori on the Left and Right states of the Riemann problem as in the HLL approach [START_REF] Harten | On Upstream Differencing and Godunov-Type Schemes for Hyperbolic Conservation Laws[END_REF]. Based on the work of Batten et al. [START_REF] Batten | On the Choice of Wavespeeds for the HLLC Riemann Solver[END_REF] for the Euler equations, we use the following estimations, as in [START_REF] Daude | On the computation of the Baer-Nunziato model using ALE formulation with HLL-and HLLC-type solvers towards fluid-structure interactions[END_REF]:

(S k ) L = min (u k ) L -(c k ) L , ûk -ĉk (S k ) R = max (u k ) R + (c k ) R , ûk + ĉk , k = 1, 2 (17)
where ûk and ĉk are the Roe-type average of the velocity and the speed of sound. We recall that the wave speeds of the LD waves (S k ) M (k = 1, 2) do not need to be estimated a priori since they are directly given by the velocities of the intermediate states. If we consider the configuration (a) of figure 1, we have:

(S 1 ) M = (u 1 ) * R and (S 2 ) M = (u 2 ) * L .
The averaged Rankine-Hugoniot relations for the mass and momentum equations across the GNL waves provide the following HLLC-type intermediate states:

         (ρ k ) * K = (ρ k ) K (S k ) K -(u k ) K (S k ) K -(u k ) * K (p k ) * K = (p k ) K + (ρ k ) K (S k ) K -(u k ) K (u k ) * K -(u k ) K (18)
for k = 1, 2 and K = L, R. They are similar to the intermediate states of the classical HLLC scheme [START_REF] Toro | Restoration of the contact surface in the HLL-Riemann solver[END_REF] as expected since α 1 is constant across the GNL waves. The intermediate states (W k ) * K can now be determined by using the preservation of the Riemann invariants of the approximated two-phase contact. The existence and uniqueness of those intermediate states are guaranteed by the following: Proposition 1. Assume that the intermediate states of the Riemann problem comply with the preservation of the Riemann invariants of the approximated two-phase contact [START_REF] Deledicque | An exact Riemann solver for compressible two-phase flow models containing non-conservative products[END_REF] and with the HLLC relations [START_REF] Dumbser | A Simple Extension of the Osher Riemann Solver to Non-conservative Hyperbolic Systems[END_REF] where the wave speed estimates [START_REF] Drew | Theory of Multicomponent Fluids[END_REF] are used. Then there exist unique intermediate states (W k ) * K (k = 1, 2 and K = L, M, R) of the Riemann problem for all initial conditions (W k ) K (k = 1, 2 and K = L, R) and for all phasic equations of state (EOS).

Proof. Thanks to the HLLC relations (18), we have (u k ) * K as a linear function of (p k ) * K :

u * kK = u kK + 1 c kK (p * kK -p kK ) with c kK = ρ kK (S kK -u kK )
where we use the short notation ϕ kK instead of (ϕ k ) K (k = 1, 2 and K = L, R) for simplicity. In order to determine the 4 intermediate pressures p * kK , we then need to use the Riemann invariants of the approximated two-phase contact from equation [START_REF] Deledicque | An exact Riemann solver for compressible two-phase flow models containing non-conservative products[END_REF]. Those Riemann invariants are also linear functions of the intermediate pressures and velocities therefore the intermediate pressures can be computed as the solutions of a linear system. We can notice that the pressure p 1 is a Riemann invariant I 4 of the approximated contact so we introduce the notation p * 1 = p * 1L = p * 1R and we end up with the following 3 × 3 linear system:

        0 1 c 2L - 1 c 2R α 1L c 1L - α 1R c 1R - α 1L c 2L α 1R c 2R α 1L -α 1R α 2L -α 2R                   p * 1 p * 2L p * 2R           =           RHS 1 RHS 2 0           with RHS 1 = p 2L c 2L -u 2L - p 2R c 2R + u 2R RHS 2 = α 1L p 1L c 1L -u 1L - p 2L c 2L + u 2L -α 1R p 1R c 1R -u 1R - p 2R c 2R + u 2R
We underline the fact that this system is exactly the same for both wave configurations, (a) or (b) in figure 1. The determinant of the matrix can be easily computed:

∆ = α 1L c 1L - α 1R c 1R α 2R c 2L - α 2L c 2R - (α 1L -α 1R ) 2 c 2L c 2R
It is important to notice that c kR > 0 and c kL < 0 for both phases (k = 1, 2) due to wave speed estimates given in equation [START_REF] Drew | Theory of Multicomponent Fluids[END_REF]. From that remark, we deduce that ∆ > 0 so the system has a unique solution for all initial conditions (W k ) K (k = 1, 2 and K = L, R). Once the intermediate pressures have been obtained, the intermediate states (W k ) * K can be fully determined thanks to equation ( 18). Then we also determine the wave configuration, (a) or (b) in figure 1, by looking at the sign of I 2 .

The existence and uniqueness of the intermediate states of the Riemann for all initial conditions is one of the major consequences of the approximation of the two-phase contact. It is important to recall that, without such an approximation, the connection through the two-phase contact is not possible for any Right and Left states as detailed in Appendix A. This approximation also has an important feature from the computational point of view since it leads to a linear system in order to determine the intermediate pressures whereas the original HLLC solver [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF] needs to use an iterative procedure to solve a non-linear system. Finally, another major consequence of that approximation is that the analytical expression of the entropy s 1 , which is an exact Riemann invariant of the wave (I 5 ), is no longer required in the computation so the scheme can deal with any EOS.

Numerical flux

Thanks to the approximate solution of the Riemann problem proposed previously, we can introduce a Finite-Volume scheme based on this solution. This scheme is really close to the original HLLC scheme [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF], the only difference lies in the intermediate states of the Riemann solver since the two-phase contact is approximated in our solver.

Let us assume a finite-volume mesh with volumes [x

i-1 2 , x i+ 1 2 ] × [t n , t n+1
] of size ∆x × ∆t. We define the cell average of the state vector U:

U n i = 1 ∆x x i+ 1 2 x i-1 2 U(x, t n )dx
The integration of equation [START_REF] Chinnayya | Modelling detonation waves in heterogeneous energetic materials[END_REF] gives us the following relation:

U n+1 i = U n i - ∆t ∆x F n i+ 1 2 -F n i-1 2 - ∆t ∆x H n i+ 1 2 ,--H n i-1 2 ,+ (19)
where we use the following discrete approximations of the fluxes:

F n i+ 1 2 ≈ 1 ∆t t n+1 t n F U(x i+ 1 2 , t) dt and H n i+ 1 2 ,--H n i-1 2 ,+ ≈ 1 ∆t t n+1 t n x i+ 1 2 x i-1 2 H U(x, t) ∂ x (α 1 (x, t)) dxdt
The computation of the fluxes is based on the approximate solution of a Riemann problem at each interface x i+ 1 2 :

F n i+ 1 2 = F U n i , U n i+1 and H n i+ 1 2 ,± = H ± U n i , U n i+1 (20)
Once the intermediate states

(W 1 ) * L , (W 1 ) * M , (W 1 ) * R , (W 2 ) * L , (W 2 ) *
R of the Riemann problem have been determined as explained previously, we use them to calculate the fluxes.

Let us first focus on the computation of the conservative part F (U L , U R ). We use the following notation:

F (U L , U R ) = 0 G k (U L , U R )
The function G 2 (U L , U R ) is defined as follows [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF]:

G 2 (U L , U R ) =              (G 2 ) L if 0 ≤ (S 2 ) L (G 2 ) * L = (G 2 ) L + (S 2 ) L (V 2 ) * L -(V 2 ) L if (S 2 ) L ≤ 0 ≤ (S 2 ) M (G 2 ) * R = (G 2 ) R + (S 2 ) R (V 2 ) * R -(V 2 ) R if (S 2 ) M ≤ 0 ≤ (S 2 ) R (G 2 ) R if (S 2 ) R ≤ 0 (21)
In the case of phase 1, G 1 (U L , U R ) depends on the wave configuration. For configuration (a) when (S 2 ) M ≤ (S 1 ) M , we have:

G 1 (U L , U R ) =                      (G 1 ) L if 0 ≤ (S 1 ) L (G 1 ) * L = (G 1 ) L + (S 1 ) L (V 1 ) * L -(V 1 ) L if (S 1 ) L ≤ 0 ≤ (S 2 ) M (G 1 ) * M = (G 1 ) * R + (S 1 ) M (V 1 ) * M -(V 1 ) * R if (S 2 ) M ≤ 0 ≤ (S 1 ) M (G 1 ) * R = (G 1 ) R + (S 1 ) R (V 1 ) * R -(V 1 ) R if (S 1 ) M ≤ 0 ≤ (S 1 ) R (G 1 ) R if (S 1 ) R ≤ 0 ( 22
)
and for configuration (b) when (S 1 ) M < (S 2 ) M :

G 1 (U L , U R ) =                      (G 1 ) L if 0 ≤ (S 1 ) L (G 1 ) * L = (G 1 ) L + (S 1 ) L (V 1 ) * L -(V 1 ) L if (S 1 ) L ≤ 0 ≤ (S 1 ) M (G 1 ) * M = (G 1 ) * L + (S 1 ) M (V 1 ) * M -(V 1 ) * L if (S 1 ) M ≤ 0 < (S 2 ) M (G 1 ) * R = (G 1 ) R + (S 1 ) R (V 1 ) * R -(V 1 ) R if (S 2 ) M ≤ 0 ≤ (S 1 ) R (G 1 ) R if (S 1 ) R ≤ 0 (23)
In equations ( 21)-( 23), we use the following notations (for k = 1, 2 and K = L, R):

(V k ) K = (α k ) K   (ρ k ) K (ρ k ) K (u k ) K (ρ k ) K (e k ) K   (G k ) K = (α k ) K   (ρ k ) K (u k ) K (ρ k ) K (u k ) K 2 + (p k ) K (ρ k ) K (e k ) K (u k ) K + (p k ) K (u k ) K  
Moreover the intermediate states are given by the HLLC relations:

(V k ) * K = (α k ) K (ρ k ) K (S k ) K -(u k ) K (S k ) K -(S k ) M     1 (S k ) M (e k ) K + (S k ) M -(u k ) K (S k ) M + (p k ) K (ρ k ) K ((Sk) K -(u k ) K )     Finally, (V 1 ) *
M is defined in the following way:

(V 1 ) * M = (α 1 ) M    (ρ 1 ) * M (ρ 1 ) * M (u 1 ) * M (ρ 1 ) * M ε 1 (ρ 1 ) * M , (p 1 ) * M + 1 2 (ρ 1 ) * M (u 1 ) * M 2   
We recall that, in the case of configuration (a), we have:

         (α 1 ) M = (α 1 ) R (ρ 1 ) * M = (ρ 1 ) * L (u 1 ) * M = (u 1 ) * R (p 1 ) * M = (p 1 ) * R
and configuration (b) gives:

         (α 1 ) M = (α 1 ) L (ρ 1 ) * M = (ρ 1 ) * R (u 1 ) * M = (u 1 ) * L (p 1 ) * M = (p 1 ) * L
Then we focus on the computation of the non-conservative part H ± (U L , U R ). The main difficulty lies in the integration of the non-conservative product p I ∂ x (α 1 ) across the two-phase contact as emphasized in [START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF]. Thanks to our approximation of this wave given in equation ( 16), we have the following relation across the wave:

p I ∂ x (α 1 ) = ∂ x (α 1 p 1 ) = -∂ x (α 2 p 2 )
where we first use the preservation of I 4 and then the preservation of I 3 . We underline the fact that the thin layer approximation proposed in [START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF] provides the same relation. Thus we treat those non-conservative products as in [START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF][START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF][START_REF] Daude | On the computation of the Baer-Nunziato model using ALE formulation with HLL-and HLLC-type solvers towards fluid-structure interactions[END_REF]:

H ± (U L , U R ) = 1 ± sign [(S 2 ) M ] 2     (S 2 ) M [(α 1 ) R -(α 1 ) L ] 0 (-1) k+1 (α 2 ) R (p 2 ) * R -(α 2 ) L (p 2 ) * L (-1) k+1 (S 2 ) M (α 2 ) R (p 2 ) * R -(α 2 ) L (p 2 ) * L     (24)
Schwendeman et al. [START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF] also proposed a second-order extension of their Godunov solver. They proposed to use a classical minmod reconstruction of the primitive variables before applying the Riemann solver. Thus, the conservative part of the flux is obtained as well as the second-order contribution of the integral of the non-conservative terms near the two-phase contact. As detailed in [START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF], the second-order contribution of the integral of the non-conservative terms away from the two-phase contact remains and is also computed based on the solution obtained using the Riemann solver. Although it won't be considered in the sequel, a second-order extension of our scheme can be proposed using a similar approach based on our approximate Riemann solver.

Numerical results

This section is dedicated to the assessment of the HLLC-type scheme introduced in the previous section. Due to the Approximation of the two-phase Contact, this solver will be referred as HLLAC in the sequel. We compare it with the Rusanov scheme [START_REF] Crouzet | Approximate solutions of the Baer-Nunziato Model[END_REF] and the original HLLC-type solver from [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF] on 1D as well as 2D test-cases. Those three schemes have been implemented in the fast transient dynamics software Europlexus [START_REF]Commissariat à l'énergie atomique et aux énergies alternatives (CEA)[END_REF][START_REF] Daude | On the computation of the Baer-Nunziato model using ALE formulation with HLL-and HLLC-type solvers towards fluid-structure interactions[END_REF] which has been used for all the following computations.

General equations of states have been used so far but we need to specify them for the numerical test-cases. We consider that both phases are governed by the Stiffened Gas EOS as in most numerical tests from the literature:

ε k (ρ k , p k ) = p k + γ k p ∞ k (γ k -1) ρ k (25)
Thus, two thermodynamic constants γ k and p ∞ k need to be defined for both phases in each test-case.

1D test-cases

Seven 1D Riemann problems have been chosen in order to assess the accuracy, the robustness and the efficiency of the present scheme. Most of them have been proposed in [START_REF] Dallet | A comparative study of numerical schemes for the Baer-Nunziato model[END_REF] where analytical solutions are provided with a maximum accuracy allowing to carry out grid convergence studies. The SG EOS parameters for all test-cases (except test 7) are gathered in table 1. The initial conditions consist of two constant states, detailed in table 2, separated by a discontinuity at x 0 = 0.5 (except for test 4).

The last 1D test-case has been chosen to assess the ability of the proposed solver to deal with complex EOS. In each test, we consider the domain [0, 1] with transmissive boundary conditions at both ends. All computations have been carried out with a Courant number C CF L = max(λ k ) ∆t ∆x = 0.5 . This test proposed in [START_REF] Dallet | A comparative study of numerical schemes for the Baer-Nunziato model[END_REF] consists of an isolated two-phase contact discontinuity. To build such a testcase, a Left state needs to be set first allowing the evaluation of the exact Riemann invariants of the wave, given in equation [START_REF] Coquel | A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model[END_REF]. Then the preservation of those Riemann invariants provides us non-linear equations that need to be solved in order to obtain the Right state (see [START_REF] Crouzet | Approximate solutions of the Baer-Nunziato Model[END_REF] for more details). The exact solution of the corresponding Riemann problem is an isolated two-phase contact discontinuity since the Riemann invariants of the wave are exactly preserved. Thus, this test-case appears to be relevant to assess the consequence of the key approximation of the HLLAC scheme. Figure 2 shows the comparison between the analytical solution and the numerical solutions obtained with the three different schemes on 1000 cells. As expected, the Rusanov scheme is very diffusive (see the results on α 2 for example) compared to the two others which seem to give a similar diffusion at the interface. We can also notice some "ghost waves" with inaccurate intermediate states on the numerical results of the HLLAC scheme linked to the approximation of the two-phase contact (particularly on ρ 1 and p 1 ). Indeed, the Riemann invariants of the approximated two-phase contact (equation ( 16)) are not preserved between the Left and Right states, unlike the exact ones, therefore additional waves appear.

In order to check if those "ghost waves" disappear when refining the grid, a convergence study has been carried out and results are provided in figure 3. For all schemes, we observed the expected rate of convergence 1/2 for first-order schemes due to the occurrence of a contact discontinuity. Moreover, the HLLC scheme appears to be the most accurate solver for all variables as expected, whereas the Rusanov scheme is the least accurate one. We notice that the accuracy of the HLLAC scheme is really close to the accuracy of the HLLC scheme for (α 2 , ρ 2 , p 2 , u 1 ). For the other variables, (ρ 1 , p 1 , u 2 ), the accuracy of the HLLAC scheme is closer to the one of the Rusanov scheme which is not surprising since those three variables correspond to some Riemann invariants of the approximated two-phase contact.

Test 2: Two-phase contact with shock waves

The second test is built on the previous one in [START_REF] Dallet | A comparative study of numerical schemes for the Baer-Nunziato model[END_REF]: one shock wave propagating to the left has been added to the two-phase contact in each phase. The Right state is the same as in the previous test. The new Left state is obtained thanks to jump conditions [START_REF] Coquel | A Positive and Entropy-Satisfying Finite Volume Scheme for the Baer-Nunziato Model[END_REF] with the Left state from the previous test-case. Thus, a classical monophasic (i.e. complying with the jump conditions of the Euler equations) shock wave propagating to the left appears in each phase in addition to the two-phase contact. Figure 4 provides the comparison between the analytical solution and the numerical ones. Convergence studies have also been carried out on this test and are presented in figure 5. The first remark is that the results of HLLAC and HLLC schemes are very close showing that the proposed approximation of the two-phase contact has a minor impact on the overall accuracy of the solver on this test-case. Differences between those two schemes can only be seen on very fine grids. Compared to those schemes, the Rusanov one appears to be very diffusive and less accurate. Except for u 2 , the expected order of convergence 1/2 is retrieved. In the case of u 2 , the observed order of convergence is 1 due to the fact that it is constant across the contact discontinuity and therefore exhibits only a shock wave. An efficiency study is also carried out on this test to assess the computational cost of the proposed scheme. Thus figure 6 shows the error on α 2 from the grid convergence study versus the CPU time of the corresponding computation for the three different schemes. At a given error, we can see that the Rusanov scheme requires a higher computational cost than the other schemes although it is the least expensive scheme on a given mesh. Among both other schemes, the HLLAC one has a lower computational cost than the HLLC one due to the fact that it does not require the use of an iterative method for each local Riemann problem. More precisely, table 3 provides the computational cost of each scheme for a given error (err(α 2 ) = 10 -4 ). The Rusanov scheme appears to be 40 times more expensive than the HLLAC scheme whereas the HLLC one is 2.5 times more expensive than the HLLAC scheme. This latter observation is linked to the fact that 2 iterations are used in the Newton's method involved in the HLLC scheme as suggested in [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF]. The third test-case consists of a stationary isolated two-phase contact proposed in [START_REF] Dallet | A comparative study of numerical schemes for the Baer-Nunziato model[END_REF]. The key point of this test is to assess that the HLLAC scheme keeps the property of the HLLC scheme to exactly preserve this wave. Whereas the two-phase contact is diffused by the Rusanov scheme, figure 7 shows that the HLLAC does preserve it as expected since the neglected terms of the Riemann invariants of this wave are equal to zero in this particular test-case.

Test 4: A general Riemann problem

This test has been used in [START_REF] Dallet | A comparative study of numerical schemes for the Baer-Nunziato model[END_REF] but was initially proposed in [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF]. The exact solution of the Riemann problem is composed by the 6 waves of the model including strong shock waves. Unlike the other tests, the initial discontinuity is located at x 0 = 0.8 . The comparison of the analytical and numerical solutions is shown in figure 8. Once again, the results of the HLLAC and HLLC scheme are very close. This test-case is particularly severe and we can observe that the accuracy of the Rusanov scheme is not sufficient to retrieve the intermediate states on coarse grids.

Test 5: Isolated two-phase contact with vanishing phases

Test 5 has been designed in the same way as test 1 so that it consists of an isolated two-phase contact. However it is also a test with vanishing phases in both states (α k close to 0 or 1) in order to assess the robustness of the schemes. Figure 9 provides the analytical solution and the numerical solutions obtained with the Rusanov and the HLLAC schemes. Indeed the HLLC scheme was not robust enough for this test since it did not preserve the positivity of pressures. We observe that the accuracy of the Rusanov scheme is very poor on this severe test, especially on phase 2. We also study the convergence of the two schemes on this test-case. Surprising results for the HLLAC scheme on coarse grids (for phase 1) are observed on figure 10. Nevertheless we retrieve the order of convergence of 1/2 for both schemes except in the case of the error on u 2 which remains stable, close to round-off errors.

Test 6: A general Riemann problem with vanishing phases

Test 6 is a general Riemann problem with vanishing phases proposed in [START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF]. The exact solution of the problem is composed by five of the six waves of the model, only the GNL wave (S 2 ) R is not present. Figure 11 provides the comparison between the analytical solution and the numerical solutions of the three schemes. We still observe similar results with the HLLAC and HLLC schemes, even on "ghost waves" in phase 2. Moreover the diffusion of the Rusanov scheme, which is striking on α 2 , leads to a poor accuracy when compared to the two other schemes.

Test 7: Complex Equations Of State

This last Riemann problem has been chosen to assess the ability of the proposed solver to deal with any Equation Of State. Thus, a test-case from [START_REF] Saurel | A Multiphase Godunov Method for Compressible Multifluid and Multiphase Flows[END_REF], which studies the interaction between gaseous detonation products (phase 1) and a copper plate (phase 2), is considered. The detonation products are governed by the Jones-Wilkins-Lee (JWL) EOS while the copper is governed by the Cochran-Chan (CC) EOS. Both EOS can be written in the Mie-Gruneisen form:

ε k (ρ k , p k ) = p k -Π k (ρ k ) Γ k ρ k + Q k (ρ k ) ( 26
)
where the functions Π k (ρ k ) and Q k (ρ k ) need to be specify. In the case of the JWL EOS, they depend on the Chapman-Jouguet state (CJ):

Π 1 (ρ) = A exp -R 1 ρ ref 1 ρ + B exp -R 2 ρ ref 1 ρ Q 1 (ρ) = A ρ ref 1 R 1 exp -R 1 ρ ref 1 ρ + B ρ ref 1 R 2 exp -R 2 ρ ref 1 ρ (27) - A ρ ref 1 R 1 exp -R 1 ρ ref 1 ρ CJ - B ρ ref 1 R 2 exp -R 2 ρ ref 1 ρ CJ - p CJ -Π 1 (ρ CJ ) Γ 1 ρ CJ + e CJ
In the case of the CC EOS, we have:

Π 2 (ρ) = A 1 ρ ρ ref 2 E1 -A 2 ρ ρ ref 2 E2 Q 2 (ρ) = A 1 ρ ref 2 (E 1 -1)   ρ ρ ref 2 E1-1 -1   - A 2 ρ ref 2 (E 2 -1)   ρ ρ ref 2 E2-1 -1   -e 0 (28)
Many thermodynamical parameters are required for both EOS and they are provided in table 4. Initially, the detonation products are at the Chapman-Jouguet state and the copper is under atmospheric conditions, as recalled in table 5. Only the HLLAC scheme is considered in this test and different meshes are used. Figure 12 provides the numerical results on the fraction α 2 and the mean variables

ρ = α 1 ρ 1 + α 2 ρ 2 , u = 1 ρ (α 1 ρ 1 u 1 + α 2 ρ 2 u 2 ), p = α 1 p 1 + α 2 p 2 at time t f = 73 × 10 -6 s.
The solution of this test-case consists of a rarefaction wave propagating to the left in the detonation gases, a two-phase contact propagating to the right and a shock wave propagating to the right in the copper. Due to the complexity of the considered EOS, no exact solution is available for this test but the comparison of the numerical solutions on the different meshes shows that the different waves and intermediate states are accurately captured by the HLLAC scheme when refining the grid. This numerical solution is in agreement with the one provided in [START_REF] Saurel | A Multiphase Godunov Method for Compressible Multifluid and Multiphase Flows[END_REF] (Figure 11 page 454). Since we use a first-order scheme, more diffusion is observed in our numerical results. We also recall that Saurel & Abgrall [START_REF] Saurel | A Multiphase Godunov Method for Compressible Multifluid and Multiphase Flows[END_REF] do not consider the Baer-Nunziato closure laws (equation ( 6)) and use instantaneous pressure and velocity relaxations. This test-case demonstrates the ability of the HLLAC to deal with complex EOS, even in the presence of vanishing phases and violent pressure waves.

α 2 ρ 1 u 1 p 1 ρ 2 u 2 p 2 (kg.m -3 ) (m.s -1 ) (Pa) (kg.m -3 ) (m.s -1 ) (Pa)

2D test-cases

Three 2D test-cases have been selected to push forward the assessment of the present HLLAC scheme and compare it with the HLLC and the Rusanov schemes. The efficiency of the three schemes is also briefly examined here.

Shock-bubble interaction

The first 2D test studies the interaction between a shock wave and a Helium bubble. It has been experimentally investigated by Haas and Sturtevant [START_REF] Haas | Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities[END_REF] and widely used to assess models and numerical methods since then [START_REF] Daude | On the computation of the Baer-Nunziato model using ALE formulation with HLL-and HLLC-type solvers towards fluid-structure interactions[END_REF][START_REF] Kokh | An anti-diffusive numerical scheme for the simulation of interfaces between compressible fluids by means of a five-equation model[END_REF][START_REF] Kreeft | A new formulation of Kapila's five-equation model for compressible twofluid flow, and its numerical treatment[END_REF]. A planar shock wave with M = 1.22 is moving in the air towards a cylindrical Helium bubble. More precisely, the computational domain is 267 mm × 89 mm with the initial shock located at 15 mm from the inlet, the radius of bubble is 25 mm and its center is vertically centered at 84.33 mm from the inlet of the domain. The computations are performed on a mesh made of 630 × 210 cells using transmissive inlet and outlet boundary conditions with walls at top and bottom. The two gases, Helium and air, are modeled as ideal gases (p ∞ k = 0 Pa) with γ 1 = 1.648 and γ 2 = 1.4 respectively. The initial conditions of the test are presented in table 6 and the Courant number is set to C CF L = 0.5 .

Figure 13 shows the comparison between numerical results of the three different schemes on the mean density ρ = α 1 ρ 1 + α 2 ρ 2 at different times. As in many of the 1D Riemann problems, the results of the HLLAC and the HLLC schemes are similar. In contrast, the Rusanov scheme is very diffusive as it can be seen on the stationary interface of the bubble before being reached by the shock wave. This diffusion prevents the Rusanov scheme from capturing important physical phenomena in the interaction unlike the other schemes. Thus, we observe with the two HLLC-type schemes that the refracted shock wave in the bubble is ahead the shock due to the higher speed of sound in Helium. Moreover, a transmitted shock wave appears at the location where the reflected rarefaction wave meets the interface and it forms a λ-shock by interacting with the incoming shock. We also retrieve the beginning of the Richtmyer-Meshkov instability due to the interaction between the shock wave and the curved interface of the bubble. Finally we compare the CPU time of the computations of 1 ms with the three different schemes. Although, they provide very different results the Rusanov and the HLLAC schemes have a similar CPU cost whereas the HLLC scheme is 2.8 times slower on this test.

Underwater explosion

The second 2D test is taken from [START_REF] Grove | Anomalous reflection of a shock wave at a fluid interface[END_REF][START_REF] Kokh | An anti-diffusive numerical scheme for the simulation of interfaces between compressible fluids by means of a five-equation model[END_REF] and focuses on an underwater explosion. A rectangular domain [-2; 2] m × [-1.5; 1] m with wall boundary conditions is considered. Initially, a circular bubble of high-pressure gas is located in the water under the air-water interface located at y = 0. The radius of the bubble is 0.12 m and its center is located at (0; -0.3) m. The air above the interface and the gas in the bubble are modeled as an ideal gas (p ∞ 1 = 0 Pa) with γ 1 = 1.4 . The Stiffened Gas EOS is used for the water (p ∞ 2 = 6 × 10 8 Pa and γ 2 = 4.4). Table 7 provides the initial conditions of the computations. This test is particularly suited to assess the robustness of the schemes due to the presence of a strong shock wave and material interfaces between vanishing phases. The initial pressure and density ratios are particularly high: around 10 4 for pressure and 10 3 for density. We use a grid with 400 × 250 cells and a Courant number C CF L = 0.5 . Moreover the pressure and velocity Rusanov scheme HLLC scheme HLLAC scheme Rusanov scheme HLLC scheme HLLAC scheme Figure 13: Shock and air/He bubble interaction: comparison between numerical solutions obtained with the Rusanov, the HLLC and the HLLAC schemes at different times.

relaxation terms [START_REF] Hérard | A fractional step method to compute a class of compressible gasliquid flows[END_REF][START_REF] Crouzet | Validation of a two-fluid model on unsteady liquid-vapor water flows[END_REF][START_REF] Lochon | Comparison of two-fluid models on steam-water transients[END_REF] are also used in this test-case. Details on those source terms can be found in appendix B. The time scales associated to those phenomena are set to very small values in order to obtain quasi-instantaneous relaxations: τ p = 10 -12 s and τ u = 10 -12 s. Figure 14 provides the numerical results of the computations with the different schemes on the mean density ρ = α 1 ρ 1 + α 2 ρ 2 at several times. First, the HLLC scheme was not robust enough for this severe test. As in test 5, the HLLC does not preserve the positivity of the pressures of the intermediate states probably due to the presence of vanishing phases. Once again, the results obtained with the Rusanov and the HLLAC schemes look very different due to the inner numerical diffusion of the Rusanov scheme. Indeed, the initial pressure difference between the bubble and the water generates a shock wave in the water and a rarefaction wave in the bubble. Then, the shock wave reaches the free surface resulting in a reflected rarefaction wave in the water and a transmitted shock wave in the air which is much weaker and cannot be observed on the present density results. In the case of the Rusanov scheme, the horizontal air-water interface is very smeared. As it can be observed on the results obtained with the HLLAC scheme, the shape of the growing bubble evolves from a circular one to an oval one. More details on the physical phenomena can be found in [START_REF] Liu | The simulation of compressible multi-medium flow: II. Applications to 2D underwater shock refraction[END_REF]. The results obtained with the Rusanov scheme show a very different shape of the bubble due to its poor accuracy. Finally we note that the CPU cost of the HLLAC scheme is only 15% higher than the one of the Rusanov scheme on this test demonstrating the efficiency of the proposed scheme.

Cylindrical underwater explosion near a planar wall

The third 2D test is similar to the previous one except that a rigid planar wall replaces the free surface. Therefore the strong interaction between a shock wave in the water and a planar wall can be investigated. This test-case was previously studied in [START_REF] Xie | Application of a one-fluid model for large scale homogeneous unsteady cavitation: The modified Schmidt model[END_REF][START_REF] Ma | A GPU based compressible multiphase hydrocode for modelling violent hydrodynamic impact problems[END_REF]. A rectangular computational domain [-6; 6] m ×[-6; 3] m is considered. The upper boundary condition located at y = 3 m is a rigid wall whereas the other boundary conditions are transmissive. Initially, a 1-m radius explosive gas bubble is located at (0; 0). This gas is modeled as an ideal gas (p ∞ 1 = 0 Pa) with γ 1 = 2 whereas the water is modeled as a Stiffened Gas with p ∞ 2 = 3 × 10 8 Pa and γ 2 = 7.15. Table 8 provides the initial conditions used for the present computations. As in [START_REF] Ma | A GPU based compressible multiphase hydrocode for modelling violent hydrodynamic impact problems[END_REF], a small amount of gas with α 1 = 0.005 is uniformly distributed in the water. A very fine mesh is used with 1200 × 900 cells. All of the computations are performed with a Courant number C CF L = 0.5. As in the previous underwater explosion test, the pressure and velocity relaxation terms [START_REF] Hérard | A fractional step method to compute a class of compressible gasliquid flows[END_REF][START_REF] Crouzet | Validation of a two-fluid model on unsteady liquid-vapor water flows[END_REF][START_REF] Lochon | Comparison of two-fluid models on steam-water transients[END_REF] are used with τ p = 10 -12 s and τ u = 10 -12 s (see appendix B for details on those source terms). density (ρ = α 1 ρ 1 + α 2 ρ 2 ) fields in conjunction with the mean pressure (p = α 1 p 1 + α 2 p 2 ) contours at different time instances. First, as in the previous underwater explosion test, the HLLC scheme is not robust enough on this challenging test. At t = 1.5 ms, the main shock generated by the explosion has been reflected from the upper planar rigid wall (see Figure 15-(a)). The excessive numerical diffusion of the Rusanov scheme smears both the density jump located at the gas-water interface and the shock waves which propagate in the water. Then, the reflected shock wave interacts with the gas-water interface resulting in a rarefaction wave propagating in the water towards the upper wall and a weak shock wave propagating in the gas bubble (with a lower celerity) which are clearly visible at t = 2 ms (see Figure 15-(b)). Notice that the rarefaction wave is strongly affected by the numerical dissipation of the Rusanov scheme. As a consequence, the obtained reflection of this rarefaction wave at the upper wall is weaker with the use of the Rusanov scheme than with the use of the HLLAC scheme. Due to this reflection, a low pressure region is generated near the wall leading to an inception of cavitation.

Notice that the cavitation is only obtained with the HLLAC scheme. Afterward, the cavitation pocket collapses as it can be seen on the numerical results obtained with the HLLAC scheme. Figure 15-(e) shows that, at t = 5.5 ms, the cavitation region near the center of the wall has collapsed whereas cavitation regions still exist on both sides of the collapse zone.

The pressure and the gas volume fraction histories at the center of the wall are presented on Figure 16. All of the computations retrieve the main pressure peak at t = 1 ms which corresponds to the reflection of the main shock wave at the upper wall. The time of occurrence and the shape of this peak are in agreement with the previous computations performed by Xie et al. [START_REF] Xie | Application of a one-fluid model for large scale homogeneous unsteady cavitation: The modified Schmidt model[END_REF] (see Figure 16 page 1190) and by Ma et al. [START_REF] Ma | A GPU based compressible multiphase hydrocode for modelling violent hydrodynamic impact problems[END_REF] (see Figure 13 page 16). The maximal amplitude of the pressure peak (6600 bar with HLLAC and 6400 bar with Rusanov) is very close to the Xie's solution (6500 bar). The generation of cavitation at t = 3.7 ms is only retrieved with the the HLLAC scheme. In contrast, no cavitation occurs during the computation using the Rusanov scheme as it can also be noticed on the time evolution of the gas volume fraction. However, the collapse of the cavitation pocket is predicted earlier in the present computation (t = 4.55 ms) in comparison with the Xie's solution (t = 5.05 ms) and the Ma's solution (t = 5.25 ms). This discrepancy may be explained by the use of first-order schemes in the present computations in contrast to the previous ones [START_REF] Xie | Application of a one-fluid model for large scale homogeneous unsteady cavitation: The modified Schmidt model[END_REF][START_REF] Ma | A GPU based compressible multiphase hydrocode for modelling violent hydrodynamic impact problems[END_REF]. In the same way, the earliest stages of the gas volume fraction history obtained with the HLLAC scheme are in agreement with the Ma's solution (see Figure 14 page 16), i.e. strong drop at t = 1 ms when the shock wave reflects at the wall and rapid increase at t = 3.7 ms when cavitation occurs. However, the gas volume fraction strongly reduces at t = 4.55 ms in the present computation which corresponds to the collapse of the cavitation pocket. Besides the maximal value of the gas volume fraction is also lower in the present computation in comparison to the Ma's computation (0.0165 vs 0.04-0.05). In order to investigate the influence of the numerical errors induced by the use of first-order schemes, a finer mesh 2400 × 1800 cells is also considered to reduce the corresponding numerical diffusion errors. The associated results show that the collapse of the cavitation pocket is postponed (t = 4.7 ms) and the maximum value of α 1 is increased (0.04).

Conclusion

Based on the solver proposed by Tokareva and Toro [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF], a new HLLC-type Riemann scheme for the computation of the Baer-Nunziato model has been proposed. The key idea of that new scheme is an approximation of the two-phase contact of the model based on the Riemann invariants of this wave. It is only valid for "subsonic" flows but physically-relevant in those cases. An important consequence of the approximation is that the solver can directly deal with any Equation Of State. Moreover the solver no longer requires the resolution of a non-linear equation using an iterative procedure since the intermediate states are obtained as the solution of a linear system unlike the original one. Finally, the approximation of the two-phase contact is in agreement with the thin layer approximation proposed by Schwendeman et al. [START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF] and therefore a similar numerical treatment of the non-conservative products is used. The proposed solver has been compared with the original one and the Rusanov scheme on both 1D and 2D test-cases. Several 1D Riemann problems have been selected in order to assess the accuracy and the robustness of the scheme. The last 1D Riemann problem also demonstrates its ability to deal with complex Equations Of State. Grid convergence studies have been carried out and particular attention has been paid to the case of vanishing phases. Numerical results show that the original HLLC scheme is the most accurate one but the new HLLC-type scheme provides similar results in most cases and appears to be more robust on challenging tests. An efficiency study has also been carried out showing that the proposed scheme is more efficient than the original HLLC one. Three 2D test-cases have then been studied, a shock-bubble interaction and two underwater explosions which are severe test-cases with violent shock waves and material interfaces. Once again, the new scheme exhibits accuracy and robustness. Moreover, it appears to be very efficient since its computational cost is similar to the Rusanov one on those tests whereas it is more accurate. Thus, the proposed scheme appears to be an excellent candidate for the computation of more realistic situations since it is a good compromise between the accuracy of the original HLLC scheme and the robustness of the Rusanov scheme, while showing great efficiency.

Only first-order schemes have been considered in the present work but the extension of the proposed solver to the second-order, using the approach of Schwendeman et al. [START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF] for example, should be assessed. Based on previous work [START_REF] Daude | On the computation of the Baer-Nunziato model using ALE formulation with HLL-and HLLC-type solvers towards fluid-structure interactions[END_REF], the solver should also be tested in the ALE framework, particularly its robustness, in order to tackle the computations of two-phase flows with fluid-structures interactions.

Appendix A. Connection through the two-phase contact

This appendix is dedicated to the connection through the two-phase contact between two states W l and W r where we use the following notation:

W K =           (α 1 ) K (ρ 1 ) K (u 1 ) K (p 1 ) K (ρ 2 ) K (u 2 ) K (p 2 ) K           , K = l, r
More precisely, we assume W l , (α 1 ) r , (ρ 2 ) r and look for W r such that:

I k (W r ) = I k (W l ) , k = 1, 5
where (I k ) k=1,5 are the exact Riemann invariants of the two-phase contact given in equation [START_REF] Coquel | A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model[END_REF]. To that end, we first define the main unknown x = (ρ 1 ) r and recast the preservation of I 4 in the following way:

I 4 (W l ) = h 1 x, (p 1 ) r + 1 2 I 2 (W l ) (α 1 ) r x 2 (29)
where h 1 (ρ 1 , p 1 ) is the specific internal enthalpy of phase 1. We also define the function P 1 such that p 1 = P 1 (ρ 1 , s 1 ). Thus, equation ( 29) can be written as a non-linear equation on x: f (x) = I 4 (W l ) with f (x) = h 1 x, P 1 x, I 5 (W l ) + 1 2

I 2 (W l ) (α 1 ) r x 2 (30)
Once this equation has been solved, W r can easily be fully determined. Indeed, (p 1 ) r = P 1 x, I 5 (W l ) , (u 1 ) r and (u 2 ) r are obtained thanks to the preservation of I 1 and I 2 . Finally the last Riemann invariant I 3 provides us the remaining unknown (p 2 ) r . Thus the connection through the two-phase contact only depends on the resolution of the non-linear equation [START_REF] Ishii | Thermo-fluid dynamic theory of two-phase flows[END_REF].

For classical EOS, the first term of f (x) is an increasing function of x. For example, the ideal gas EOS gives:

h 1 x, P 1 x, I 5 (W l ) = γ 1 I 5 (W l ) γ 1 -1 x γ1-1
Therefore, for non-zero values of I 2 (W l ), the function f (x) admits a minimum for x = x 0 . Moreover f (x) is decreasing for x in ]0, x 0 ] and then increasing for x in [x 0 , +∞[. As a consequence, the equation ( 30) has no solution if f (x 0 ) > I 4 (W l ), one solution if f (x 0 ) = I 4 (W l ) and two solutions if f (x 0 ) < I 4 (W l ) which corresponds to figure 17. It must be underlined that, when f (x 0 ) > I 4 (W l ), no state W r can be connected to W l by the two-phase contact.

Figure 17: Resolution of equation [START_REF] Ishii | Thermo-fluid dynamic theory of two-phase flows[END_REF]

Figure 1 :

 1 Figure 1: Structure of the exact solution of the Riemann problem for the Baer-Nunziato model in the "subsonic" situation. On the left, configuration (a) where (S 2 ) M ≤ (S 1 ) M . On the right, configuration (b) where (S 1 ) M < (S 2 ) M .

Figure 2 :

 2 Figure 2: Comparison between analytical and numerical solutions on Test 1 with 1000 cells at time t f = 0.1s . From the top to the bottom right: α 2 , ρ 1 , ρ 2 , p 1 , p 2 , u 1 , u 2 vs x.

Figure 3 :

 3 Figure 3: Convergence study on Test 1 at time t f = 0.1s . From the top to the bottom right: L1 norm of the error on α 2 , ρ 1 , ρ 2 , p 1 , p 2 , u 1 , u 2 vs ∆x.

Figure 4 :

 4 Figure 4: Comparison between analytical and numerical solutions on Test 2 with 100 cells at time t f = 0.1s . From the top to the bottom right: α 2 , ρ 1 , ρ 2 , p 1 , p 2 , u 1 , u 2 vs x.

Figure 5 :

 5 Figure 5: Convergence study on Test 2 at time t f = 0.1s . From the top to the bottom right: L1 norm of the error on α 2 , ρ 1 , ρ 2 , p 1 , p 2 , u 1 , u 2 vs ∆x.

Figure 6 : 24 4. 1 . 3

 62413 Figure 6: Efficiency study on Test 2 at time t f = 0.1s : L1 norm of the error on α 2 vs CPU time of the computation.

Figure 7 :Figure 8 :

 78 Figure 7: Comparison between analytical and numerical solutions on Test 3 with 100 cells at time t f = 0.05s . From the top to the bottom right: α 2 , ρ 1 , ρ 2 , p 1 , p 2 , u 1 , u 2 vs x.

Figure 9 :

 9 Figure 9: Comparison between analytical and numerical solutions on Test 5 with 1000 cells at time t f = 0.01s . From the top to the bottom right: α 2 , ρ 1 , ρ 2 , p 1 , p 2 u 1 , u 2 vs x.

Figure 10 :Figure 11 :

 1011 Figure 10: Convergence study on Test 5 at time t f = 0.01s . From the top to the bottom right: L1 norm of the error on α 2 , ρ 1 , ρ 2 , p 1 , p 2 , u 1 , u 2 vs ∆x.

Figure 12 :

 12 Figure 12: Numerical solutions of Test 7 obtained with the HLLAC scheme on different meshes at time t f = 73 × 10 -6 s. From the top left to the bottom right: α 2 , ρ, u, p vs x.

Figure 15 Figure 14 :

 1514 Figure 15 provides the numerical results of the computation with the different schemes on the mean

Figure 15 :

 15 Figure 15: Cylindrical underwater explosion near a planar wall: comparison between numerical solutions obtained with the Rusanov (left) and HLLAC (right) schemes at different instants: (a) t = 1.5 ms, (b) t = 2 ms, (c) t = 3 ms, (d) t = 4 ms, and (e) t = 5.5 ms.

Figure 16 :

 16 Figure 16: Cylindrical underwater explosion near a planar wall: comparison of the time evolutions of the pressure (a) and the gas volume fraction (b) obtained with the HLLAC and Rusanov schemes.

Table 1 :

 1 SG EOS parameters for the 1D Riemann problems

		γ 1	γ 2	p ∞ 1	p ∞ 2
	Test 1	3	2	0.	0.
	Test 2	3	2	0.	0.
	Test 3	1.4	3	0.	10.
	Test 4	1.4	3	0.	100.
	Test 5	3	2	0.	0.
	Test 6	1.4	3	0.	0.

Table 2 :

 2 Initial conditions for the 1D Riemann problems

			α 2	ρ 1	u 1	p 1	ρ 2	u 2	p 2
	Test 1 [14]	Left state Right state	0.8 0.2	1.25 1.2539534	0.1837722 0.2793675	1. 1.0095183	1.5185185 0.3110988 0.5 0.3110988 1.7771419 1.2
	Test 2 [14]	Left state Right state	0.8 0.2	1. 1.2539534	0.5 0.2793675	0.5 1.0095183	1. 0.5	0.8 0.3110988 1.7771419 0.5
	Test 3 [14]	Left state Right state	0.6 0.3	1.4 1.	0. 0.	1. 1.	1.4 1.	0. 0.	2. 3.
	Test 4 [42]	Left state Right state	0.7 0.2	1. 1.	-19.5975 -19.5975	1000. 0.01	1. 1.	-19.5975 -19.5975	1000. 0.01
	Test 5	Left state Right state	1 -10 -6 10 -6	1.25 1.25421479 0.311098673 1.01014963 0.1837722 1.	1.5185185 0.3110988 0.5 0.3110988 189851.199 1.2
	Test 6 [41]	Left state Right state	0.5 10 -6	1. 1.8	0. 0.	2. 4.	2. 1.7829	0. 0.2972	5. 3.5422

Table 4 :

 4 JWL and CC EOS parameters for Test 7

	JWL EOS	CC EOS
	Γ 1	0.25	Γ 2	2
	ρ ref 1 A (Pa) (kg.m -3 )	1840 8.545 × 10 11	ρ ref 2 A 1 (Pa) (kg.m -3 )	8900 1.45667 × 10 11
	B (Pa)	2.05 × 10 10	A 2 (Pa)	1.47751 × 10 11
	R 1	4.6	E 1	2.994
	R 2	1.35	E 2	1.994
	e CJ (J.kg -1 )	2.610793 × 10 6 e 0 (J.kg -1 )	1.17900 × 10 5
	p CJ (Pa)	37 × 10 9
	ρ CJ (kg.m -3 )	2485.37

Table 5 :

 5 Initial conditions for Test 7

Table 6 :

 6 Initial conditions for the shock-bubble interaction test-case

		α 1	ρ 1	u 1	p 1	ρ 2	u 2	p 2
			(kg.m -3 ) (m.s -1 ) (Pa) (kg.m -3 ) (m.s -1 )	(Pa)
	Right state	10 -6	0.2546	-114.42 10 5	1.92691	-114.42 1.5698 × 10 5
	Left state	10 -6	0.2546	0.	10 5	1.4	0.	10 5
	Bubble	1 -10 -6	0.2546	0.	10 5	1.4	0.	10 5

Table 7 :

 7 Initial conditions for the underwater explosion test-case

		α 1	ρ 1	u 1	p 1	ρ 2	u 2	p 2
			(kg.m -3 ) (m.s -1 )	(Pa)	(kg.m -3 ) (m.s -1 )	(Pa)
	Water	10 -6	1.225	0.	1.01325 × 10 5	1000	0.	1.01325 × 10 5
	Air	1 -10 -6	1.225	0.	1.01325 × 10 5	1000	0.	1.01325 × 10 5
	Gas bubble 1 -10 -6	1250	0.	10 9	1000	0.	10 9

Table 8 :

 8 Initial conditions for the computation of the 2D cylindrical underwater explosion near a

	planar wall
		α 1	ρ 1	u 1	p 1	ρ 2	u 2	p 2
			(kg.m -3 ) (m.s -1 )	(Pa)	(kg.m -3 ) (m.s -1 )	(Pa)
	Water	0.005	1.27	0.	10 5	1000	0.	10 5
	Gas bubble 0.995	1270	0.	8.29 × 10 8	1000	0.	8.29 × 10 8

Acknowledgements The first author received a financial support through the EDF-CIFRE contract 0459-2013. Computational facilities were provided by EDF. Numerical simulations have been performed with the Europlexus software.

Appendix B. Pressure and velocity relaxations

The Baer-Nunziato model with the source terms including pressure and velocity relaxations can be written as follows:

where k = 1, 2 and j = 3 -k. The scalars K p and K u are positive functions given by [START_REF] Hérard | A fractional step method to compute a class of compressible gasliquid flows[END_REF]:

where τ p and τ u are respectively the time scales of pressure and velocity relaxation. System (31) is approximated using a fractional step method proposed in [START_REF] Hérard | A fractional step method to compute a class of compressible gasliquid flows[END_REF]. More details on the different relaxation substeps can be found in [START_REF] Hérard | A fractional step method to compute a class of compressible gasliquid flows[END_REF][START_REF] Crouzet | Validation of a two-fluid model on unsteady liquid-vapor water flows[END_REF][START_REF] Lochon | Comparison of two-fluid models on steam-water transients[END_REF].