Nicolas Grenier

Virginie Daru

Patrick Le Quéré

Keywords: sharp interface, compressible two-phase flow, low Mach correction, phase change vi interface raide, écoulement diphasique compressible, correction à faible Mach, changement de phase sharp interface, compressible two-phase flow, low Mach correction, phase change

The work of this thesis would not be possible without the help and support around me, both profes-

List of Symbols

Two-phase flows are present in a wide range of industrial processes and could impact production efficiency and security. In the chemical industry, bubble columns are used to increase the exchange surfaces, either to promote heat exchange or chemical reactions. In petroleum engineering, gas is injected into an oil field to increase the rate of oil extraction. In thermal power plants, the steam production (vaporization in liquid phase) intensifies heat transfer by latent heat. In the storage of natural gas, to optimize the amount of gas stored into a limited volume, natural gas can be cooled down to its condensing temperature. At this temperature, the liquid-vapor phase change transforms the gas into liquid form and thus increases its density. Moreover, thermal control of electronic components of satellites is ensured by liquid-vapor cooling systems such as heat pipes or loop heat pipes. Heat release is becoming a critical issue with the increase of the power of electronic components. For all these configurations, the density ratio between the two-phases could be very important. Phase change intensifies the heat transfer and in the liquid, there exists phase change and surface tension effects at the interface. Moreover, the material velocity could be very small compared to the sound velocity (low Mach flows). It is therefore essential to understand and anticipate the behavior of these two-phase flows for the design, optimization and control of industrial equipment. Moreover, in accidental scenarios, a good knowledge of two-phase flows is important in order to find an adequate solution in order to deal with the problem encountered and thus guarantee safety. A numerical model to predict with confidence the physical behavior of low-Mach flows involving a liquid-vapor interface that occurs around saturation temperature (at constant thermodynamic pressure) is necessary. However, many difficulties and problems exist in the mathematical modeling of two-phase flows, the mathematical analysis of established models and numerical methods, and also in the numerical resolution of these mathematical models.

Scientific context

In this section, we will first recall the methods used for interface description during the simulation of two-phase flows. A list of the different existing solvers to simulate compressible two-phase flows will follow. Finally, a short presentation of the strategy adopted for the development of our numerical model will be made.

The Interface

We hereafter recall the definition of an interface. For any multiphase study, zones of small thicknesses in which the characteristic gradients of certain physical quantities such as viscosity and density are much more important than in the rest of the field are always encountered. This zone of small thicknesses is called the interface.

Since the original and pioneering work of Harlow and Welch [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF], a substantial number of interface description methods have been developed. There are two main classes of interface description methods: sharp/diffuse approaches and Lagrangian/Eulerian approaches.

Sharp and diffuse interface methods

x x

Field variable Field variable The diffuse method has a mixing zone around the interface as shown in Fig. 1.1, this makes the presence of mixture of different material at the interface. This makes it possible to obtain a continuous function for fluid variables. Different models exist to treat this mixing zone. Abgrall [START_REF] Abgrall | How to prevent pressure oscillations in multicomponent flow calculations: a quasi conservative approach[END_REF] and Karni [START_REF] Karni | Multicomponent flow calculations by a consistent primitive algorithm[END_REF] propose a four-equation model while Allaire et al. [START_REF] Allaire | A five-equation model for the simulation of interfaces between compressible fluids[END_REF] and Murrone and Guillard [START_REF] Murrone | A five equation reduced model for compressible two phase flow problems[END_REF] develop models with five equations. The five-equation model proposed by Allaire et al. [START_REF] Allaire | A five-equation model for the simulation of interfaces between compressible fluids[END_REF] combines the momentum conservation equation, the energy conservation equation, two mass conservation equation (each phase) as well as the transport of volume fraction equation. To close the system, isobaric or isothermal conditions are used. Within the diffused numerical interface non-physical mixing states may occur, for which an artificial equation of state has to be defined [START_REF] Fechter | A sharp interface method for compressible liquid-vapor flow with phase transition and surface tension[END_REF]. In order to control the thickness of the mixing zone, Shyue and Xiao [START_REF] Shyue | An eulerian interface sharpening algorithm for compressible two-phase flow: the algebraic thinc approach[END_REF] use a hyperbolic tangent profile to represent the interface while Kokh and Lagoutiere [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] utilize an anti-diffusive numerical scheme. The sharp interface method treats the interface as a contact discontinuity, there is no mixture at the interface. In the present study, we concentrate on this sharp interface method.

To avoid the presence of mixture at the interface, we assume that the thickness of this zone will be negligible and that the gradients will be integrated and replaced by a jump condition (sharp interface). This thin zone is reduced by one dimension. Thus, mathematically, the interface Γ is an object of dimension d -1 which evolves in a space D ⊂ R d of dimension d and separates the global region into two-separated regions as presented in Fig. 1.2. Some physical quantities will undergo a discontinuity which equals the jump condition across the interface. In general, we use this notion of interface to evoke the separation surface that exists between two phases. in Fig. 1.3. Their advantage resides in their capacity to take into account very precisely the conditions of jump through the interface (these conditions of jump are in fact boundary conditions between the two domains separated by the mesh-interface). On the other hand, these methods lack flexibility when one wishes to follow complex interface movements since it is necessary to reconstruct the mesh at each time step. Background fixed mesh based method Another type of Lagrangian approach has been developed with more success on fixed meshes, for example, the Front-tracking method [102]. The interface is represented by a set of points connected to each other moving on a background fixed grid as presented in Fig. Although this approach has some advantages such as explicit representation and Lagrangian transport of the interface, this method requires frequent markers rearrangements because of interface motion that is an obvious shortcoming. This procedure could be complex especially when the interface suffers from significant deformations. Keeping some homogeneity of distribution of markers on the interface as well as tackling interface break-up and coalescence could also be a big challenge. The reconstruction of the interface and the calculation of its geometric characteristics (curvature, normal vector) is not immediate and often requires additional complex and expensive algorithms. All these drawbacks are greatly amplified in three-dimensional simulations.

There exist other Lagrangian methods where a set of points are not connected. In Marker And Cell (MAC) method [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF], a cloud of points presents the interface while the equations of motion are solved on background fixed mesh. In meshless methods, such as Smoothed Particle Hydrodynamics (SPH) [START_REF] Hu | A multi-phase sph method for macroscopic and mesoscopic flows[END_REF], Moving Particle Semi-implicit (MPS) [START_REF] Koshizuka | Numerical analysis of breaking waves using the moving particle semi-implicit method[END_REF], Particle Finite Element Method (PFEM) [START_REF] Idelsohn | The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves[END_REF],

• • • , the interface and Navier-Stokes equation are solved on Lagrangian points, known as particles.

Eulerian approach

The Eulerian approach considers the interface as a scalar field, which will be transported by the velocity field thanks to a transport equation. This approach can avoid mesh or marker rearrangement issues and overcome difficulties of the Lagrangian approach. There are mainly two different algorithms for the Eulerian approach: the Volume Of Fluid method [START_REF] Hirt | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF] and the Level Set method [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations[END_REF].

Volume Of Fluid method The Volume Of Fluid (VOF) method, specially designed for simulating two-phase flows, is based on the principle of mass conservation. This method is an Eulerian approach with a scalar field representing the volume fraction of a fluid inside a control volume as shown in Fig. 1.5. The interface is represented implicitly, defined at the instant t > 0 by the points x ∈ R d where a volume function equals to 0.5. The interface reconstruction is necessary as depicted in Fig. 1.5. This [START_REF] Lafaurie | Modelling merging and fragmentation in multiphase flows with surfer[END_REF][START_REF] Puckett | A high-order projection method for tracking fluid interfaces in variable density incompressible flows[END_REF][START_REF] Gueyffier | Volume-of-fluid interface tracking with smoothed surface stress methods for three-dimensional flows[END_REF][START_REF] Sato | A sharp-interface phase change model for a mass-conservative interface tracking method[END_REF].

Despite this good property of conservation, VOF suffers from a serious problem: insufficient information to represent the interface. As presented in Fig. 1.5, valuable information to capture the interface only exists in a narrow band of width about 2 or 3 times the grid size, and therefore the location of the interface is diffuse, making it difficult to accurately calculate the geometry quantities of the interface (normal vector, curvature). Numerical diffusion in the transport scheme could cause non-physical smearing of the interface. Thus, parasitic current phenomena occur when problems of high-density jump are studied (for example water-air interface), making the numerical scheme inaccurate. More importantly, mesh refinement can not significantly decrease the impact of these phenomena. In addition, for three-dimensional problems, this algorithm becomes complex and time-consuming, because it is necessary to perform a reconstruction of the interface from the volume fraction field.

Level Set method The Level Set method is based on a scalar field of distance to the interface. Being the points x ∈ R d , the value of the Level Set function [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations[END_REF] φ (x x x,t) is, at each instant of time t > 0, the minimum distance of the point from the interface Γ, with a sign depending on which side it is contained as shown in Fig. 1.6. This function (signed distance function [START_REF] Fedkiw | Level set methods and dynamic implicit surfaces[END_REF]) φ (x x x,t) can be seen as an infinite set of contour lines each of them at a fixed distance from the interface. The interface Γ(t) that is defined at the instant t > 0 by the points x ∈ R d where a function (x,t) → φ (x,t) vanishes, i.e. Γ(t) = {x ∈ D | φ (x,t) = 0}. This is one possible option to describe a Level Set function. It is however the simplest and has the advantage to allow the point x to know how far is the interface. The domain where the fluid i = 1, 2 lies will be noted

D i (t) = {x ∈ D | (-1) i φ (x,t) < 0}. Any quantity Z becomes t > 0, x ∈ D i (t) Z (x x x,t) = Z i (x x x,t)
The so formulated Level Set has some distinct advantages. The first is that the interface is exactly known at each time step, no reconstruction is needed. The interface geometry can be directly obtained from the values of φ , the normal vector

x → y, n n n = ∇φ |∇φ | , for y ∈ Γ(t), (1.1)
and the interface curvature

x → y, κ = ∇ • ∇φ |∇φ | , for y ∈ Γ(t). (1.2)
The second advantage is its hyper regularity property |∇φ

(x x x)| = 1 ∀ x ∈ D.
Being a signed distance function, no discontinuity appears on the interface, with the obvious advantages of enhanced accuracy, stability and diffusion in the numerical schemes.

Compared with Lagrangian approaches and the VOF method, the Level Set method avoids the problem of mesh or maker rearrangement and the reconstruction of interface. Numerical diffusion in the scheme of interface transport is significantly reduced, providing an accurate estimation of the interface location and geometries. Considering the good properties presented by the Level Set method, we selected this method to implicitly represent the interface.

Simulation of low Mach two-phase flows

Flows of interest contain one liquid phase and one gas phase separated by an interface. As the Mach number (ratio between the material velocity and the sound velocity of the medium) in a liquid is generally small, there are generally two different approaches for liquid modeling: pressure-based approach [START_REF] Fedkiw | A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF][START_REF] Lalanne | On the computation of viscous terms for incompressible two-phase flows with Level Set/Ghost Fluid method[END_REF][START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF][START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF][START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF] and density-based approach [START_REF] Fechter | A sharp interface method for compressible liquid-vapor flow with phase transition and surface tension[END_REF][START_REF] Houim | A ghost fluid method for compressible reacting flows with phase change[END_REF][START_REF] Hu | On the HLLC riemann solver for interface interaction in compressible multi-fluid flow[END_REF][START_REF] Lin | Simulation of compressible two-phase flows with topology change of fluid-fluid interface by a robust cut-cell method[END_REF].

Pressure-based approaches are the most natural approach to simulate low Mach or incompressible flows. However, these methods require solving Poisson's equation [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF]. In other words, a linear system needs to be resolved even without impliciting the system which will cost too many computational resources. Compared with the pressure-based methods, the density-based methods do not require solving Poisson's equation which is a better choice for high speed or high Mach flows. Classical density-based methods used in the low Mach regime face a known issue: when the flow is slow compared to sound speed, common Godunov type compressible solvers lose accuracy due to excessive numerical dissipation [START_REF] Dellacherie | Analysis of godunov type schemes applied to the compressible euler system at low mach number[END_REF]. Common Godunov type Riemann solver can not capture the features of low Mach flows, we call this feature as low Mach capturing. In this work, we decide to implement the density-based methods with low Mach capturing.

With the low Mach capturing feature and interface description method (sharp or diffuse), the density-based approaches in the literature could be mainly divided into three categories: diffuse interface with low Mach capturing [START_REF] Chalons | An all-regime Lagrange-Projection like scheme for 2D homogeneous models for two-phase flows on unstructured meshes[END_REF][START_REF] Fuster | An all-Mach method for the simulation of bubble dynamics problems in the presence of surface tension[END_REF][START_REF] Pelanti | Low Mach number preconditioning techniques for Roe-type and HLLCtype methods for a two-phase compressible flow model[END_REF][START_REF] Peluchon | A robust implicit-explicit acoustic-transport splitting scheme for two-phase flows[END_REF], diffuse interface without low Mach capturing [START_REF] Jain | A conservative diffuse-interface method for compressible two-phase flows[END_REF][START_REF] Ansari | Numerical simulation of compressible two-phase flow using a diffuse interface method[END_REF][START_REF] So | Anti-diffusion interface sharpening technique for two-phase compressible flow simulations[END_REF][START_REF] Shyue | An eulerian interface sharpening algorithm for compressible two-phase flow: the algebraic thinc approach[END_REF] and sharp interface without low Mach capturing [START_REF] Deng | Simulating compressible two-medium flows with sharp-interface adaptive runge-kutta discontinuous galerkin methods[END_REF][START_REF] Fechter | A sharp interface method for compressible liquid-vapor flow with phase transition and surface tension[END_REF][START_REF] Lin | Simulation of compressible two-phase flows with topology change of fluid-fluid interface by a robust cut-cell method[END_REF][START_REF] Shen | A 3d conservative sharp interface method for simulation of compressible two-phase flows[END_REF]. Unfortunately, there are very few studies concern a sharp interface with low Mach capturing issues. However, the sharp interface approach with low Mach capturing can better describe the compressible effects, better capture low Mach feature, more importantly, there is no artificial mixture across the interface.

According to all these analyses in the literature, in the present work, based on the density-based "Lagrange-projection" type scheme [START_REF] Chalons | An all-regime lagrange-projection like scheme for the gas dynamics equations on unstructured meshes[END_REF] for single fluid flows, we will develop a sharp interface approach with low Mach capturing for two-phase low Mach flows with phase change and surface tension effects at the interface. To achieve such a goal, each phase will be treated as a pure phase without mixture and described by the density-based "Lagrange-projection" type solver with low Mach capturing feature. The interface will be captured by the Level Set method and the Level Set function should be updated with a high-order numerical scheme to get a precise description of the interface. Since the interface is considered sharp, based on the Level Set method, each cell is completely treated as a cell of fluid 1 or of fluid 2 to avoid the small cut-cell problem. The phase interface is shifted to the grid cell boundary of the nearest neighboring grid cell. This is the most efficient way to reconstruct the interface. To take into account phase change and surface tension effects at the interface, with pure states beside the interface, an approximate Riemann solver account for jumps conditions (related to phase change and surface tension) is developed. Across the interface, the density ratio could be very important, which will impose difficulties on the interface coupling. Special treatments and analysis are necessary to recover the good property of convergence.

Outline of the thesis

Chapter 2 is dedicated to the physical model of the fluid solver which contains the description of the equations and associated hypothesis, in particular the physical effects which are considered in the computation, such as the gravity, viscosity, heat transfer, surface tension and phase change.

Chapter 3 describes the principle of operator splitting and high-order methods. With the acoustictransport splitting, the global model given in Chapter 2 is decomposed into three subsystems: acoustic, transport and diffusion.

In Chapter 4, different schemes for interface description issues are presented. For the Level Set advection problem, we first propose the coupled time-space approach called One-Step scheme [START_REF] Daru | High order one-step monotonicity-preserving schemes for unsteady compressible flow calculations[END_REF] and compare its performance with classical schemes. As for the Level Set redistancing problem, we evaluate the performance of two different approaches: PDE based and Hopf-Lax based methods. As for the interface curvature approximation, we assess the performance of the classical Level Set method and the Height function method. By evaluating the performance of the numerical scheme on related issues with well-known test cases, we determinate the adopted schemes for the interface description.

Chapter 5 gives the approximation of the global system except the Level Set advection equation. As for the phase change model, to take the interface temperature as a boundary condition, in the literature [START_REF] Sato | A sharp-interface phase change model for a mass-conservative interface tracking method[END_REF][START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF], an implicit method is always applied for the stability issue. To avoid this implicit procedure, in this chapter we propose an explicit method to approximate the velocity jump. In the present work, we are interested in low Mach flows. For a low Mach regime simulated by a compressible solver, excessive numerical diffusion is involved, the conventional compressible solvers lose accuracy. To get a precise prediction, the low Mach correction is necessary. As the interface is sharp, without a smooth region at the interface, such discontinuity existing at the interface presents a huge challenge for the design of a numerical scheme. The well-known low Mach fix in literature could lead to significant truncation error, especially for two-phase flows with large density and sound speed ratios. To recover a good asymptotic-preserving property, in Chapter 6, we propose a new low Mach correction with rigorous asymptotic analysis.

Chapter 7 is devoted to the validation of the present numerical model.

Chapter 2

Global Model

This chapter is devoted to present the physical model that has been used to describe two-phase flows including effects of viscosity, surface tension, external forces, heat transfer and phase change. We give general governing equations to predict the fluid behavior in each phase. The jump condition at the interface is finally presented to give the coupling conditions between different phases.

Governing equations

We suppose that both fluids are separated by a sharp interface Γ(t) = {x ∈ D | φ (x,t) = 0} at the instant t > 0 and consider pure phase (without mixture) on each side of the interface. The evolution of each phase could be described by a compressible model. We shall suppose that each fluid i = 1, 2 verifies the compressible Navier-Stokes equations in D i (t) that is to say:

t > 0, x ∈ D i (t)        ∂ t ρ + ∇ • (ρu) = 0, ∂ t (ρu) + ∇ • (ρu ⊗ u) + ∇p = ∇ • S + ρ f f f , ∂ t (ρE) + ∇ • ((ρE + p)u) = ∇ • (Su) + ρ f f f • u u u + ∇ • (K ∇T ), (2.1a) (2.1b) (2.1c)
where ρ is the fluid density, p is the the fluid pressure, K is the thermal conductivity, T is the

fluid temperature, S = µ[∇u + (∇u) T ] - 2 3
µ∇ • uI d is the viscous stress tensor following the Stokes hypothesis, µ is the viscosity, and I is the identity tensor of size d. The vectorial function f f f contains the volume forces, like the gravity g g g. In absence other volume forces, f f f equals to g g g. E denotes the total energy per unit of mass which is related to the specific internal energy e via E = e + 1 2 |u| 2 .

This system is closed by adding an equation of state (EOS). We consider two compressible fluids i = 1, 2 are both equipped with a Mie-Gruneisen Equation of State (EOS) of the form [START_REF] Saurel | A multiphase godunov method for compressible multifluid and multiphase flows[END_REF][START_REF] Goncalvès | Constraints on equation of state for cavitating flows with thermodynamic effects[END_REF]:

(ρ, e) → p EOS i (ρ, e) = (γ i -1)ρe -π ∞ i γ i , (2.2)
and

ρ(γ i -1)C v,i T = p + π ∞ i . (2.
3)

The sound speed c i is then defined for the fluid i by

c 2 i = (∂ p EOS i /∂ ρ) + (p EOS i /ρ 2 )(∂ p EOS i /∂ e) = γ i p i + π ∞ i ρ i .
γ being the heat capacity ratio (γ i = C p,i /C v,i ), where C p and C v are the heat capacities at constant pressure and volume respectively, π ∞ is a constant pressure, representing the molecular attraction between fluid molecules. We want to underline that other EOS such as barotropic EOS (see details in Appendix. A) could be also employed to close the system. In multiphase flows, the interface is transported by the fluid velocity field u. The Level Set φ function is then passively advected by the local velocity field, it is therefore the solution of the transport which follows the linear hyperbolic advection equation [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations[END_REF]:

∂ t φ + u Γ • ∇φ = 0, ∀x ∈ D,t > 0 (2.4)
where u Γ is the velocity of the interface.

Interface jump conditions

The jump conditions at the interface are obtained from conservation principles and are provided by the terms appearing in the surface integrals. As presented in the precedent section, the interface Γ is considered with zero-thickness and treated as a contact discontinuity.

Jump condition on mass

The mass conservation can be expressed as:

ρ (u -u Γ ) • n Γ = 0, (2.5)
Z Γ = lim x→y x∈D 1 Z -lim x→y x∈D 2
Z for y ∈ Γ(t), t > 0 is the jump operator across the interface in the normal direction n. With the phenomenon of phase change, by introducing the interface mass flux ṁ, the mass conservation can be rewritten as:

u • n Γ = ṁ 1 ρ Γ . (2.6)
Without phase change, the conservation equation can be expressed as:

u • n Γ = 0, (2.7)
the continuity of the normal velocity at the interface is ensured.

Jump condition on stress

According to the equation of momentum, the jump condition on pressure [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF] can be expressed as:

p Γ = σ κ + 2 µ ∂ u n ∂ n Γ - 1 ρ Γ ṁ2 , (2.8)
where ∂ u n ∂ n is the normal derivative of the normal velocity component. For flows at equilibrium, (2.8) verifies the Laplace's law:

p Γ = σ κ.

Jump condition on energy

While phase change occurs at the interface, a jump could be observed on the thermal gradient as presented in [START_REF] Juric | Computations of boiling flows[END_REF]. By assuming that the temperature at the interface equals to the saturation temperature (T Γ = T sat ), we can get:

K ∇T Γ • n = ṁL heat , (2.9)
where L heat is the latent heat of vaporization, T Γ and T sat are the temperature at the interface and the saturation temperature respectively.

Global model

In this chapter, we determined the global model which consists of the combination of compressible Navier-Stokes equations (2.1) and the Level Set advection equation (2.4) for x ∈ D i (t), t > 0:

             ∂ t ρ + ∇ • (ρu) = 0, ∂ t (ρu) + ∇ • (ρu ⊗ u) + ∇p = ∇ • S + ρ f f f , ∂ t (ρE) + ∇ • ((ρE + p)u) = ∇ • (Su) + ρ f f f • u u u + ∇ • (K ∇T ), ∂ t φ + u Γ • ∇φ = 0, (2.10a) (2.10b) (2.10c) (2.10d)
the EOS to close the system (2.10) as well as jump conditions across the interface Γ(t) to couple sub-domains D 1 (t) and D 2 (t). The following chapters will be devoted to numerically resolve this system with a satisfactory prediction.

Chapter 3

Operator splitting

To resolve any complicated problem, an efficient way could be decomposing it into simpler parts with specific characteristics. According to related characteristics, each part could be resolved separately with the most adapted methods. By combining separate solutions, this complicated problem is solved with a satisfactory result. This technique is called operator splitting which could optimize and simplify the global resolution. Based on different motivations, such as directional splitting, split linear from nonlinear, and split terms related to different physical phenomena, splitting technique can be extended to an arbitrary number of steps [START_REF] Macnamara | Operator splitting[END_REF]. We present in the following sections the basic idea of the splitting technique. This basic idea is then extended by Strang [START_REF] Strang | On the construction and comparison of difference schemes[END_REF] to recover high-order accuracy. To decompose the global system, we first introduce two splitting techniques for the Euler system: Lagrangian-Projection method and acoustic-transport splitting. We split the global system into three subsystems: acoustic with source terms, transport and diffusion.

Basic idea

The idea of operator splitting could be traced to the 1870s and introduced by Lie and Engel [START_REF] Lie | Theorie der transformationsgruppen[END_REF]. For any linear ordinary differential equation (ODE) [START_REF] Glowinski | Splitting methods in communication, imaging, science, and engineering[END_REF][START_REF] Holden | Splitting methods for partial differential equations with rough solutions: Analysis and MATLAB programs[END_REF]:

∀ t > 0, x ∈ D i ,      dU dt + (A 1 + A 2 )U = 0, U(x,t = 0) = U 0 (x). (3.1a) (3.1b)
where U 0 is the initial function, A 1 and A 2 are two matrices. The solution of this time-depend problem can be expressed as:

U(x,t) = e -(A 1 +A 2 )t U 0 (x). (3.2)
With the relation lim

k→+∞ e -A 1 t k e -A 2 t k k = e -(A 1 +A 2 )t
, by defining t n = n∆t , n ∈ N + , the approximate update from t n to t n+1 with Lie's splitting can be given as:

U(x,t n+1 ) = e -A 1 ∆t e -A 2 ∆t U(x,t n ), (3.3)
or

U(x,t n+1-) = e -A 1 ∆t U(x,t n ), U(x,t n+1 ) = e -A 2 ∆t U(x,t n+1-). (3.4a)
(3.4b)

High-order splitting

The Taylor series on the solution of the ODE problem (3.2) is that:

lim e ∆t(A 1 +A 2 ) U(x,t n ) = U(x,t n ) I + ∆t(A 1 + A 2 ) + 1 2! ∆t 2 (A 1 + A 2 ) 2 + O(∆t 3 ) , (3.5)
(3.5) shows that (3.3) is the splitting of first-order accuracy. To recover higher-order accuracy, Strang [START_REF] Strang | On the construction and comparison of difference schemes[END_REF] introduced a splitting called "Strang splitting" . From t n to t n+1 , the second-order Strang splitting can be given as:

n ∈ N + , x ∈ D i U(x,t n+1 ) = U(x,t n )(e -1 2 A 2 ∆t e -A 1 ∆t e -1 2 A 2 ∆t ), (3.6)
or

n ∈ N + , x ∈ D i U(x,t n+1 ) = U(x,t n )(e -1 2 A 1 ∆t e -A 2 ∆t e -1 2 A 1 ∆t ). (3.7)
(3.6) and (3.7) provide second-order accuracy, while this version of Strang splitting involves the half time step 1 2 ∆t. To avoid the half time step, a sequential symmetrically weighted splitting is developed which is a combination of two Lie splittings in different ordering [START_REF] Csomós | Weighted sequential splittings and their analysis[END_REF]. For each two time-steps from t n to t n+2 , its algorithm reads:

n ∈ N + , x ∈ D i U(x,t n+2 ) = U(x,t n ) (e -A 1 ∆t e -A 2 ∆t ) t n → t n+1 (e -A 2 ∆t e -A 1 ∆t ) t n+1 → t n+2 . ( 3.8)
Every two time steps, (3.8) recovers second-order accuracy. Higher-order accuracy splitting methods are theoretically available and have been introduced in many work [START_REF] Blanes | Optimized high-order splitting methods for some classes of parabolic equations[END_REF][START_REF] Jia | A third accurate operator splitting method[END_REF][START_REF] Yoshida | Construction of higher order symplectic integrators[END_REF]. However, they are difficult to implement, moreover, third or higher-order schemes involve sub-steps going backward in time which could be a source of solution instability, usually due to intermediate boundary conditions. Second-order splitting is still the most implemented method.

Euler system splitting

Before involving the full global system, we concentrate first on the Euler system in each pure phase:

t > 0, x ∈ D i (t)        ∂ t ρ + ∇ • (ρu) = 0, ∂ t (ρu) + ∇ • (ρu ⊗ u) + ∇p = 0, ∂ t (ρE) + ∇ • (u(ρE + p)) = 0. (3.9a) (3.9b) (3.9c)
As presented in Section 1.2.1, Lagrangian approaches have their own advantages such as explicit tracking of the interface, while the Eulerian approaches could avoid the problem of mesh degeneration or grid reconnection. To treat the Euler system in Eulerian grid and give consideration to the advantages of Lagrangian approaches (such as without mixture and explicit treatment of jump condition), two different splitting methods are presented to solve the Euler system in Eulerian grid.

Lagrange-Projection method

The Lagrange-Projection method could be divided into two steps: a Lagrangian evolution step followed by a projection step. We start from the Eulerian coordinate on which the physical quantities are known at t = 0. We carry out the Lagrangian step which describes the evolution of physical quantities in a Lagrangian formalism. The second step is then carried out which is the projection. This step consists in projecting the physical quantities after the Lagrangian update onto the initial Eulerian grid.

Lagrangian step

In the Lagrangian step, the evolution of the physical system considered is described by a set of equations which translates the evolution of physical quantities by adding evolution of the geometry of the system. Denoting by u = u(x,t) the fluid velocity of position x at t > 0, we consider the differential system:

for t > 0 dx dt = u(x,t). (3.10)
By defining ξ ξ ξ ∈ R d as the Lagrangian coordinate, we denote by t → x(ξ ξ ξ ,t) and the initial condition of (3.10) reads:

x(0) = ξ ξ ξ . (3.11)
(ξ ξ ξ ,t) are the Lagrangian coordinates associated with the velocity u(x,t) and (3.11) means that the the initial Lagrangian coordinate ξ ξ ξ coincides with Eulerian coordinate x. To follow the material, the Lagrangian grid undergoes a transformation:

J(ξ ξ ξ ,t) = ∂ x ∂ ξ ξ ξ (ξ ξ ξ ,t). (3.12)
We obtain, at the end of this step, a deformed coordinate ξ * . (3.9) can be set into the Lagrangian formulation [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF]:

t > 0, ξ ξ ξ ∈ R d        ∂ t (ρJ) = 0, ∂ t (ρuJ) + ∇ ξ ξ ξ p = 0, ∂ t (ρEJ) + ∇ ξ ξ ξ • (pu) = 0. (3.13a) (3.13b) (3.13c)
with ∇ ξ ξ ξ = J∇.

Projection step

By defining Ω E j ⊂ D E and Ω L k ⊂ D L as arbitrary volumes in Eulerian and Lagrangian space, the projection step for any quantity Z could be written as:

Z E j = Z L dΩ E j = 1 |Ω E j | ∑ ∀k, Ω E j ∩Ω L k ̸ =∅ Z L k |Ω E j ∩ Ω L k |, (3.14)
where Z L and Z E are the quantity Z after the Lagrangian and projection update. For example, as depicted in Fig. 3.1, at time t n = n∆t, the Lagrangian coordinates ξ ξ ξ n et Eulerian coordinates x coincide. In the Lagrangian zone, the Eulerian coordinates are transport by the material velocity and deformed to x n+1 * . Then with a projection step, we project Lagrangian quantities onto the Eulerian grid. The projection to a cell j could be expressed:

Z E j = Z L j+1 |Ω E j ∩ Ω L j+1 | + Z L j |Ω E j ∩ Ω L j | |Ω E j | . (3.15) j + 1 j x j-1/2 x j+1/2 x j+3/2 ξ ξ ξ n j-1/2 ξ ξ ξ n j+1/2 ξ ξ ξ n j+3/2
x n+1 * j+3/2

x n+1 * j+1/2

x n+1 * j-1/2 ∆t Fig. 3.1 Grids evolution of Lagrangian-Projection method. green area: Ω E j ∩ Ω L j , blue area:

Ω E j ∩Ω L j+1 .
We can observe that: the mixture of physical quantities in different volumes Ω is only associated with the projection step. This splitting method leaves the mixture into the projection procedure. Compared to direct Eulerian methods (One-Step Eulerian, such as the Roe scheme [START_REF] Roe | Approximate riemann solvers, parameter vectors, and difference schemes[END_REF]), the use of the Lagrangian phase allows more precise treatment of interface problems and facilitates interface coupling.

Acoustic-transport splitting

To have similar functionality as the Lagrange-Projection and avoid the possible complex geometry problem in the projection step, here we introduce the acoustic-transport splitting as given in [START_REF] Chalons | An all-regime lagrange-projection like scheme for the gas dynamics equations on unstructured meshes[END_REF]. The goal is to identify an acoustic and a transport operator to do the splitting.

We change variables of the Euler system (3.9) from conservative variables to primitive variables (ρ, u, p). Then one can show that the Euler system is equivalent to the following system, easier to manipulate:

t > 0, x ∈ D i (t)          ∂ t ρ + (u u u • ∇)ρ + ρ∇ • u u u = 0, ∂ t u u u + (u u u • ∇)u u u + 1 ρ ∇p = 0, ∂ t p + (u u u • ∇)p + ρc 2 ∇ • u u u = 0, (3.16a) (3.16b) (3.16c)
where c is the speed of sound in the fluid. We here consider a two dimensional problem u u u = (u, v), and an extension to three dimensions is straightforward. A more compact form can be written:

∂ t W + A (W, n ′ ) ∇W = 0 where W =       ρ u v p       , A (W ) =       u u u • n ′ ρ|n ′ 1 | ρ|n ′ 2 | 0 0 u u u • n ′ 0 |n ′ 1 |/ρ 0 0 u u u • n ′ |n ′ 2 |/ρ 0 ρc 2 |n ′ 1 | ρc 2 |n ′ 2 | u • n ′      
, where n ′ is a unit vector and n ′ 1 and n ′ 2 are its two components with

|n ′ 1 | 2 + |n ′ 2 | 2 = 1.
The matrix A is diagonalizable in R and the eigenvalues are: [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF]. The acoustic splitting of the Euler system consists in noticing that we can decompose A into two operators, A = A acoustic +A transport with A acoustic = diag (-c, 0, c) and A transport = (u u u • n ′ )I 3 . As A = A acoustic + A transport , the Euler system can be rewritten in the following non-conservative form:

u u u • n ′ -c, u u u • n ′ , u u u • n ′ + c
∂ t W + A acoustic ∇W + A transport ∇W = 0.
For the Euler equations, we can get two subsystems by the splitting technique, the acoustic subsystem is given by:

t > 0, x ∈ D i (t)          ∂ t ρ + ρ∇ • u u u = 0, ∂ t u + 1 ρ ∇p = 0, ∂ t p + ρc 2 ∇ • u u u = 0. (3.17a) (3.17b) (3.17c)
By the same way, the transport subsystem is given by:

t > 0, x ∈ D i (t)        ∂ t ρ + (u u u • ∇)ρ = 0, ∂ t u u u + (u u u • ∇)u u u = 0 0 0, ∂ t p + (u u u • ∇)p = 0. (3.18a) (3.18b) (3.18c)
For the acoustic subsystem, we make a change of variable (ρ, u u u, p) to (τ, u u u, E), with τ = 1 ρ . We obtain the hyperbolic system:

t > 0, x ∈ D i (t)        ∂ t τ -τ∇ • u u u = 0, ∂ t u u u + τ∇p = 0 0 0, ∂ t E + τ∇ • (pu u u) = 0. (3.19a) (3.19b) (3.19c)
Similar to the Lagrange-Projection, (3.19) involves only the acoustic phenomenon, and the mixture effects are splitted into the transport step. Moreover, the eigenvalues of the acoustic and transport operators are associated with the sound speed c and fluid material velocity u respectively. The acoustic-transport splitting provides an efficient way to construct a semi-implicit scheme for flows in the low Mach regime where the material velocity is much slower than the sound speed.

Decomposition of the global system

According to different phenomena presented in the global system (2.10), one possible approach is to split the acoustic, advection, external force, diffusion effects during a short time-interval [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF], so the global system can be divided into the acoustic subsystem:

t > 0, x ∈ D i (t)              ∂ t ρ + ρ∇ • u u u = 0, ∂ t (ρu) + ρu∇ • u + ∇p = 0 0 0, ∂ t (ρE) + ρE∇ • u + ∇ • (pu u u) = 0, ∂ t φ = 0, (3.20a)
(3.20b) (3.20c) (3.20d)
the transport subsystem:

t > 0, x ∈ D i (t)              ∂ t ρ + u u u • ∇ρ = 0, ∂ t (ρu) + (u • ∇)ρu = 0 0 0, ∂ t (ρE) + u • ∇(ρE) = 0, ∂ t φ + u Γ • ∇φ = 0, (3.21a) (3.21b) (3.21c) (3.21d)
the diffusion subsystem:

t > 0, x ∈ D i (t)              ∂ t ρ = 0, ∂ t (ρu) = ∇ • S, ∂ t (ρE) = ∇ • (Su) + ∇ • (K ∇T ), ∂ t φ = 0, (3.22a) (3.22b) (3.22c) (3.22d)
and the external force subsystem:

t > 0, x ∈ D i (t)              ∂ t ρ = 0, ∂ t (ρu) = ρ f f f , ∂ t (ρE) = ρ f f f • u u u, ∂ t φ = 0. (3.23a) (3.23b) (3.23c) (3.23d)
When we want to study a dynamic phenomenon near equilibrium, we need special treatment of the source term. Indeed, the problem of the proposed splitting is that a continuous equilibrium is not the solution of the scheme. Then the background state that should be at equilibrium is not preserved. A way to overcome this issue is to use a well-balanced scheme. A well-balanced scheme aims at preserving discrete versions of some continuous equilibrium states. It means that the scheme has its own version of equilibrium and it is able to preserve it around machine precision. We follow the work of Chalons et al. [START_REF] Chalons | A large time-step and wellbalanced lagrange-projection type scheme for the shallow water equations[END_REF] which derived a well-balanced scheme for the Saint-Venant equations and add well-chosen terms in the scheme to exactly compensate pressure gradients at equilibrium.

We saw that the scheme we are using is obtained by a splitting involving the acoustic, transport, diffusion and the external force subsystems. To obtain an equilibrium, we need the pressure to balance the external source term. It is then natural to put the external source term in the acoustic subsystem:

t > 0, x ∈ D i (t)              ∂ t ρ + ρ∇ • u u u = 0, ∂ t (ρu) + ρu∇ • u + ∇p = -ρ∇Ψ, ∂ t (ρE) + ρE∇ • u + ∇ • (pu u u) = -ρu • ∇Ψ, ∂ t φ = 0, (3.24a) (3.24b) (3.24c) (3.24d)
where Ψ is the gravitational potential, f f f = -∇Ψ. For hydrostatic balance, the flow is at rest and the external force is balanced by the pressure. This hydrostatic balance can be expressed as:

∇p = -ρ∇Ψ, u = 0. (3.25)
With the techniques of splitting, the compressible Navier-Stokes equations are decoupled to three subsystems: the acoustic subsystem (3.24), the transport subsystem (3.21) and the viscous subsystem (3.22). In the next subsection, we will discretize these subsystems and gives a numerical resolution for the fluid system.

Given a fluid state (ρ, ρu u u, ρE, φ ) n j in a cell j , the splitting algorithm decomposes the global system (2.10) as follows:

1. Update the fluid state (ρ, ρu u u, ρE, φ ) n j to the value (ρ, ρu u u, ρE, φ ) n+ j by approximating the solution of the acoustic subsystem (3.24).

2. Update the fluid state (ρ, ρu u u, ρE, φ ) n+ j to the value (ρ, ρu u u, ρE, φ ) n+1- j by approximating the solution of the transport subsystem (3.21).

3. Update the fluid state (ρ, ρu u u, ρE, φ ) n+1- j to the value (ρ, ρu u u, ρE, φ ) n+1 j by approximating the solution of the viscous subsystem (3.22).

We want to mention that, by omitting the external force and the Level Set advection, the acoustic and transport subsystems recover a conservative Euler system. To guarantee a conservative update for the conservative variables (ρ, ρu, ρE), the acoustic subsystem is discretized with fluid states at t n , while the transport subsystem is discretized with the fluid state at t n+ , that is the first-order Lie's splitting. High-order Strang type splitting is not implemented to the acoustic-transport splitting. The reason for this choice will be justified later. To avoid long-chain splitting and large stencil, the diffusion subsystem is discretized with fluid states at t n .

Chapter 4

Interface description

In the precedent chapter, we defined the global model for two-phase flows and split it into three subsystems. As the interface description is an essential issue for two-phase flows, a precise numerical resolution is indispensable. To have a precise description of the evolution of the interface, we here apply high-order schemes for Level Set advection. This section follows by implementing different numerical methods to describe the temporal evolution of the Level Set function and comparing their performances in terms of both accuracy and efficiency. During Level Set advection, different contour lines are advected variously. The distance property of the Level Set function is not guaranteed which leads to the inaccuracy problem to curvature estimation, a procedure called "redistancing" is necessary to preserve the distance property. In this chapter, we will focus on numerical schemes for the Level Set advection and redistancing problem.

Numerical resolution of Level Set advection

We consider a two-dimensional advection problem, and extension to 3D will be straightforward. The linear hyperbolic advection equation's discrete form in two dimensions looks like:

for t > 0, ∂ t φ + uφ x + vφ y = 0. (4.1)
u and v are two components of the velocity in the direction of x and y respectively. The numerical methods for the advection equation can be divided into two families: coupled time-space approach and separate time-space approach. In this work, the retained schemes for the first family will be OS (One-Step) type schemes with different constraints [START_REF] Daru | High order one-step monotonicity-preserving schemes for unsteady compressible flow calculations[END_REF]. The separate WENO (Weighted Essentially Non-Oscillatory) [START_REF] Jiang | Efficient implementation of weighted eno schemes[END_REF] or HOUC (High-order Upstream Central Scheme) [START_REF] Nourgaliev | High-fidelity interface tracking in compressible flows: Unlimited anchored adaptive level set[END_REF] type spatial discretization and Runge-Kutta type time discretization will be employed for the second family. Numerical resolution of this equation will be given by both the separate time-space approaches and the coupled time-space approach.

Separate time-space approach

The separate time-space approach is usually based on multi-stage time integration. The most recent popular time-space methods use a Runge-Kutta type time discretization. In each stage of time integration, a high-order spatial discretization is applied. Considering the computational accuracy and efficiency, usually we apply a scheme equipped with a 5 th or 7 th order spatial discretization and a 3 rd -order Runge-Kutta time discretization.

WENO type scheme for spatial discretization

Considering the semi-discrete equation and leaving the temporal derivative latter, we now concentrate on the discretization of spatial derivatives. With upwind schemes, for a grid point (x i , y j ) at t n , in the x direction the first derivative of φ can be written as:

(φ x ) i, j =    (φ x ) + i, j , if u i, j < 0, (φ x ) - i, j , if u i, j > 0.
(4.2) (φ x ) + i, j and (φ x ) - i, j signify the forward and backward choices of the first derivative in the x-direction facing positive and negative velocity u. The 5 th -order of accuracy WENO scheme (WENO5) [START_REF] Jiang | Efficient implementation of weighted eno schemes[END_REF] offers an approximation of the spatial derivative (φ x ) + i, j , (φ x ) - i, j by combining three upwind decentered derivatives based on three different stencils of 3 rd -order ENO (Essentially Non-Oscillatory) scheme:

(φ x ) ± i, j = 3 ∑ k=1 ω k (φ x ) ±,k i, j , (4.3)
for the values of k = 1, 2, 3, the derivatives are defined as

         (φ x ) ±,1 i, j = + 1 3 q ± 1 -7 6 q ± 2 + 11 6 q ± 3 , (φ x ) ±,2 i, j = -1 6 q ± 2 + 5 6 q ± 3 + 1 3 q ± 4 , (φ x ) ±,3 i, j = + 1 3 q ± 3 + 5 6 q ± 4 -1 6 q ± 5 , with    q + k = φ i-3+k, j -φ i-4+k, j △x , q - k = φ i+4-k, j -φ i+3-k, j △x . (4.4)
Based on the expression of α ± , the coefficients ω k give a convex combination of three different stencils of ENO scheme (ω 1 + ω 2 + ω 3 = 1) that gives a (2r -1) th order of accuracy (r is the order of the ENO scheme)

ω ± k = α ± k / 3 ∑ i=1 α ± i , (4.5)
where α can be expressed as:

           α ± 1 = 1 10 1 ε+IS ± 1 2 , α ± 2 = 6 10 1 ε+IS ± 2 2 , α ± 3 = 3 10 1 ε+IS ± 3 2 , (4.6)
where ε ensures the denominator to be different from zero and the IS are indicators of regularity can be given as:

         IS ± 1 = 13 12 (q ± 1 -2q ± 2 + q ± 3 ) 2 + 1 4 (q ± 1 -4q ± 2 + 3q ± 3 ) 2 , IS ± 2 = 13 12 (q ± 2 -2q ± 3 + q ± 4 ) 2 + 1 4 (q ± 2 -q ± 4 ) 2 , IS ± 3 = 13 12 (q ± 3 -2q ± 4 + q ± 5 ) 2 + 1 4 (3q ± 3 -4q ± 4 + q ± 5 ) 2 .
(4.7)

HOUC type scheme for spatial discretization

The high-order upstream central scheme for Level Set advection was proposed by Nourgaliev and Theofanous [START_REF] Nourgaliev | High-fidelity interface tracking in compressible flows: Unlimited anchored adaptive level set[END_REF]. The spatial discretization is based on Taylor series expansion and can be constructed up to 11 th -order of accuracy. This scheme is efficient and easy to implement. The spatial discretization of HOUC is upwind, thus (φ x ) i, j can be expressed by (φ x ) + i, j and (φ x ) - i, j

for u positive and u negative respectively as in (4.2). 5 th -order of accuracy spatial discretization can be given as follows:

(φ x ) ± i, j = 1 60 [∓2φ i∓3, j ± 15φ i∓2, j ∓ 60φ i∓1, j ± 20φ i, j ± 30φ i±1, j ∓ 3φ i±2, j ]. (4.8)
7 th -order of accuracy spatial discretization can be expressed by:

(φ x ) ± i, j = 1 420 [ ± 3φ i∓4, j ∓ 28φ i∓3, j ± 126φ i∓2, j ∓ 420φ i∓1, j ± 105φ i, j ± 252φ i±1, j ∓ 42φ i±2, j ± 4φ i±3, j ]. (4.9)
The discretization of φ x in x-direction for separate time-space approaches has been given in the precedent section. With the same principle, we can derive spatial discretization in other directions.

Temporal discretization

To complete the design of separate time-space approaches, we need to discretize the temporal derivative. In terms of accuracy and efficiency, here, we employ a 3 rd -order Runge-Kutta scheme from Gottlieb and Shu [START_REF] Gottlieb | Total variation diminishing runge-kutta schemes[END_REF]. If the spatial discretization of the advection term in (2.4) can be written as

∂ t φ = J (φ ). (4.10)
According to the work of Gottlieb and Shu [START_REF] Gottlieb | Total variation diminishing runge-kutta schemes[END_REF] and knowing the solution at time nδt, the solution with 3 rd -order temporal discretization at (n + 1)δt can be estimated as:

         φ 1 = φ n + ∆tJ (φ n ), φ 2 = 3 4 φ n + 1 4 φ 1 + ∆tJ (φ 1 ) , φ n+1 = 1 3 φ n + 2 3 φ 2 + ∆tJ (φ 2 ) .
(4.11)

Coupled time and space approach

To present the coupled time-space approach, here we apply a Strang directional splitting to the equation and consider the one-dimensional problem:

∂ t φ + u∂ x φ = 0, (4.12)
the Cauchy-Kowalewski procedure [START_REF] Del Pino | Arbitrary high-order schemes for the linear advection and wave equations: application to hydrodynamics and aeroacoustics[END_REF] applied to (4.12) that reads:

∂ k t φ = (-u) k ∂ k x φ , 0 ≤ k ≤ n. (4.13)
Supposing u is positive, for a point at (x i ,t n ), Taylor series analysis of φ lead to:

(∂ x φ ) i = φ n i -φ n i-1 ∆x - N ∑ k=2 (-∆x) k-1 k! ∂ k x φ + O(∆x N ), ∂ t φ i = φ n+1 i -φ n i ∆t - N ∑ k=2 ∆t k-1 k! ∂ k t φ + O(∆t N ) = φ n+1 i -φ n i ∆t + u N ∑ k=2 (-ν i ∆x) k-1 k! ∂ k x φ + O(∆t N ), (4.14)
where ν i is the CFL number ν i = u i ∆t/∆x and the discretization is upwind. Then by combining the two sub-equations in (4.14), we can get

∂ t φ i + u i (∂ x φ ) i = φ n+1 i -φ n i ∆t + u n i φ n i -φ n i-1 ∆x -u n i N ∑ k=2 (1 -ν k-1 i )(-∆x k-1 ) k! ∂ k x φ + O(∆x N ) + O(∆t N ). (4.15)
We can get the explicit second-order scheme for a positive u i :

φ n+1 i = φ n i -ν i (φ n i -φ n i-1 ) + ν i -ν 2 i 2 (φ n i+1 -2φ n i + φ n i-1 ) + O(∆x 2 ) + O(∆t 2 ). (4.16)
(4.16) can be given by an explicit Lax-Wendroff form:

φ n+1 i = φ n i -ν i (F lw i+1/2 -F lw i-1/2 ). (4.17)
The second-order "flux" of Lax-Wendroff type can be given by:

F lw i+ = 1 2 (φ n i + φ n i+1 ) - 1 -(1 -ν i ) 2 (φ n i+1 -φ n i ),
F lw i-= 1 2 (φ n i + φ n i-1 ) - 1 -(1 -ν i ) 2 (φ n i -φ n i-1 ). (4.18)
Here F i+ and F i-could be estimated as flux at the right and left of the the grid point x i . As the Level Set advection equation is not in a conservative form, the flux is not conservative at "face" F i+ ̸ = F i+1-.

The modified scheme reads:

∂ t φ + u∂ x φ = u ∆x 2 6 (ν 2 -1)∂ 3 x φ + O(∆x 3 ) + O(∆t 3 ). (4.19)
The explicit third-order upwind scheme with a numerical flux can be given by:

F 3 i+ = 1 2 (φ n i + φ n i+1 ) - 1 -(1 -ν i ) 2 (φ n i+1 -φ n i ) + ν 2 i -1 6 (φ n i+1 -2φ n i + φ n i-1 )
,

F 3 i-= 1 2 (φ n i + φ n i-1 ) - 1 -(1 -ν i ) 2 (φ n i -φ n i-1 ) + ν 2 i -1 6 (φ n i -2φ n i-1 + φ n i-2 ). (4.20)
The flux can be written by the following form:

F 3 i+ = 1 2 (φ n i + φ n i+1 ) -(1 -Φ 3 i+ (1 -ν i )) (φ n i+1 -φ n i ) 2
,

F 3 i-= 1 2 (φ n i + φ n i-1 ) -(1 -Φ 3 i-(1 -ν i )) (φ n i -φ n i-1 ) 2 , (4.21)
and

Φ 3 i+ = 1 - 1 + ν i 3 (1 -r i+1/2 ), Φ 3 i-= 1 - 1 + ν i 3 (1 -r i-1/2 ), (4.22)
where r is the gradient ratio and r i+1/2 is defined as

r i+1/2 = (φ n i -φ n i-1 )/(φ n i+1 -φ n i ).
Following such successive corrections of the higher-order error terms, one can construct schemes of arbitrarily high N th -order of accuracy, whose numerical flux can be written in the generic form:

F N i+ = 1 2 (φ n i + φ n i+1 ) -(1 -Φ N i+ (1 -ν i )) (φ n i+1 -φ n i ) 2
,

F N i-= 1 2 (φ n i + φ n i-1 ) -(1 -Φ N i-(1 -ν i )) (φ n i -φ n i-1 ) 2 , (4.23)
where Φ N i+ is a function which gives N th -order of accuracy and it is also not continuous at the cell face. The modified scheme of third-order accuracy can be given by:

∂ t φ + u∂ x φ = u (1 -ν)(1 + ν)(ν -2) 4! ∆x 3 ∂ 4 x φ + O(∆x 4 ). (4.24)
To obtain a fourth-order scheme, Φ 4 can be given by:

Φ 4 i+ = Φ 3 i+ + 1 + ν i 3 ν i -2 4 (1 -2r i+1/2 + r i+1/2 r i-1/2 ), Φ 4 i-= Φ 3 i-+ 1 + ν i 3 ν i -2 4 (1 -2r i-1/2 + r i-1/2 r i-3/2 ). (4.25)
With the same idea, Φ at fifth-order of accuracy can be given by:

Φ 5 i+ = Φ 4 i+ - 1 + ν i 3 ν i -2 4
ν i -3 5 ( 1 r i+3/2 -3 + 3r i+1/2 -r i+1/2 r i-1/2 ), Φ 5 i-= Φ 4 i-- 1 + ν i 3 ν i -2 4
ν i -3 5 ( 1 r i+1/2 -3 + 3r i-1/2 -r i-1/2 r i-3/2 ). (4.26)
Φ at sixth-order of accuracy:

Φ 6 i+ =Φ 5 i+ + 1 + ν i 3 ν i -2 4
ν i -3 5
ν i + 2 6 • ( 1 r i+3/2 r i+5/2 - 4 r i+3/2 + 6 -4r i+1/2 + r i+1/2 r i-1/2 ), Φ 6 i-=Φ 5 i-+ 1 + ν i 3 ν i -2 4
ν i -3 5
ν i + 2 6 • ( 1 r i+1/2 r i+3/2 - 4 r i+1/2 + 6 -4r i-1/2 + r i-1/2 r i-3/2 ). (4.27)
Φ at seventh-order of accuracy:

Φ 7 i+ = Φ 6 i+ - 1 + ν i 3 ν i -2 4
ν i -3 5
ν i + 2 6
ν i + 3 7 • ( 1 r i+3/2 r i+5/2 - 5 r i+3/2 + 10 -10r i+1/2 + 5r i+1/2 r i-1/2 -r i+1/2 r i-1/2 r i-3/2 ), Φ 7 i-= Φ 6 i-- 1 + ν i 3 ν i -2 4
ν i -3 5
ν i + 2 6
ν i + 3 7 • ( 1 r i+1/2 r i+3/2 - 5 r i+1/2 + 10 -10r i-1/2 + 5r i-1/2 r i-3/2 -r i-1/2 r i-3/2 r i-5/2 ). (4.28)
with the principle of Taylor series analysis, the One-Step scheme can be developed up to eleventh-order accuracy [START_REF] Daru | High order one-step monotonicity-preserving schemes for unsteady compressible flow calculations[END_REF]. The case u < 0 can however be treated by symmetry relative to each cell interface.

Extension to multidimensional problem

To present the principle of the Strang directional splitting, we give an illustration in 2D. For a two-dimensional problem, the Level Set advection can be given by:

∂ t φ + u∂ x + v∂ y φ = 0. (4.29)
For each two time steps (t n to t n+2 ), according the principle given in Section 3.1.1, the Strang directional splitting strategy can be expressed as:

               φ n+ i, j = φ n i, j -ν n i,x (F i+, j (φ n ) -F i-, j (φ n )), φ n+1 i, j = φ n+ i, j -ν n i,y (G i, j+ (φ n+ ) -G i, j-(φ n+ )), φ n+1+ i, j = φ n+1 i, j -ν n+1 i,y (G i, j+ (φ n+1 ) -G i, j-(φ n+1 )), φ n+2 i, j = φ n+1+ i, j -ν n+1 i,x (F i+, j (φ n+1+ ) -F i-, j (φ n+1+ )), (4.30)
where F and G are the flux in different directions given by the OS type schemes. In such a way, the second-order accuracy of directional coupling is recovered every two time-steps if F and G does not comute.

Numerical validation of Level Set advection

In order to validate the interface description method, a significant amount of simulation cases have been devised over the years of development. Not only can these simulations quantify the performance of the methods, but also help to put the numerical schemes into practice. Simulation cases called "sin 4 ", "circle translation", "Zalesak disk" [START_REF] Zalesak | Fully multidimensional flux-corrected transport algorithms for fluids[END_REF] and "Severe interface deformation" are the most often encountered in the literature, so logically these cases should be implemented to validate our approach. The "sin 4 " is devoted to test the implementation of numerical schemes even sin 4 is not a signed distance function. With this test case, we can study the convergence with a smooth initial function. The circle translation problem are applied for testing the simple function advection with a periodic boundary condition. The test cases "Zalesak disk" and "Severe interface deformation" are calibrated for particular problems, the generation of the ligament in the case of "Severe interface deformation" and the contact angle transport in the case of Zalesak disk. Numerical results obtained from the OS family are compared to those from more classical schemes: HOUC-RK and WENO-RK.

To evaluate the performance of these numerical schemes, we use two mean errors: global L 1 error computed on the entire domain and local L 1 error computed in a narrow band of width 3∆x (three times of stencil size) around the interface. This narrow band is of higher interest, as it largely influences the interface precision and the computation of geometric quantities which are important for the hydrodynamic solver. Computations are performed on a personal computer with a 2.80 GHz Intel(R) i7-7600U CPU and use double floating-point values (64 bits) without multi-core parallelism.

Interface description CPU time t cpu is measured based on an average of at least 100 simulations. Presented methods are implemented in C++.

Test case: sin 4 (πx)

This test case is a one-dimensional advection of the function sin 4 (πx) with a periodic boundary condition in the domain of [-1, 1], as presented in Fig. 4.1. This problem can help us to test the implementation of numerical schemes, as the derivatives of the function sin 4 (πx) are always continuous, and the order of accuracy should reach the theoretical value.

The velocity of advection u equals to 1, the initial scalar function φ is defined as:

φ (x,t = 0) = sin 4 (πx). (4.31)
For all numerical schemes, the CFL number equals to 0.5. The numerical results of the separate and coupled time-space approaches are presented in Table 4.1 and Table 4. L 1 is measured at t = 1 since the period of this function equals to 1. Comparing the WENO5-RK3 scheme and the HOUC5-RK3 scheme, we can find that the HOUC5-RK3 scheme gives a better result with the rough meshes. Since the errors of WENO5-RK3 scheme is two times higher than those of HOUC-RK scheme, the WENO5-RK3 scheme converges faster, its numerical results are a little better than those given by the HOUC5-RK3 scheme for the fine meshes. The order of convergence for the separate time-space approaches is always limited to around 3. As u equals to 1 and CFL = 0.5, ∆t and ∆x are of the same order of magnitude. The errors of the separate time-space approaches are dominated by the temporal discretization of 3 rd -order of accuracy.

WENO5-RK3 HOUC5-RK3 HOUC7-RK3 Mesh L 1 error L 1 order L 1 error L 1 order L 1 error L 1 order 20 × 20 8.28e-2 - 4.58e-2 - 3.58e-2 - 40 × 40 1.
The schemes of OS families presented in Table 4.2 are much more precise than the separate time-space approaches presented in Table 4 4 (πx) by using the coupled time-space approaches Comparing numerical results of One-Step family to those of the separate time-space approaches, we can observe a significant improvement in terms of accuracy.

Test case: Circle advection

In order to study the convergence in two dimensions, here we set the initial function as φ (x i ) = 2 ∑ i=1 (x i -0.5) 2 -0.2 whose zero-level is a circle of radius 0.2 centered at (0.5, 0.5). This circle is advected in a simple stationary field, non-aligned with Cartesian mesh:

   u = 1, v = 1. (4.32)
The computation domain is

(x × y) ∈ [0, 1] × [0, 1]
with a periodic boundary condition in each direction. The CFL number is set to 0.5 and the errors are measured at t = 100 after 100 translations over the diagonal of the domain. As presented in the left panel of Fig. 4.2, all numerical schemes converge, but their order of accuracy of the global error is limited to around 1.6 even with a 7 th order accurate scheme. Although φ is a continuous distance function, its gradient is not always continuous.

Such discontinuities of ∇φ which exist at the boundary and the circle center for the initial function can lower numerical precision of high-order scheme. If we are interested in errors close to the interface where there is no discontinuity, local convergence is much faster than the global errors. In the right panel of Fig. 4.2, we can observe that for the coupled time-space approaches, they reach and almost exceed the theoretical order of accuracy. But for separate time-space, they are limited to an order of accuracy around 3. Numerical errors of separate time-space approach are constrained by O(∆t 3 ) which is the order of accuracy of the Runge-Kutta scheme while errors of OS scheme are constrained by O(∆x N ). By reducing time-steps much lower than the criterion of stability condition with RK schemes, we could recover the theoretical order of accuracy in time but at a significantly larger computational cost. When considering the errors versus the computation time, the OS scheme has the best performance. As presented in Fig. 4.3 where the error is plotted versus the CPU time, curves of both OS5 and OS7 are always below other schemes (for same order of accuracy) which signifies that OS schemes can give better accuracy for the same computational effort compared to separate time-space approaches.

-1.2 -0.9 -0.6 -0.3 0 As presented in Fig. 4.4, the advected zero-level contour of separate time-space approaches is deformed and becomes an ellipse while the zero-level contour of coupled time-space approach remains almost a circle as we paid attention of the operator symmetry in the Strang splitting procedure.

x/R -1.2 0.9 -0.6 -0.3 0 0.3 0.6 0.9 1.2 y/R OS5 WENO5-Rk3 HOUC5-RK3 Exact 1.2 0.9 0.6 0.3 0 x/R OS7 WENO5-RK3 HOUC7-RK3 Exact

Test case: Zalesak disk

This test case consists in the rotation of a rigid body (a disk with a rectangular slot located inside as the exact zero-level contour as presented in Fig. 4.5). This test presented by [START_REF] Zalesak | Fully multidimensional flux-corrected transport algorithms for fluids[END_REF] is demanding since it could reveal excessive numerical diffusion of schemes, as it presents some very strong local gradients around the slot. The computational domain is a square [0, 100] × [0, 100]; the disk of radius r = 15 is centered on the coordinates [START_REF] Idelsohn | The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves[END_REF][START_REF] Murrone | A five equation reduced model for compressible two phase flow problems[END_REF], the slot is a rectangle [START_REF] Aslam | A partial differential equation approach to multidimensional extrapolation[END_REF][START_REF] Dellacherie | Analysis of godunov type schemes applied to the compressible euler system at low mach number[END_REF] located on the vertical diameter of the disk, on its lower part. The rotating velocity field is stationary and defined as

   u = π(50 -y)/314, v = π(x -50)/314. (4.33)
Disk performs in this field one complete rotation in a dimensionless time t = 628, allowing the superimposition of the numerical solution with its initial state and thus a direct comparison with a subsequent error evaluation. The CFL number is set to 0.5. Fig. 4.5 shows the results of simulation after a full rotation of disk with different schemes with a mesh of 50 × 50 that corresponds to 2.5 cells in the slot width. The dissipation effects of numerical schemes are shown by the rounding angles of the slot corners. We compare the zero-level after a full rotation with the initial contours. We can observe that the WENO5-RK3 method smears the slot severely. For HOUC-RK and OS schemes, corners of the slot are significantly better resolved. When we look at the bottom of the slot, both HOUC-RK and OS schemes have almost the same performance (when comparing schemes at the same order of accuracy of spatial discretization). But when carefully considering the top of the slot, we can observe that this area is slightly better resolved by OS schemes (closer to initial solution).

To quantify the accuracy of these numerical schemes, we use global and local L 1 (narrow band of 3∆x) errors. Results are presented in Fig. 4.6. We can observe that both global and local errors of disk advection decrease with mesh refinement but that rate of convergence of global errors also decreases with mesh spacing. It could be explained by the complex geometry of the disk that exhibits sharp corners that lower the accuracy of numerical schemes. For the same type of scheme, the 7 th order scheme does not significantly improve solution compared to a 5 th order scheme. When comparing methods, we can observe that the WENO5-RK3 scheme results are always less accurate than other

4⋅10 -2 6⋅10 -2 8⋅10 -2 1⋅10 -1 2⋅10 -1 4⋅10 -1 32 64 128 256 global L 1 error number of grid points WENO5-RK3 HOUC5-RK3 OS5 HOUC7-RK3 OS7 1⋅10 -2 2⋅10 -2 4⋅10 -2 1⋅10 -1 2⋅10 -1 4⋅10 -1

Test case: Severe interface deformation

Contrary to previous tests, the velocity field selected in this problem induces the deformation of the initial Level Set condition. Starting from a circle shape, interface deforms itself in the generation of a long ligament. This case is of high interest since ligament stretches itself indefinitely as time advances: numerical simulation will face an under resolution situation. Numerical methods must be accurate enough to maintain ligament even when its width approaches cell size, and by consequence to maintain initial mass. A circle of radius r = 0.15 is initially centered at point (0.5, 0.75) inside a square domain [0, 1] 2 . The stationary rotating velocity field is defined by the potential function:

P = 1 π sin 2 (πx)sin 2 (πy), (4.34)
so that u, v components are defined by

   u(x, y) = -2sin 2 (πx)sin(πy)cos(πy), v(x, y) = 2sin 2 (πy)sin(πx)cos(πx).
The CFL number is set to 0.5. Interface is largely deformed at time t = 3 as presented in Fig. 4.7. The initial circle rolls itself around a central point and transforms into a long ligament. When comparing to the analytical zerolevel contours, we can observe that the WENO5-RK3 scheme does not capture interface as much as other schemes do. Both OS and HOUC-RK schemes with the same order of accuracy in space have equivalent interface description accuracy. In this test case, the velocity field verifies a free divergence. To assess the mass conservation property of each scheme, temporal surface evolution is computed and presented in Fig. 4.8. Compared to other schemes, surface conservation is less accurate for the WENO5-RK3 scheme. For all schemes, this surface conservation converges as the resolution is increased which will guarantee mass conservation of fluids. For lower resolution, poor mass conservation is a known drawback of the Level Set method and can be circumvented by coupling to the Volume Of Fluid method [START_REF] Ménard | Coupling level set/VOF/ghost fluid methods: Validation and application to 3D simulation of the primary break-up of a liquid jet[END_REF].

The redistancing

During the advection of the Level Set function, if the velocity field does not impose a rigid translation or rotation, the different contour lines are advected differently, as the local velocity gradient is not uniform: they are no more parallel and they no more represent a constant distance from the interface as shown in Fig. 4.9, |∇φ | ̸ = 1, causing numerical issues particularly in the evaluation of the normal and the curvature [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF]. In order to circumvent these problems, an additional treatment called "redistancing" is required to impose the distance property to φ . The distance property can be determined as the Fig. 4.9 Evolution of some iso-contours of Level Set function for the test case at different times: Severe interface deformation. red lines: iso-contours; black lines: interface; blue region:

D 1 (t) = {x ∈ D | φ (x,t) > 0} .
viscosity solution of the following eikonal equation:

   |∇φ | = 1, sign(φ ) = sign(φ 0 ). (4.35)
where φ 0 is the Level Set function before redistancing. To impose the distance property to φ , there are well-known algorithms for solving (4.35) such as the fast marching method [START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF] and fast sweeping method [START_REF] Zhao | A fast sweeping method for eikonal equations[END_REF]. Although these methods are simple and efficient, they are limited to first-order accuracy which can affect the accuracy of computing important geometric quantities. For this reason, here we introduce two different high-order numerical procedures for redistancing: PDE (Partial differential equation) and Hopf-Lax formula based redistancing procedures.

PDE based redistancing procedure

In [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF], Sussman et al. developed an algorithm for the Level Set redistancing problem. The main idea is that the zero contour line is correct, and the others should be corrected. The steady-state of the following advection equation solved into a fictitious time is the corrected Level Set:

   ∂ φ ∂t ′ = sign(φ 0 (x x x) 1 -|∇ φ | , φ (x x x,t ′ = 0) = φ 0 (x x x,t). (4.36)
The steady solution is a distance function |∇ φ | = 1, and no modification should be imposed on the interface. This equation can be solved with the extension of the WENO scheme done in Jiang and Peng [START_REF] Jiang | Weighted eno schemes for hamilton-jacobi equations[END_REF] to the Hamilton-Jacobi equations. Numerical resolution of (4.36) can be given by:

∂ φ ∂t ′ = H(φ 0 , φ , φ + x , φ - x , φ + y , φ - y ), (4.37)
where φ + x , φx , φ + y and φy are spatial derivatives of φ . H is the Hamiltonian and can be presented as a Godunov type flux:

H =    [max((∂ x φ + ) -, (∂ x φ -) + )] 2 + max((∂ y φ + ) -, (∂ y φ -) + )] 2 -1 if φ > 0, [max((∂ x φ + ) + , (∂ x φ -) -)] 2 + max((∂ y φ + ) + , (∂ y φ -) -)] 2 -1 otherwise, (4.38)
where (a) + = max(a, 0) and (b) -=min(b, 0). H is presented as a Godunov type flux. We should mention that the time-space coupled approach (One-Step) is difficult to adapt to (4.36), therefore, the spatial discretization of φ is given by (4.3) with the scheme of WENO5 as presented in precedent section of Level Set advection. (4.36) is then updated using:

φ l+1 i, j = φ l i, j -sign(φ 0 i, j )∆t ′ H, (4.39)
where sign(φ 0 i, j ) can be defined as:

sign(φ 0 i, j ) = φ 0 i, j (φ 0 i, j ) 2 + ∆x 2 . (4.40)
∆t ′ = 0.3∆x is a fictitious time step and proportional to ∆x. After certain number of iterations, the numerical results will converge. We set a global stopping criterion as given in [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF]:

E PDE = ∑ N | φ l+1 i, j -φ l i, j | N < ∆t ′ ∆x 2 . (4.41)
Here, N is the total cell number, and E PDE is a global average increment between two iterations.

Hopf-Lax formula based redistancing procedure

In [START_REF] Lee | Revisiting the redistancing problem using the hopf-lax formula[END_REF], Lee et al. proposed an algorithm for redistancing based on the Level Set method and the Hopf-Lax formula. We use the Hamilton-Jacobi formulation of the eikonal equation:

   ∂ ∂t ′ φ (x x x,t ′ ) + |∇ φ (x x x,t ′ )| = 0, φ (x x x, 0) = φ 0 (x x x). (4.42)
To know the distance to the interface, we simply need to know at which time t ′ the zero-isocontour of φ has progressed to the point x x x. For a point x x x with a negative initial value φ 0 (x x x), the time t ′ of propagation is also negative. In order to find t ′ (x x x) satisfying φ (x x x,t ′ (x x x)) = 0, here we follow the work of Royston et al. [START_REF] Royston | Parallel redistancing using the hopf-lax formula[END_REF]: the problem is treated as root finding by using secant method. At the first time, we consider the initial function φ 0 is convex.

Secant method for roots for φ (x x x i i i ,t ′ ) = 0 To find the correct t ′ (x x x i i i ), we use the secant method to solve for the root and we use the iterative update:

t ′, k+1 = t ′, k -φ (x x x i i i ,t ′ ) t ′, k -t ′, k-1 φ (x x x i i i ,t ′, k ) -φ (x x x i i i ,t ′, k-1 ) . (4.43)
Here, we take t ′,0 = 0 and the φ (x x x i i i , 0) = φ 0 (x x x) as initial guess, and set t ′,1 = ε∆x. In practice, we take ε equal to 0.5 and set ∆t ′ max = ±5∆x.

Hopf-Lax formula for φ (x x x i i i ,t)

To complete the iterative update for secant method, here we use the Hopf-Lax formula to compute

φ (x x x i i i ,t ′, k ) [64]: φ (x x x i i i ,t ′, k ) = min y y y∈R n φ 0 (y y y) + t ′, k H * x x x i i i -y y y t ′,k , (4.44)
where or equivalently φ (x

H * is the Fenchel-Legendre Transform of H * = ∥•∥ H * (x x x) =    0 |x x x| ≤ 1, ∞ otherwise, (4.45)
x x i i i ,t ′, k ) = argmin φ 0 (y y y), ∀ y y y ∈ B(x x x i i i ,t ′,k ), (4.46)
while B(x x x i i i ,t ′,k ) is the ball of radius t ′, k around the point x x x i i i . Then this problem becomes finding the minimum absolute value of the initial φ 0 over a ball. Here we follow the work of Royston et al. [START_REF] Royston | Parallel redistancing using the hopf-lax formula[END_REF], we use the projected gradient descent. Using the argmin y y y k-1 of t ′, k-1 as the initial guess of φ 0 over the ball B(x x x i i i ,t ′, k ), we can update the approximation iteratively from ȳ y y j+1 k = y y y j k -ε∆x ∇φ 0 (y y y j k ) |∇φ 0 (y y y j k )| sign(φ 0 (y y y j k )),

y y y j+1 k = PROJ(ȳ y y j+1 k ), (4.47)
where

PROJ(ȳ y y k ) =    y y y k |x x x i i i -y y y k | ≤ t ′, k , x x x i i i -t ′, k x x x i i i -y y y k |x x x i i i -y y y k | otherwise. (4.48)
Here, we take the same ε as the secant method and the function sign(φ 0 (y y y j k )) is the sign of the initial value of φ 0 (-1 or 1). Here, φ 0 (y y y) is defined by its interpolated values for grid nodes as

φ 0 (y y y) = ∑ i φ 0 i N i (x x x)
where N i is the bi-cubic interpolation kernel associated with grid node x x x i i i . So the gradients ∇φ 0 (y y y) can be expressed by: ∇φ 0 (y y y) = ∑ i φ 0 i ∇N i (y y y). In order to give a criterion to the convergence, here we set a tolerance number tol = 10 -8 . When the absolute value of φ (x x x i i i ,t ′, k ) is less than tol, the redistancing procedure for point x x x i i i converges, the regulated value for x x x i i i equals to ±t ′, k .

As each point can be updated individually, the Hopf-Lax approach naturally tends itself to narrow banding strategies.

Non-convexity of initial function φ 0

When we use the projected gradient descent to find the minimum value over a ball B(x x x i i i ,t ′, k ), when the initial function is not convex, the projected gradient descent will only converge to a local minimum, as presented in Fig. 4.10. The initial function φ 0 (x, y) = 4x 4x 2 + 4y 4y 2 has a local minimum at point (0, 0) for points with a negative value of φ 0 as presented in Fig. 4.10a. When we use the projected gradient descent method, we can find that there are some points that converge to this false local minimum, thus the result shows a terrible discontinuity as presented in Fig. 4.10c. Failure to converge to a global minimum can lead to large errors in the approximation of φ (x x x i i i ,t ′, k ). In order to resolve this problem, here we add a procedure called "Check procedure".

In [START_REF] Royston | Parallel redistancing using the hopf-lax formula[END_REF], Royston et al. solve (4.46) multiple times with different initial guesses y y y j k and then take the minimum over these solutions to come up with a final answer that is likely to be close to a global minimum. As this procedure should solve the projected gradient descent for every initial guess, it costs too much. Here, we only take the argmin y k-1 of t ′, k-1 as the initial guess, and use the projected gradient descent procedure to find a minimum. When the projected gradient descent procedure finds a minimum in the ball of B(x x x i i i ,t ′, k ), we prioritize random guesses in the shell B(x i ,t ′, k ) -B(x i ,t ′, k-1 ) between balls of t ′, k-1 and t ′, k as shown in Fig. 4.12. Then we compare these guesses with the minimum found by the projected gradient descent procedure to obtain the global minimum in the ball of B(x i ,t ′, k ).

In practice, on the order of one guess per grid cell in the shell

B(x i ,t ′, k ) -B(x i ,t ′, k-1
) is sufficient to find a reliable global minimum. With these random guesses, a more accurate global minimum can be found. Numerical results with check procedure are presented in Fig. 4.10d, the discontinuity disappears, this phenomenon means that the check procedure could help us to find the global minimum and compute the correct distance to the interface. The complete algorithm of the Hopf-Lax based method can be found in Algorithm 1.

y k y k-1
x i Fig. 4.12 Representation the check procedure. Blue dash line: ball of B(x i ,t ′, k ); red dash line: ball of B(x i ,t ′, k-1 ); blue point: argmin of t ′, k given by projected gradient decent procedure; red point: global argmin of t ′, k-1 ; black point: the initial point x i ; green points: random guesses of argmin of t ′, k ; gray area: shell B(x i ,t ′, k ) -B(x i ,t ′, k-1 ) between balls of t ′, k-1 and t ′, k . When d dt ′ φ (x x x i i i ,t ′ ) = 0, we can just set ∆t ′ = ±∆t ′ max , and jump out this non-convex area with the help of the check procedure. Then based on the projected gradient descent procedure, we can quickly find the correct time t of propagation or the regulated φ value as presented in Fig. 4.11. The redistancing of the Zalesak disk [START_REF] Zalesak | Fully multidimensional flux-corrected transport algorithms for fluids[END_REF], points in a narrow band of width 4∆x are well defined, the rest points have a value of ±4∆x.

Algorithm 1: Hopf-Lax based method for redistancing

while | φ (x x x i i i ,t ′, k )| > tol do // secant method: start ∆t ′ = -φ (x x x i i i ,t ′, k ) t ′, k -t ′, k-1 φ (x x x i i i ,t ′, k ) -φ (x x x i i i ,t ′, k-1 )
; k ); y y y j k = y y y j+1 k end φ P jd (x x x i i i ,t ′, k ) = φ 0 (y y y j k ); // projected gradient descent method: end // check procedure: start random guesses φ 0 (y

if |∆t ′ | > ∆t ′ max or φ (x x x i i i ,t ′, k ) = φ (x x x i i i ,t ′, k-1 ) then if φ (x x x i i i ,t ′, k-1 ) > 0 then ∆t ′ = ∆t ′ max ; else ∆t ′ = -∆t ′ max ; end end t ′, k-1 = t ′, k ; t ′, k = t ′, k + ∆t ′ ; //
y y ′ ) ∀ y y y ′ ∈ B(x x x i i i ,t ′, k ) -B(x x x i i i ,t ′, k-1 ); φCk (x x x i i i ,t ′, k ) = min(φ 0 (y y y ′ )); φ (x x x i i i ,t ′, k ) = min( φ ck (x x x i i i ,t ′, k ), φ P jd (x x x i i i ,t ′, k ));
// check procedure: end end

Numerical results of Level Set redistancing

In this section, we will present some numerical results based on these two redistancing algorithms presented in the precedent sections. Three test cases are employed in this section. In the first two test cases "Circle" and "Square", the initial functions are convex and simple, the check procedure is not necessary to implement, these two cases are implemented to show the advantage of the Hopf-Lax based method for initializing simple functions. While the test case "Ellipse" is closer to the practical application where the check procedure should be activated. We examine the performance of the Hopf-Lax method with the check procedure and compare it with the PDE based method. Let the procedure converge and end automatically by the stopping criterion. To evaluate these two methods, here we define the L 2 global error L G 2 , the L 2 local error L L 2 and the computation time t cpu . The local error is defined in a narrow band of width 3∆x around the interface.

Test case: Circle

In this test case, we set the initial function φ 0 (x than PDE based method, while the computational time t cpu is always one order of magnitude smaller than the PDE based method. Hopf-Lax method is accurate for the locating interface as it can reach to an error of 10

x x i i i ) = 1 2 ( 2 ∑ k=1 x 2 k -( 1

Test case: Square

In this test case, we perform redistancing to the initial φ 0 (x, y), whose zero-level is a square with side length 0.5 centered on (0.5, 0.5). Its computation domain is set to be [0, 1] × [0, 1], the graphs and The same problem is also presented in Fig. 4.14b. For the Hopf-Lax based method, for the center point, the initial direction of gradient descent is difficult to choose. Finally, no matter which direction it chooses, it can converge to the interface and find a correct value, as the interface is centrally symmetric. Numerical errors and computational time are presented in

Test case: Ellipse

In this test case, the initial function is set to be φ 0 (x

x x i i i ) = 1 2 ( 2 ∑ k=1 x 2 k B 2 k -1
) with a zero level of an ellipse.

We set B 1 = 0.8 and B 2 = 0.4 for the test as the representation and isocontours given in Fig. 4.17 We should mention that for this test case that the check procedure should be activated. Because the path given by the projected descent method is not always the fastest even if the initial function is convex. The first derivatives might vary a lot in the domain, direction of the fastest descent direction may be followed by an area of small gradients, so choosing the projected gradient descent direction only by the local gradient is not always correct. The check procedure gives several random guesses in the ball B(x x x i i i ,t k ), a relatively correct global minimum can be found. Results after redistancing are presented in Fig. 4.18. Form the graphs and isocontours, we cannot observe the difference. With the errors and CPU time presented in Table 4.5, we can find that the Hopf-Lax based method doesn't have an advantage compared with PDE based method as the test case "Circle" and "Square" on CPU time since to ensure that the global minimum is not omitted, a large number of random guesses are employed. The local and global errors given by Hopf-Lax based method is always one or two orders of magnitude smaller than the PDE based method and can reach an accuracy of 10 From these test cases, we can get a conclusion that the PDE based method is a classical method and always gives a stable and relatively accurate result, the Hopf-Lax based method can give a more accurate result and has a good performance for simple and convex initial function. When the function becomes complex and non-convex, its computation time increases. But generally, it can promise excellent accuracy for the interface, as it uses the projected gradient descent method to find the closest point at the interface that imposes zero change to the interface.

Effect of redistancing on Level Set advection

In the precedent section, several different numerical schemes for Level Set advection and redistancing are presented and studied. As in practice, the redistancing procedure should be imposed on an advected field of Level Set to grantee the signed distance property, here we want to study the effect of redistancing on Level Set advection. In the test case "Severe interface deformation", the redistancing is necessary to improve the numerical results.

As presented in the section of Level Set advection, the WENO5-RK3 scheme has a poor performance while and the OS5 scheme can relatively well predict the interface motion. In order to observe the effect of redistancing procedure, here we present numerical results of WENO5-RK3 and OS5 Level Set advection schemes with redistancing procedure in Fig. 4.19, we can find the redistancing procedure can help us to suppress the break-up at the end of the ligament. Comparing these two different redistancing methods, the Hopf-Lax based method has a visible effect of increasing the length. However, to keep this length in a coarse mesh, the width of the ligament is numerically enlarged, leading to a significant increase in surface as presented in Fig. 4.21. For the PDE based method, we can find that the length is slighted increased, but smaller than the theoretical solution. The surface evolution trend is not obvious as presented in Fig. 4.20.

The surface evolution of different couplings of numerical transport and redistancing schemes are presented in Fig. 4.20 and Fig. 4.21. Both of these two different redistancing methods suffer the problem of surface conservation. In terms of surface conservation, all of these numerical methods converge. To obtain a precise method, an alternative option might be adaptive mesh refinement.

Conclusion on Level Set advection and redistancing

With numerical tests for Level Set advection, we can find that the OS5 scheme has an important advantage in terms of accuracy and efficiency, therefore, the retained scheme for Level Set advection is the OS5 scheme. Based on the discretization of the OS scheme in Section 4.1.2, we deduce an efficient implementation of the OS scheme. The discrete form of the one-dimensional OS scheme recovers a simple expression:

φ n+1 i = φ n -u∂ N,± x,t φ + O(∆x N , ∆t N ), (4.49)
where ∂ N x,t is a numerical discretization of first derivative to ensure odd N th order of accuracy in time and space, and ∂ N,+

x,t and ∂ N,- x,t are upwind choices applied for positive and negative u respectively. An arbitrary odd N th order of accuracy can be obtained by using the odd order accuracy discretization:

∂ N,± x,t φ = 1 ∆x N ∑ k=0 C N,± k φ i±k∓ N+1 2 .
In the present work, N = 5 order of accuracy should be employed. Coefficients C for the fifth accuracy order are:

C 5,± 0 = ±(ε 4 + ε 5 ), C 5,± 1 = ±(ε 3 -4ε 4 -5ε 5 ), C 5,± 2 = ±(-1 + ε 2 -3ε 3 + 6ε 4 + 10ε 5 ), C 5,± 3 = ±(1 -2ε 2 + 3ε 3 -4ε 4 -10ε 5 ), C 5,± 4 = ±(ε 2 -ε 3 + ε 4 + 5ε 5 ), C 5,± 5 = ±(-ε 5 ).
The ε coefficients based on the CFL condition ν are the followings:

ε 2 = 1 -ν 2 , ε 3 = ε 2 1 + ν 3 , ε 4 = ε 3 ν -2 4 , ε 5 = ε 4 ν -3 5 . (4.50)
For the Level Set redistancing problem, compared with the PDE based method, the Hopf-Lax method has an important advantage on simple geometry. When the interface is distorted and the Level Set function is non-smooth and non-convex, a check procedure should be implemented. Based on random guesses, this procedure is time-consuming and difficult to parallelize. With coarse meshes, the problem of surface conservation is encountered for both these two redistancing methods. To overcome this problem, a common strategy is to refine the mesh. Besides, the PDE based method is more adapted to the Level Set advection procedure and easy to parallelize. As a reason, the retained scheme for Level Set redistancing is the PDE based method. Numerical results of "Severe interface deformation" presented in Section 4.5 show that the coupling of OS5 advection scheme and PDE based redistancing method can give a prediction with satisfactory with a resolution of 128 × 128.

Curvature estimation

To investigate capillary effects, the estimation of interface curvature is mandatory. As present in Section 1.2.3, the Level Set method has an important advantage on interface geometry estimation. The classical method of interface curvature estimation is to discretize (1.2) with the Level Set function.

The estimation of interface curvature in the framework of VOF method is not as direct as the Level Set method, the Height function method is then developed for curvature estimation [START_REF] Cummins | Estimating curvature from volume fractions[END_REF]. This method is then applied to Level Set method [START_REF] Abadie | On the combined effects of surface tension force calculation and interface advection on spurious currents within volume of fluid and level set frameworks[END_REF] even if the interface curvature could be derived from (1.2).

In this section, we will introduce these two methods and compare their numerical performance.

Classical Level Set (CLS) method

As presented in (1.2), the interface curvature can be deduced directly from the signed distance function φ and can be evaluated as:

κ = ∇ • (∇φ ) |∇φ | = 2φ x φ y φ xy -φ 2 x φ yy -φ 2 y φ xx (φ 2 x + φ 2 y ) 3/2 . ( 4.51)
The first and second derivatives are discretized with a centered second-order scheme. The curvature obtained with (4.51) is the local curvature rather than the curvature at the interface which depends on its exact position. To estimate the curvature at interface, κ should be evaluated at the interface between two adjacent cells j and k separated by the interface, its expression can be given as:

(κ) jk = κ j |φ k | + κ k |φ j | |φ j | + |φ k | . (4.52)

Height function method

The Height function method is based on the collection of heights h (parallel distances to the interface) of actual and neighboring cells and the curvature is then computed using finite-difference operators on those heights:

κ = h ′′ (1 + h ′2 ) 3/2 . ( 4.53)
In the framework of VOF, heights h are computed by integrating the liquid volume fraction. With Level Set method, h are distances to the interface in a certain direction as presented in Fig. 4.22.

Based on choices of direction, here we introduce two sub-methods: pseudo-normal direction method (PND), and rotated normal direction (RND).

Pseudo-normal direction (PND) method

The pseudo-normal direction [START_REF] Cummins | Estimating curvature from volume fractions[END_REF] is the Cartesian direction x, y, or z with the largest component of the interface normal vector as presented in Fig. 4.22 for κ 1 estimation. The largest component of the normal vector to the interface is in the direction of y, so the heights are the distance to the interface in direction of y.

In order to discretize (4.53), we set three consecutive heights h n with n = (0, ±1), centered around h 0 . The curvature is then calculated using finite difference operators on those heights. The discrete form of (4.53) could be written as [START_REF] Bornia | On the properties and limitations of the height function method in two-dimensional cartesian geometry[END_REF]:

κ = 2d(h -1 -2h 0 + h +1 ) (d 2 + ∆h 2 (0,-1) )(d 2 + ∆h 2 (0,+1) )(4d 2 + ∆h 2 (-1,+1) ) h 0 Γ h -1 h +1 x y κ 1 d κ 2
Fig. 4.22 PND method Heights h representation for κ 1 estimation where d is the distance between adjacent heights and ∆h 2 (l,m) = (h lh m ) 2 . We employ the same method as the redistancing Hopf-Lax method to find the height or distance to interface in a certain direction. Compared to the redistancing procedure, we only need to modify the projected gradient descent and apply it in the direction n n n ′ following the largest component of the normal vector to the interface. The complete algorithm is available in Algorithm 2.

Rotated normal direction (RND) method

For highly under-resolved interface where κ∆x ≈ 1 as κ 2 presented in Fig. 4.22, an adequate number of well-defined heights in Cartesian directions may not be available. This problem could be resolved by defining h in the normal direction to interface [START_REF] Owkes | A mesh-decoupled height function method for computing interface curvature[END_REF], which is the RND method. Heights for κ 1 and κ 2 estimation are then defined in the normal direction to the interface as presented in Fig. 4. [START_REF] Daru | High order one-step monotonicity-preserving schemes for unsteady compressible flow calculations[END_REF].

Similar to the PND method, to estimate the curvature of a point x x x i i i , we compute different heights to the interface in its normal direction to the interface by the Hopf-Lax method as described in Algorithm 2 with n n n ′ = n n n(x x x i i i ), and take the same discrete form (4.54) to estimate the interface curvature. 4.23 RND method Heights h representation for κ 1 and κ 2 estimation Algorithm 2: Hopf-Lax based method for computing heights in any given direction n ′ while | φ (x x x i i i ,t ′, k )| > tol do // secant method: start

h 0 h -1 h +1 x y κ 1 n 1 d κ 2 Γ n 2 h 0 h -1 h +1 d Fig.
∆t ′ = -φ (x x x i i i ,t ′, k ) t ′, k -t ′, k-1 φ (x x x i i i ,t ′, k ) -φ (x x x i i i ,t ′, k-1 )
;

if |∆t ′ | > ∆t ′ max or φ (x x x i i i ,t ′, k ) = φ (x x x i i i ,t ′, k-1 ) then if φ (x x x i i i ,t ′, k-1 ) > 0 then ∆t ′ = ∆t ′ max ; else ∆t ′ = -∆t ′ max ; end end t ′, k-1 = t ′, k ; t ′, k = t ′, k + ∆t ′ ;
// secant method: end // projected directional descent method: start // in any given direction n n n ′ y y y j k = argmin φ (y y y), ∀y y y ∈ L(x x x i i i , n n n ′ ,t ′, k-1 ) ; while |y y y j kx x x i i i | < |t ′, k | do ȳ y y j+1 k = y y y j k -ε∆xn n n ′ sign(φ 0 (y y y j k )); y y y j+1 k = PROJ(ȳ y y j+1 k ); y y y j k = y y y j+1 k end φ P jd (x x x i i i ,t ′, k ) = φ 0 (y y y j k ); // projected directional descent method: end // check procedure: start random guesses φ 0 (y

y y ′ ) ∀ y y y ′ ∈ L(x x x i i i ,t ′, k ) -L(x x x i i i ,t ′, k-1 ); φCk (x x x i i i ,t ′, k ) = min(φ 0 (y y y ′ )); φ (x x x i i i ,t ′, k ) = min φ ck (x x x i i i ,t ′, k ), φ P jd (x x x i i i ,t ′, k ) ; // check procedure: end end Note: L(x x x i i i ,t ′, k
) is a set of points on the line passing point x x x i i i in the given direction n n n ′ with a distance to point x x x i i i less than |t ′, k |

Numerical results for curvature estimation

In this section, we will perform a test on interface curvature estimation for a 2D circle of radius R = 1. From (4.54), we can find that the interface curvature precision depends not only on the way that we compute heights but also on the distance between adjacent heights. With numerical experiments, we will study the influence of different Height function methods as well as the distance between heights d. From Fig. 4.24, we can find that all these methods converge and the classical Level Set method can reach theoretical second-order of accuracy, but Height function methods are limited to first-order of accuracy. Though we adopt three consecutive heights, in theory, second-order of accuracy could be reached. One possible reason for lowering the order of accuracy might be related to errors of height estimation.

Comparing PND method with RND method, we can get a conclusion that: with the same value of d, the RND method is able to measure a higher curvature, but L 1 error of the PND method is always slightly smaller than the RND method. As for the influence of the value of d, a smaller d corresponds to the ability to estimate higher curvature and the error is more important than a big height distance d.

Conclusion on interface curvature estimation

For the interface curvature estimation method, considering that the classical Level Set method is easy to implement and converges more quickly than Height function methods, the retained scheme for estimation is the classical Level Set method.

Summary

In this chapter, by comparing different numerical schemes to resolve the problem of Level Set advection, redistancing and curvature estimation, we made a choice for the employed methods: they are respectively OS5, PDE based redistancing method and the classical Level Set method.

With a precise tool to describe the interface, the next chapters will be devoted to develop a hydrodynamic solver with an accurate prediction for two-phase compressible flows.

Chapter 5

Numerical solver

In this chapter, we will present the numerical scheme to solve the global system within the Eulerian framework. Based on the acoustic-transport splitting, the global system is splitted into acoustic, transport and diffusion subsystems. By following the works of Chalons et al. [START_REF] Chalons | An all-regime lagrange-projection like scheme for the gas dynamics equations on unstructured meshes[END_REF][START_REF] Chalons | A large time-step and wellbalanced lagrange-projection type scheme for the shallow water equations[END_REF] on the Lagrangeprojection type scheme, we resolve the acoustic and transport subsystems with a well-balanced discretization for the bulk flow. The diffusion subsystem is discretized with a classical second-order centered scheme. The coupling across the interface is realized by the ghost-fluid method [START_REF] Fedkiw | A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF], an approximated Riemann solver that accounts for jump conditions is developed. With the presence of phase change, an explicit approximation of heat flux to the interface is introduced.

Definition of the approximate Riemann solver for the acoustic subsystem

By taking jump conditions across the interface Γ(t) into account, the acoustic subsystem is set by accounting for phase change and the surface tension boundary condition. The jump conditions are transformed into source source term at the interface with Dirac mass. We consider the following acoustic subsystem defined on the whole domain D:

t > 0, x ∈ D              ∂ t ρ + ρ∇ • u = ρS v j , ∂ t (ρu) + ρu∇ • u + ∇p = S c ∇φ + ρS v j u -ρ∇Ψ, ∂ t (ρE) + ρE∇ • u + ∇ • (pu) = S c u • ∇φ -ρu • ∇Ψ + (ρE + p)S v , ∂ t φ = 0, (5.1a) (5.1b) (5.1c) (5.1d) with S c = p Γ δ (φ ) /|∇φ |, S v = ∇ • ( u n Γ n)δ (φ ) , (5.2)
where δ (φ ) is the Dirac measure whose support is Γ(t). We present an approximate Riemann solution to solve the acoustic subsystem within the Eulerian framework.

Firstly, by omitting source terms related to external force and velocity jump, the acoustic subsystem (5.1) is a quasilinear system. In order to derive the resolution of this system, we will perform several approximations.

We notice that for a smooth solution (5.1), the pressure verifies: ∂ t p + (ρc) 2 τ∇ • u = 0. We thus perform a Suliciu-type relaxation adapted to vacuum [START_REF] Bouchut | Nonlinear stability of finite Volume Methods for hyperbolic conservation laws: And Well-Balanced schemes for sources[END_REF] of (5.1) by introducing a surrogate pressure π and considering following relaxed system by using non-conservative variables with τ = 1/ρ:

                           ∂ t τ -τ∇ • u = 0, ∂ t u + τ∇π = τS c ∇φ , ∂ t E + τ∇ • (πu) = τS c u • ∇φ , ∂ t π + a 2 τ∇ • u = χ(p -π), ∂ t a = 0, ∂ t φ = 0.
(5.3) a > 0 is the acoustic impedance and its definition will be given later. In regime χ → ∞, we formally recover (5.1). In our numerical solver context, we classically mimic χ → ∞ regime by enforcing at each time step π n i = p EOS (τ n i , e n i ) (given by (2.2)) and then solving (5.3) with χ = 0. Assuming a one-dimensional problem in x-direction, we have ρ(x,t)∂ t ≈ ρ(x,t n )∂ t , then if we define a mass variable m by dm(x) dx = ρ(x,t n ), we obtain up to a slight abuse of notation, a system that is written in vector form: ∂ t W + A∇ m W = 0, where

W =                 τ u π E a v w φ                 , A =                
0 -1 0 0 0 0 0 0 0 0 1 0 0 0 0 -S c 0 a 2 0 0 0 0 0 0 0 π u 0 0 0 0 -S c u 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

                . The matrix A is diagonalizable det(A -λ Id) = -λ 6 (λ 2 -a 2 ),
and its eigenvalues are: (-a, 0, 0, 0, 0, 0, 0, +a). Fields involved in this system are all linearly degenerated. Considering a discontinuity that propagates at celerity D, let b be the jump of the variable b across the cell interface. Following the Rankine-Hugoniot relationships we get the jump conditions:

                                     -D τ -u = 0, -D u + π -q φ = 0, -D π + a 2 u = 0, -D E + πu -q ′ φ = 0, -D a = 0, -D v = 0, -D w = 0, -D φ = 0, (5.4a)
(5.4b)

(5.4c)

(5.4d)

(5.4e)

(5.4f)

(5.4g)

(5.4h)

with D has the same value as eigenvalues of matrix A. q and q ′ are weights associated with the Dirac masses M and Mu, their definition will be given latter. D = 0 yields that u is an independent Riemann invariant for the wave λ = 0 and brings in the jump conditions π = q φ and a ∈ R in (5.4c) and (5.4e). Thus a is not necessary to be continuous across the wave λ = 0, and the jump condition on π should conform the jump relation (2.8).

The structure of the acoustic system (5.1) shares similarities with that has been studied in [START_REF] Chalons | A large time-step and wellbalanced lagrange-projection type scheme for the shallow water equations[END_REF]. This allowed [START_REF] Zou | An accurate sharp interface method for two-phase compressible flows at low-mach regime[END_REF] to propose an approximate Riemann solver W RP that is consistent in the integral sense with the full system (5.1) including the surface tension terms.

Let W T = (ρ, u, v, w, E, φ , π) where π acts as an approximated pressure value for the system, we briefly recall hereafter the definition of W RP and have

W RP (ζ ; W l , W r ) =                W l , if ζ < -a l , W * l , if -a l < ζ < 0, W * r , if 0 < ζ < a r , W r , if a r < ζ .
(5.5)

First the jump conditions across the wave λ = 0 imply that π * rπ * l = q(φ rφ l ) and u * = u * l = u * r . In W RP , the discrete expression of the interface pressure jump can be expressed as:

π * l -π * r = p lr (5.6)
where p lr = lim between Ω l and Ω r , this jump could be not zero when φ l φ r < 0. Then the jump conditions (5.4b) across the waves -a l and a r give:

a l (u * -u l ) +π * l -π l = 0, a r (u r -u * ) -π r -π * r = 0.
(5.7)

By solving (5.7) with (5.6), intermediate states W * l and W * r can be defined by

π * l =
(a l π r + a r π l ) a l + a r -a l a r (u ru l ) a l + a r + p lr a l a r + a l , (5.8a)

π * r =
(a l π r + a r π l ) a l + a r -a l a r (u ru l ) a la r -p lr a r a l + a r , (5.8b)

1 ρ * l = 1 ρ l + a r (u r -u l ) + π l -π r -p lr a l (a l + a r ) , (5.8c) 1 ρ
* r = 1 ρ r + a l (u r -u l ) + π r -π l + p lr a r (a l + a r ) , (5.8d
)

u * = u * l = u * r = a l u l + a r u r a l + a r - π r -π l + p lr a l + a r , (5.8e)
v * l = v l , v * r = v r , (5.8f)
w * l = w l , w * r = w r , (5.8g)

E * l = E l - π * l u * -π l u l a l , (5.8h
)

E * r = E r + π * r u * -π r u r a r , (5.8i)
where κ is the local value of interface curvature. The parameters a l and a r are user-chosen constants that need to be sufficiently large approximations for ρ l c l and ρ r c r so as to satisfy stability constraints for the numerical scheme. Indeed, according to Bouchut [START_REF] Bouchut | Nonlinear stability of finite Volume Methods for hyperbolic conservation laws: And Well-Balanced schemes for sources[END_REF], if one notes {b} + = max(b, 0), it is possible to show that defining a l and a r as follows if π rπ l + p lr ≥ 0,

       a l = ρ l c l + γ l +1 2 ρ l π r -π l + p lr ρ r c r + u l -u r + a r = ρ r c l + γ r +1 2 ρ r π l -π r -p lr a l + u l -u r + , if π r -π l + p lr ≤ 0        a r = ρ r c l + γ r +1 2 ρ r π l -π r -p lr c l ρ l + u l -u r + a l = ρ l c l + γ l +1 2 ρ l π r -π l + p lr a r + u l -u r + ,
ensures positivity for the density. In the present work we will consider two alternate choices: a first simple choice that allows to recover the classic Suliciu relaxation solver [START_REF] Bouchut | Nonlinear stability of finite Volume Methods for hyperbolic conservation laws: And Well-Balanced schemes for sources[END_REF][START_REF] Chalons | Navier-stokes equations with several independent pressure laws and explicit predictor-corrector schemes[END_REF][START_REF] Chalons | Relaxation approximation of the euler equations[END_REF][START_REF] Suliciu | On the thermodynamics of rate-type fluids and phase transitions. i. rate-type fluids[END_REF] by setting

a l = a r = K max(ρ l c l , ρ r c r ), (5.9)
where K > 1 is a constant and a second heuristic choice obtained by defining

a l = K ′ ρ l c l , a r = K ′ ρ r c r , K ′ = 1.1.
(5.10)

This approximate Riemann solver naturally takes the pressure jump into account. By following the work of Chalons et al. [START_REF] Chalons | A large time-step and wellbalanced lagrange-projection type scheme for the shallow water equations[END_REF], with the same wave structure, we modify this approximate solver by considering the external force terms and velocity jump in a consistent way.

Approximate Riemann solver for the bulk flows

When we want to keep the interface sharp and use the ghost-fluid method to realize the coupling across the interface, the surface tension and phase change do not have a direct impact on bulk flows. As a reason, in bulk flows, effects related to the surface tension and phase change are not taken into consideration. To use the pressure gradient to balance the external force, we modify the approximated solver W RP , the intermediate states W * l and W * r for the bulk flow approximate Riemann solver can be given by:

π * l = π * r = (a l π r + a r π l ) a l + a r - a l a r (u r -u l ) a l + a r + (a l -a r )(ρ∆Ψ) lr a l + a r , (5.11a
)

u * = u * l = u * r = a l u l + a r u r a l + a r - π r -π l -(ρ∆Ψ) lr a l + a r , (5.11b) 1 ρ
* l = 1 ρ l + a r (u r -u l ) + π l -π r + (ρ∆Ψ) lr a l (a l + a r ) , (5.11c) 1 ρ
* r = 1 ρ r + a l (u r -u l ) + π r -π l -(ρ∆Ψ) lr a r (a l + a r ) , (5.11d)
v * l = v l , v * r = v r , (5.11e)
w * l = w l , w * r = w r , (5.11f)

E * l = E l - π * l u * -π l u l a l , (5.11g
)

E * r = E r + π * r u * -π r u r a r . (5.11h)
To get a well-balanced scheme, the external force term is discretized at the cell face as:

(ρ∆Ψ) lr = ρ l + ρ r 2 (Ψ l -Ψ r ).
(5.12)

Here we recall that a well-balanced scheme aims at preserving discrete versions of some continuous equilibrium states.

Approximate Riemann solver at the interface

For the approximate Riemann solver at the interface, effects related the surface tension, phase change and external force effects should be taken into account. Considering a one-dimensional problem in the direction x, we define the discrete expression of the pressure and velocity jumps as:

p * l -p * r = p Γ lr , u * l -u * r = u • ∂ x φ Γ lr = u Γ lr , (5.13)
where

Z Γ lr = lim x→y x∈Ω l Z -lim x→y x∈Ω r
Z for y ∈ Γ(t n ), t n > 0 means the jump of the variable Z across the interface between Ω l and Ω r .

To eliminate numerical oscillations and preserve equilibrium for cases where the interface does not coincide with the cell face, the potential energy ρΨ is estimated by using the interface position as a reference that takes the density discontinuity into account. We want our scheme to strictly preserve the equilibrium steady solutions, that are the states satisfying (3.25). The approximated system of the Riemann problem being solved at the interface instead of at the cell face. To be consistent with the bulk flow, the intermediate pressure is then extrapolated constantly to the cell face by taking the potential energy into account to satisfy the condition (3.25). Then we have intermediate states at the interface:

π * l = a l π r + a r π l + p Γ lr a l + a l δ p l a l + a r - a l a r (u r -u l + u Γ lr ) a l + a r , (5.14a)
π * r = a l π r + a r π l -p Γ lr a r + a r δ p r a l + a r - a l a r (u r -u l + u Γ lr ) a l -a r , (5.14b) 1 ρ
* l = 1 ρ l + a r (u r -u l + u Γ lr ) + π l -π r -p Γ lr + ρ lr ∆Ψ lr a l (a l + a r ) , (5.14c) 1 ρ
* r = 1 ρ r + a l (u r -u l + u Γ lr ) + π r -π l + p Γ lr -ρ lr ∆Ψ lr a r (a l + a r ) , (5.14d)
u * l = a l u l + a r u r + a r u Γ lr a l + a r -π rπ l + p Γ lr a l + a r + ρ lr ∆Ψ lr a r + a l , (5.14e)

u * r =
a l u l + a r u ra l u Γ lr a l + a r -π rπ l + p Γ lr a l + a r + ρ lr ∆Ψ lr a r + a l , (5.14f)

v * l = v l , v * r = v r , (5.14g)
w * l = w l , w * r = w r , (5.14h)
E * l = E l - π * l u * -π l u l a l , (5.14i)
E * r = E r + π * r u * -π r u r a r , (5.14j)
with

∆Ψ lr = (Ψ l -Ψ r ), ρ lr = |φ l |ρ l + |φ r |ρ r |φ r -φ l | , δ p l = -ρ lr ∆Ψ lr + ρ l (a l + a r )∆Ψ lr 2a l , δ p r = ρ lr ∆Ψ lr -ρ r (a l + a r )∆Ψ lr 2a r .
We should note that, the approximation of the jump relation associated with surface tension has been given in Chapter 4, while the approximation of jump relation associated with phase change will be given latter in the discretization of heat transfer. (5.10) can also satisfy the condition of positivity of densities in (5.11) and (5.14) in low Mach regime, for more details please refer to the work of Bouchut [START_REF] Bouchut | Nonlinear stability of finite Volume Methods for hyperbolic conservation laws: And Well-Balanced schemes for sources[END_REF].

Finite volume discretization

We now recall the finite volume approximation for each of acoustic, transport and viscous effect systems. We suppose that the computational domain is discretized over Cartesian grid. If Ω j and Ω k are two cells of the mesh, we note ∂ Ω jk = Ω k ∩ Ω j . We define N ( j) = {k | ∂ Ω jk ̸ = / 0} and we note n jk the unit normal vector to ∂ Ω jk oriented from Ω j to Ω k . The center of the cell Ω j will be denoted by x j and for k ∈ N ( j) the center of ∂ Ω jk will be noted x jk = x k j as depicted in figure 5.1. With a slight abuse of notation, for a function x → b, when k ∈ N ( j) we will note b(x jk ) = lim

x→x jk x∈Ω j b(x), b(x k j ) = lim x→x jk x∈Ω k b(x).
We define N ± ( j) = {k ∈ N ( j) | ± φ n k φ n j > 0} and we also note We should note that, in the present work, the acoustic and transport subsystems are discretized with the finite volume method of first-order accuracy while the viscous and diffusion subsystems are discretized with the finite volume method of second-order accuracy.

∆x jk = x j -x k = ∆x k j , so that ∆x jk /2 = x j -x jk = x k j -x k and ∆x k j /2 = x k -x k j = x k j -x j . x k j x j cell j x k cell k

Approximation of the acoustic subsystem

By adapting the lines of Chalons et al. [START_REF] Chalons | A large time-step and wellbalanced lagrange-projection type scheme for the shallow water equations[END_REF] thanks to the rotational invariance of acoustic system, we use the approximate Riemann solver described in Section 5.1 to build a Finite Volume approximation of the acoustic system.

It reads

L j ρ n+ j = ρ n j , L j (ρu) n+ j = (ρu) n j - ∆t |Ω j | ∑ k∈N ( j) |∂ Ω jk |π * jk n jk -∆t{ρ∇Ψ} j , L j (ρE) n+ j = (ρE) n j - ∆t |Ω j | ∑ k∈N ( j)
|∂ Ω jk |π * jk u * jk -∆t{u * ρ∇Ψ} j ,

φ n+ j = φ n j , L j = 1 + ∆t |Ω j | ∑ k∈N ( j)
|∂ Ω jk |u * jk .

(5.15a)

(5.15b)

(5.15c)

(5.15d)

(5.15e)

With Cartesian mesh, terms related to external forces can be approximated as:

{ρ∇Ψ} j = 1 2 ∑ k∈N ( j) {ρ∇Ψ} jk , {u * ρ∇Ψ} j = 1 2 ∑ k∈N ( j)
u * jk n jk • {ρ∇Ψ} jk .

(5.16a)

(5.16b)

∀k ∈ N + ( j) :                      u * jk = n jk • (a j u # j + a k u # k ) a j + a k - π # k -π # j -(ρ∆Ψ) n jk (a j + a k ) , π * jk = a k π # j + a j π # k a j + a k - a j a k a j + a k n jk • (u # k -u # j ) + (a l -a r )(ρ∆Ψ) jk a j + a k , {ρ∇Ψ} jk = ρ n j + ρ n k 2|∆x jk | (Ψ k -Ψ j )n jk .
(5.17a) Z for y ∈ Γ(t n ), t n > 0 means the jump of the variable Z across the interface between Ω j and Ω k , ∀k ∈ N -( j). We take the same acoustic impedance definition as (5.10):

(5.17b) (5.17c) ∀k ∈ N -( j) :                              u * jk = n jk • (a j u # j + a k u # k + a k u Γ jk ) a j + a k - π # k -π # j + p Γ jk -(ρ∆Ψ) n jk (a j + a k ) , π * jk = a k π # j + a j π # k + a j p Γ jk + a j δ p jk a j + a k - a j a k n jk • (u # k -u # j + u Γ jk ) a j + a k , ∆Ψ jk = (Ψ j -Ψ k ), ρ jk = |φ j |ρ j + |φ k |ρ k |φ k -φ j | , {ρ∇Ψ} jk = ρ n j |∆x jk | (Ψ k -Ψ j )n jk , δ p jk = -ρ n jk ∆Ψ jk + ρ n j (a j + a k )∆Ψ jk 2a j , (5.18a)
a j = Kρ j c j , a k = Kρ k c k , K = 1.1. (5.19)
In this work, we employ an explicit update for the acoustic subsystem with # = n. However, The superscript # introduced can take another value, # = n+, then the acoustic scheme is implicit. Whatever the choice made, we can update the acoustic variables ρ , ρu and ρE.

A Courant-Friedrichs-Lewy (CFL) constraint on the time step choice for (5.15)-(5.19) pertaining to the acoustic effects reads

∆t acoustic |∂ Ω jk | |Ω j | max(τ j a j , τ k a k ) ≤ 1 2 , (5.20)
while surface tension effects are stable on the condition that ∆t verifies [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF]:

∆t capillary ≤ (ρ j + ρ k ) 4Πσ |Ω j | |∂ Ω jk | 3 , (5.21)
where Π is the Archimedes' constant.

Approximation of the transport subsystem

As a precise interface description is highly requested for two-phase flow simulations, the Level Set advection has been previously treated independently with a high-order scheme. We write the transport subsystem with conservative variables which causes conservative terms to appear:

∂ t b + ∇ • (bu) -b∇ • u = 0, b = (ρ, ρu, ρv, ρw, ρE) T .
We discretize this system by using result of the previous acoustic step and with an upwind scheme. The speed used in the upwind scheme needs to be defined at the cell face, so we choose to reuse u * jk . For an active cell j and neighbor cells k ∈ N ( j) of the same fluid, we have the discrete form of the set of equations:

b n+1- j = b n+ j - ∆t |Ω j | ∑ k∈N ( j) |∂ Ω jk |u * jk b n+ jk + b n+ j ∆t |Ω j | ∑ k∈N ( j) |∂ Ω jk |u * jk , (5.22)
where

b n+ jk =    b n+ j if u * jk > 0, b n+ k otherwise.
(5.23)

Coupling for transport subsystem

The transport subsystem is discretized with conservative terms that are treated with an upwind scheme. In order to avoid mixed cells containing both fluid 1 and 2, we define a ghost state based on cells belonging to another fluid next to the interface and we use the ghost state to discretize conservative terms in the transport step.

For example, as presented in Fig. 5.2, cell j and cell k with real states of fluid 1 and of fluid 2 respectively are separated by an interface Γ. If the interface velocity u * jk is opposite to the direction of n jk , the upwind flux is then computed with the ghost state in cell k. Then (5.23) can be modified as:

b n+ jk =    b n+ j u * jk > 0, b n+
k,ghost otherwise.

(5.24)

We briefly recall how the ghost values are computed across the interface and refer the reader to [START_REF] Aslam | A partial differential equation approach to multidimensional extrapolation[END_REF] for a more detailed presentation. For the sake of simplicity we shall suppose here that the computational • compute an approximate value (b n ) j = (n • ∇b) j for (b n ) in the cell Ω j for j such that φ j < -∆x thanks to second order central difference formula,

• extrapolate the value of b n in the region φ > -∆x by solving the evolution equation

∂ t b n + H(φ + ∆x)n • ∇b n = 0, (5.25)
thanks to a WENO5 scheme as presented in [START_REF] Jiang | Weighted eno schemes for hamilton-jacobi equations[END_REF],

• propagate the value of b in the ghost region φ > 0 by solving the PDE

∂ t b + H(φ )(n • ∇b -b n ) = 0, (5.26)
with a method of lines.

Stability is ensured under the following classic CFL that only involves u * jk :

∆t transport |∂ Ω jk | |Ω j | max(|u j |, |u * jk |) ≤ 1 2 .
(5.27)

Let us note that in the low Mach regime, (5.27) is less restrictive than the acoustic CFL (5.20).

Note:

By combining the discretization of acoustic and transport subsystems, we recover a conservative form of the single-fluid Euler system.

ρ n+1- j = ρ n j - ∆t |Ω j | ∑ k∈N ( j) |∂ Ω jk |ρ n+ jk u * jk , (ρu) n+ j = (ρu) n j - ∆t |Ω j | ∑ k∈N ( j) |∂ Ω jk |π * jk n jk - ∆t |Ω j | ∑ k∈N ( j) |∂ Ω jk |(ρu) n+ jk u * jk , (ρE) n+ j = (ρE) n j - ∆t |Ω j | ∑ k∈N ( j) |∂ Ω jk |π * jk u * jk - ∆t |Ω j | ∑ k∈N ( j)
|∂ Ω jk |(ρE) n+ jk u * jk .

(5.28a)

(5.28b)

(5.28c)

This conservative form shows the necessity of updating order between the acoustic and transport subsystem as given in Section 3.3.

Approximation of the diffusion subsystem

In this section, we give the discretization of the diffusion subsystem and introduce the phase change model based on energy jump condition model [START_REF] Kharangate | Review of computational studies on boiling and condensation[END_REF]. The diffusion subsystem (3.22) can be discretized as:

t > 0, x ∈ D              ρ n+1 j = ρ n+1- j , (ρu) n+1 j = (ρu) n+1- j + ∆t{∇ • S} j , (ρE) n+1 j = (ρE) n+1- j + ∆t{∇ • (Su)} j + ∆t{∇ • (K ∇T )} j , φ n+1 j = φ n+1- j .
(5.29a)

(5.29b)

(5.29c)

(5.29d)

Terms related to the diffusion phenomena in (5.29) are approximated by the finite volume discretization with variables at step t n .

Discretization of terms related to viscous diffusion

The viscous subsystem is approximated by using a classical centered second-order finite-volume method. To eliminate errors related to the velocity jump across the interface, when the viscous stress tensor approximation involves the velocity of another fluid, this velocity should be corrected by the velocity jump given by the mass transfer in Section 5.2.3. The centered second-order discretization of viscous stress tensor of cells near the interface could involve a cell not close to the interface belonging to another fluid. As presented in Fig. 5.5, the viscous stress tensor discretization for cell j (colored cell) involves the velocity of cell k (white cell).

Similar to the interface velocity extrapolation, we employ a constant extrapolation to propagate the velocity jump to undefined cells. The time step constraint for the viscous step reads:

∆t viscosity ≤ |∂ Ω j | 2 ρ j 2|Ω jk | 2 µ j .
(5.30)

Discretization of term related to heat diffusion

The heat diffusion term can be discretized as:

{∇ • (K ∇T )} j = 1 |Ω jk | ∑ k∈N ( j)
|∂ Ω jk |n jk • {K ∇T } jk .

(5.31)

In the framework of sharp interface, the heat flux n jk • {K ∇T } jk can be rewritten as:

n jk • {K ∇T } jk =      K j + K k 2 T k -T j ∆x jk ∀k ∈ N + ( j), K j ∇T Γ j • n jk ∀k ∈ N -( j).
(5.32a)

(5.32b)

where

∇T Γ j = lim x→y x∈Ω j
∇T for y ∈ Γ(t n ), t n > 0 is the temperature gradient at the interface of Ω j . The estimation of heat flux across the interface could be divided into two cases: without and with phase change.

Without phase change

Without phase change, heat flux is continuous across the interface. For Ω j ⊂ D 1 and Ω k ⊂ D 2 as depicted in Fig. 5.3, the continuity of heat flux can be expressed as:

∀k ∈ N -( j), K j T j -T Γc ∆x Γ j = K k T Γc -T k ∆x Γ k , (5.33)
where T Γc is the temperature at the interface without the effect of phase change. With (5.33), we can fully determinate T Γc , that is

T Γc = K j T j ∆x Γ k + K k T k ∆x Γ j K j ∆x Γ k + K k ∆x Γ j , (5.34)
with

∆x Γ j = ∆x jk |φ j | |φ j -φ k | ∆x Γ k = ∆x jk |φ k | |φ j -φ k | . (5.35)
Without phase change, heat flux across the interface can be given as:

∀k ∈ N -( j), K j n jk • ∇T Γ j = K j T Γc -T j ∆x Γ j .
(5.36)

T x j k k ′ Γ T Γ ∆x Γ j ∆x Γ k

With phase change

In the present work, we are interested in low Mach regime. Here we simplify the thermodynamic model and consider a known and constant saturation temperature at the interface, that is a uniform Dirichlet condition on the interface temperature. More complex thermodynamic [START_REF] Faccanoni | Numerical simulation with finite volume of dynamic liquid-vapor phase transition[END_REF] could be always possible which will give varying saturation temperatures according to the fluid states close to the interface.

The most straight way to approximate the heat flux K n jk • ∇T Γ jk near the interface is that:

∀k ∈ N -( j), K j n jk • ∇T Γ j = K j T sat -T j ∆x Γ j .
(5.37)

Approximation given by (5.37) suffers from the problem of stability, In the precedent works [START_REF] Sato | A sharp-interface phase change model for a mass-conservative interface tracking method[END_REF][START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF][START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF], an implicit approach is realized to overcome the problem of stability. In this work, we want to propose an explicit and simple algorithm to estimate the heat flux at the interface with the presence of phase change.

For any cell j near the interface with a neighbor k ∈ N -( j) , to avoid the problem of stability related to the interface heat flux estimation, we use another neighbor cell in the direction n k j belonging to the same fluid to estimate the heat flux K j n jk • ∇T Γ j . To clearly illustrate the algorithm, we start from a one-dimensional problem. As depicted in Fig. 5.3, the cell k ′ belongs to the same fluid as cell j and in the direction of n k j , the interface heat flux in the direction of n jk can be estimated as:

K j n jk • ∇T Γ j = 1 2 K j (3 T sat -T k ′ ∆x k ′ j + ∆x Γ j - T j -T k ′ ∆x k ′ j
).

(5.38)

The discretization given by (5.38) does not suffer from extra constraint on stability compared with the bulk flow. The interface temperature is taken into account as a boundary condition, this approximation has a good ability of convergence when we refine the mesh.

Regarding the multi-dimension problem, here we show the principle of 2D cases, the extension to 3D will be straightforward. To be clear, we first define "explicit interface heat flux availability", which means that the interface heat flux of cell j can be discretized in the direction n jk for k ∈ N -( j) by using (5.38). This availability corresponds to the existence of a neighbor cell belonging to the same fluid in the opposite direction. With this availability, we can separate all possibilities of the interface geometry into 3 cases:

(Case1) Available in all directions

As depicted in Fig. 5.4a, the interface heat flux in all directions can be approximated by using (5.38). The interface heat flux in the normal direction can be approximated by:

K j ∇T Γ j • n = K j (∂ x T Γ j ∂ x φ + ∂ y T Γ j ∂ y φ ) (5.39)
where

∂ y T Γ j = lim x→y x∈Ω j ∂ y T and ∂ y T Γ j = lim x→y x∈Ω j ∂ y T for y ∈ Γ(t n ), t n > 0 are temperature gradient
for Ω j at the interface in the direction of x and y respectively.

(Case2) Available in at least one direction but not all directions

As shown in Fig. 5.4b, heat flux K j ∇T Γ j • n y in y direction can be estimated by using our explicit model, while the interface heat flux in x direction is not available by employing (5.38). In this case, the heat flux in the normal direction can be approximated by:

K j ∇T Γ j • n = K j ∂ y T Γ j /∂ y φ . (5.40)
Then the interface heat flux in x direction can be estimated by:

K j ∂ x T Γ j • n x = K j ∇T Γ j • n∂ x φ • n x .
(5.41)

(Case3) Not available in any direction

This case is generally related to the separation or coalescence of bubbles. For this problem, our model has some difficulties to give an explicit prediction of the interface temperature gradient with the presence of phase change.

For cells near the interface, when the interface heat flux is not available in all directions, its interface heat flux can be defined by propagating information in cells of (Case1) and (Case2) with a constant extrapolation.

The time-step accounting for the stability conditions on heat diffusivity follows the classical constraint:

∆t heat ≤ ρ j C p, j |∂ Ω jk | 2 2|Ω j | 2 K j .
(5.42)

Mass transfer estimation

For the approximate Riemann solver at the interface, the discretization of velocity jump which is related to the mass transfer rate is still unknown. The present work approximates the mass transfer by energy jump condition model: the jump of heat flux acts as the absorbed energy or emitted energy in the vaporisation or condensation procedure. According to (2.9), the approximation of mass transfer rate ṁ between two neighbor cells Ω j ⊂ D 1 and Ω k ⊂ D 2 :

ṁ jk = K j ∇T Γ j -K k ∇T Γ k L heat . (5.43)
Following (2.6), the velocity jump across the interface u • n Γ jk between two neighbor cells Ω j ⊂ D 1 and Ω k ⊂ D 2 can be approximated as:

u • n Γ jk = ṁ jk ( 1
ρ j - 1 ρ k ).
(5.44)

The mass conservation in the procedure of phase change can be described as:

ṁ jk = ρ j (u Γ jk -u j ) • n j = ρ k (u Γ jk -u k ) • n k . (5.45)
where u Γ jk = lim x∈y u for y ∈ Γ(t n ), t n > 0 is the velocity of the interface between Ω j and Ω k for k ∈

N -( j).
As shown in Fig. 5.5, the interface velocity is initialized only in cells near to the interface (colored cells). In the procedure of Level Set advection, a smooth velocity field could provide a more accurate description of the interface (Level Set advection) and decrease the influence of the redistancing step. It is better to use the interface velocity rather than the local fluid cell to update the Level Set function. A constant extrapolation is applied to propagate the interface velocity to undefined cells, for more details please refer to the Ghost cell definition in Section 5.2.2 and the work of Aslam [START_REF] Aslam | A partial differential equation approach to multidimensional extrapolation[END_REF].

Approximation of system with barotropic EOS

As mentioned in Chapter 2, there is always possible to close the system (2.1) with other EOS. Approximation of system (2.1) with barotropic EOS is given in Appendix. A.

Time Step Constraint

By combining (5.20) (5.21) (5.27) (5.42) and (5.30), the global time step restriction accounting for acoustic, advection, diffusion and capillary effects can be given as: ∆t = min(∆t acoustic , ∆t capillary , ∆t transport , ∆t heat , ∆t viscosity ).

(5.46)

Interface movement

As the interface is advected by the fluid velocity field, the value of φ evolves with time and may change its sign. For the cell that changes its sign, its real state is replaced by the ghost state in that cell.

Chapter 6

Low Mach correction

In the present work, we are interested in flows that involve two phases both represented by two compressible materials separated by an infinitely thin interface. The interface is considered sharp and treated as a contact discontinuity. The Ghost Fluid Method (GFM) [START_REF] Fedkiw | A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF][START_REF] Lalanne | On the computation of viscous terms for incompressible two-phase flows with Level Set/Ghost Fluid method[END_REF][START_REF] Caiden | A numerical method for two-phase flow consisting of separate compressible and incompressible regions[END_REF][START_REF] Liu | Ghost fluid method for strong shock impacting on material interface[END_REF][START_REF] Liu | The ghost fluid method for compressible gas-water simulation[END_REF][START_REF] Zou | An accurate sharp interface method for two-phase compressible flows at low-mach regime[END_REF] is a very popular family of methods that propose discretization strategies for simulation such flows. They rely on a Level Set function to describe the position of the interface [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations[END_REF][START_REF] Sethian | Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science[END_REF][START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF]. A common approach in the GFM consists in coupling single fluid numerical schemes like Godunov-type methods across the interface in order to accurately approximate the fluid motion in the vicinity of the interface. Unfortunately, in the case of single-fluid flows it is now well-established that standard Godunovtype methods may dramatically lose their accuracy when the material velocity is small compared to the sound velocity of the medium and when the computation domain is discretized over a quadrangle (resp. hexahedral) grid in 2D (resp. 3D). This low-velocity regime is referred to as the low Mach regime. This flaw has been extensively studied in the literature and several possible fixes are available (see for example [START_REF] Turkel | Preconditioned methods for solving the incompressible and low speed compressible equations[END_REF][START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit[END_REF][START_REF] Li | An all-speed Roe-type scheme and its asymptotic analysis of low Mach number behaviour[END_REF][START_REF] Murrone | Behavior of upwind scheme in the low mach number limit: Iii. preconditioned dissipation for a five equation two phase model[END_REF][START_REF] Dellacherie | Analysis of godunov type schemes applied to the compressible euler system at low mach number[END_REF][START_REF] Rieper | A low-Mach number fix for Roe's approximate Riemann solver[END_REF][START_REF] Dellacherie | Construction of modified Godunov-type schemes accurate at any Mach number for the compressible Euler system[END_REF][START_REF] Chalons | An all-regime lagrange-projection like scheme for the gas dynamics equations on unstructured meshes[END_REF]7,[START_REF] Xie | An accurate and robust HLLC-type Riemann solver for the compressible Euler system at various Mach numbers[END_REF]). This question has also been addressed concerning two-phase flows for diffuse interface models [START_REF] Chalons | An all-regime Lagrange-Projection like scheme for 2D homogeneous models for two-phase flows on unstructured meshes[END_REF][START_REF] Peluchon | A robust implicit-explicit acoustic-transport splitting scheme for two-phase flows[END_REF] where the interface is considered as a numerically diffused zone which can reduce the density and sound speed ratio. Comparing with the sharp interface method, the challenge associated with large discontinuities across the interface could be reduced as these discontinuities are smoothed by the mixing layer.

In the case that we consider here, the properties of the medium may experience a severe jump across the interface. Therefore one cannot assume that the orders of magnitude associated with the fluid parameters are the same for both phases. Consequently, it is necessary to define a Mach regime that pertains to each fluid separately. In this context, the classic low Mach fixes are no longer sufficient to propose an accurate discretization of the flow equations when both Mach numbers of each fluid are low.

In the present work, we propose a new low Mach regime analysis that is adapted to compressible flows with sharp interfaces including surface tension effects. Moreover, we also propose a new low Mach correction of the GFM presented in [START_REF] Zou | An accurate sharp interface method for two-phase compressible flows at low-mach regime[END_REF] that accounts for the jump of characteristic magnitude with respect to density and sound velocity that provides good accuracy in the low Mach regime.

The chapter is structured as follows: With the global system given in Chapter 2, we define the low Mach regime in Section 6.1. The low Mach behavior of numerical model developed in Chapter 5 far from the interface is studied in Section 6.2 suggesting that the low Mach fixes proposed in [START_REF] Peluchon | A robust implicit-explicit acoustic-transport splitting scheme for two-phase flows[END_REF][START_REF] Chalons | An all-regime Lagrange-Projection like scheme for 2D homogeneous models for two-phase flows on unstructured meshes[END_REF] will enhance the behavior of the scheme far from the interface. In Section. 6.3 we investigate these low Mach schemes in the vicinity of the interface and we expose truncation errors that are driven by the jump of magnitude with respect to both the density and the sound velocity of the material that dramatically affects the accuracy of the low Mach schemes. A new low Mach scheme is then proposed in Section 6.4 to address this issue. This new numerical method is tested in Section. 6.5 against several 2D tests that show significant improvement in the simulation results over the previous low Mach fixes.

Definition of low Mach regime

In this chapter, for the sake of simplicity, we neglect the source term related to the external force f and heat transfer. If n is the unit normal vector to the interface oriented from fluid 1 to fluid 2, we suppose that the motion equations for the fluid 1 and the fluid 2 are coupled across the interface thanks to the following two hypotheses: By following the work of [START_REF] Majda | Compressible fluid flow and systems of conservation laws in several space variables[END_REF][START_REF] Klein | Semi-implicit extension of a godunov-type scheme based on low mach number asymptotics i: One-dimensional flow[END_REF][START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit[END_REF], the behavior of this system of equation will be investigated with respect to the Mach regime at continuous level. We will perform an asymptotic analysis in the regime M i ≪ 1 that will provide a set of conditions that must be fulfilled by the flow parameters during the motion and also at the initial conditions [START_REF] Rieper | A low-Mach number fix for Roe's approximate Riemann solver[END_REF][START_REF] Cordier | An asymptotic-preserving all-speed scheme for the euler and navier-stokes equations[END_REF][START_REF] Haack | An all-speed asymptotic-preserving method for the isentropic euler and navier-stokes equations[END_REF].

Definition of non-dimensional variables

In order to study the behavior of numerical schemes in the low Mach regime, we first need to characterize the flow with respect to the Mach number. To that end, we need characteristic values for length and time that will be respectively noted x, t, û = x/ t. This allows to define the following non-dimensional quantities and nabla operator that are common to both fluids

t = t/ t, x = x/ x, ∇ = ∂ ∂ x1 , . . . , ∂ ∂ xd T = x ∇, ũ = u/ û, ũ Γ = u Γ / û.
The non-dimensional domain occupied by the fluid i and the interface will be respectively noted Di and Γ. We now define a set of non-dimensional variables for each fluid i = 1, 2 by setting

⌊ ρ⌋ i = ρ/ ρi , ⌊ p⌋ i = p/ pi , ⌊ c⌋ i = c/ ĉi , ⌊ ẽ⌋ i = e/ êi , ⌊ p Γ ⌋ i = p Γ / pi , S i = 1 Re i S/( ρi û2 ), (6.2)
where ρi , ĉi , pi = ρi ĉ2 i , ûi , êi = pi / ρi , and Re i denote respectively a characteristic density, sound velocity, pressure, fluid velocity, internal energy and Reynolds number associated with the fluid i = 1, 2. In the following, we will enclose an expression between the signs ⌊•⌋ i in order to indicate that it has to be considered as a non-dimensional expression with respect to the characteristic quantities of the fluid i. For example, ⌊ ρ⌋ 1 is the dimensionless density variable with respect to the characteristic density ρ1 of fluid 1. We note M i the Mach number associated with the fluid i that is defined by

M i = û/ ĉi .
The dimensionless total energy is then be defined as:

Ẽ i = E/ êi = ⌊ ẽ⌋ i + M 2 i . (6.3)
Furthermore, we introduce the following non dimensional parameters α i ∈ R + and β i ∈ R + that characterize the magnitude of characteristic densities, sound velocities and fluid velocities with respect to the fluid 1

α i = ρi / ρ1 , β i = ĉi / ĉ1 .
And finally, we also note

α m = min(α 1 , α 2 )/ max(α 1 , α 2 ) ≤ 1, β m = min(β 1 , β 2 )/ max(β 1 , β 2 ) ≤ 1.
(6.4)

Dimensionless evolution equation and low Mach regime

For the sake of readability, when there is no ambiguity, we shall replace ⌊•⌋ i by • to refer to a non-dimensional expression with respect to the fluid i. With this simplified notation, the governing equation (2.1) for the bulk flow in D i can be written in a dimensionless form:

t > 0, x ∈ Di (t)                ∂ t ρ + ∇ • ( ρ ũ) = 0, ∂ t ( ρ ũ) + ∇ • ( ρ ũ ⊗ ũ) + 1 M 2 i ∇ p = 1 Re i ∇ • S, ∂ t ( ρ Ẽ) + ∇ • (( ρ Ẽ + p) ũ) = M 2 i Re i ∇ • ( S ũ).
(6.5a) (6.5b) (6.5c)

The non-dimensional form of the evolution (2.4) is

∂ t φ + ũΓ • ∇φ = 0, t > 0, x ∈ D. (6.6)
Concerning the boundary conditions across the interface Γ: (H1) implies that (u + u Γ jk ) • n is continuous across Γ and the relation (6.1) yields a non-dimensional jump relation that reads

    lim x→ỹ x∈ D1 ũ -lim x→ỹ x∈ D2 ũ    1 =⌊ ũ Γ ⌋ 1 (ỹ, t), for ỹ ∈ Γ( t).     lim x→ỹ x∈ D1 p -lim x→ỹ x∈ D2 p    1 =⌊ p Γ ⌋ 1 (ỹ, t), for ỹ ∈ Γ( t). (6.7)
Before going any further, let us underline that the overall two-phase system formed by (6.5), (6.6) and the boundary conditions on Γ depends on both M 1 and M 2 .

In this context, we will now characterize what is intended in this paper as low Mach regime. We suppose that M 1 ≪ 1 and M 2 ≪ 1, and for a dimensionless variable Ãi of fluid i, we consider an asymptotic expansion with respect to M i as follows:

à i (x,t) = M 0 i à (0) i (x, t) + M 1 i à (1) i (x, t) + M 2 i à (2) i (x, t) + • • • , x ∈ Di ( t), t > 0. ( 6.8)
We now examine the behavior of à (k) i in the limit regime M i → 0: introducing expression (6.8) in (6.5) and collecting the terms with equal power of M i , we obtain:

• Order M -2 i : ∇ p(0) = 0, x ∈ Di ( t), t > 0, (6.9)
• Order M -1 i : ∇ p(1) = 0, x ∈ Di ( t), t > 0, (6.10)

• Order M 0 i : t > 0, x ∈ Di (t)            ∂ t ρ(0) + ∇ • ( ρ(0) ũ(0) ) = 0, ∂ t ( ρ(0) ũ(0) ) + ∇ • ( ρ(0) ũ(0) ⊗ ũ(0) ) + 1 M 2 i ∇ p(2) = 1 Re i ∇ • S(0) , ∂ t ( ρ(0) Ẽ(0) ) + ∇ • (( ρ(0) Ẽ(0) + p(0) ) ũ(0) ) = 0. (6.11a) (6.11b) (6.11c)
where p(2) is the dynamic pressure related to the fluid motion, and p(0) is the thermodynamic pressure. The thermodynamic pressure p(0) can be expressed as a function of ρ(0) and Ẽ(0) thanks to the EOS

(2.3): p(0) = p EOS i ( ρ(0) , Ẽ(0) ),
in Di ( t), t > 0. (6.12)

Thanks to (6.3), we can see that Ẽ(0) i = ẽ(0) i and then using (6.12) and (6.11c) we obtain

∂ t p(0) + ũ • ∇ p(0) + ρ c2 ∇ • ũ(0) = 0, in Di ( t), for t > 0. (6.13)
We can now further characterize the limit regime M i → 0: equations (6.9) and (6.10) that pertain respectively to terms of order M -2 i and M -1 i imply that the pressure is homogeneous with respect to the space variable up to a fluctuation of magnitude of M 2 i :

∇ p = O(M 2 i ), x ∈ Di ( t), t > 0. (6.14)
We now add two supplementary hypotheses. First, we assume that the boundary conditions are chosen such that the space-homogeneous p(0) is also independent of t. Consequently, (6.13) implies that ∇ • ũ(0) = 0 in Di ( t), for t > 0. (6.15)

Let us emphasize that if such hypothesis is classic in the case of single material flow, it is much stronger in our case due to the presence of the two-phase boundary. Second, we suppose that the above property is also true across the interface, that is to say:

∇ • ( ũ(0) + ũ(0) Γ ) 1,2 = 0, on Γ. ( 6
         ∇ • ũ = O(M i ), t > 0, x ∈ D i (t), ∇ p = O(M 2 i ), t > 0, x ∈ D i (t), ∇ • ( ũ + ũ Γ ) 1,2 = O(M i ), t > 0, x ∈ Γ(t).
(6.17a) (6.17b) (6.17c)

An accurate scheme for low Mach two-phase flows should be able to preserve the defined low Mach regime. In the following sections, we will recall the sharp interface numerical scheme for compressible two-phase flows presented in [START_REF] Zou | An accurate sharp interface method for two-phase compressible flows at low-mach regime[END_REF], and examine its behavior for initial data in the low Mach regime defined by (6.17).

Low Mach behavior of the numerical scheme

In this section we will present a set of conditions that the discretization strategy of Chapter 5 needs to verify in order to accurately approximate the two-phase flow presented in Chapter 2 when the flow matches the low Mach definition (6.17). For the sake of brevity we will not examine the viscous stage of the numerical method as it did not have a significant impact on our numerical results in practice.

We first recast the first two subsystem we considered in Section 3.3 using non-dimensional form in the each domain Di . By omitting the external sources, the acoustic system reads

t > 0, x ∈ Di                  ρ∂ t (1/ ρ) -∇ • ũ = 0, ρ∂ t ũ + 1 M 2 i ∇ p = 0, ρ∂ t Ẽ + ∇ • ( p ũ) = 0, ∂ t φ = 0. (6.18a) (6.18b) (6.18c) (6.18d)
As the Mach number M i is explicitly involved in (6.18), it is also interesting to consider the evolution of ρ(0) , ũ(0) and Ẽ(0) through the acoustic phase. We get

t > 0, x ∈ Di              ∂ t ρ(0) = 0, ρ(0) ∂ t ũ(0) + ∇ p(2) = 0, ∂ t Ẽ(0) = 0, ∂ t φ = 0. (6.19a) (6.19b) (6.19c) (6.19d)
The transport subsystem in non-dimensional form becomes

t > 0, x ∈ Di ∂ t b + ∇ • ( b ũ) -b( ∇ • ũ) = 0, b = ( ρ, ρ ũ, ρ ṽ, ρ w, φ ) T . (6.20)
In order to examine the behavior of the numerical scheme, we will follow the lines of [START_REF] Chalons | An all-regime lagrange-projection like scheme for the gas dynamics equations on unstructured meshes[END_REF] by considering smooth solutions of (6.5) that fulfill the low Mach regime requirements (6.17). We note (x, t) → b any rescaled fluid parameter associated with this solution. We commit an abuse of notation by setting b(x j , t) = b j . Then we shall inject such solutions into the rescaled numerical scheme and study the resulting truncation error.

Low Mach behavior of the acoustic step in the bulk

We start by examining the acoustic step by considering a cell Ω j ⊂ D i . Let us note a i = max

{r | Ω r ⊂D i } (ρ r c r ),
for the sake of simplicity we suppose that the artificial acoustic impedance are defined using (5.9) with the uniform choice a j = a i for all cell Ω j ⊂ D i . Then, the non-dimensional form of (5.15) is

                         ρn j 1 ρn+ j - 1 ρn j -∆ t ∑ k∈N ( j) |∂ Ω jk | |Ω j | ũ * jk = 0, ρn j ( ũn+ j -ũn j ) + ∆ t M 2 i ∑ k∈N ( j) |∂ Ω jk | |Ω j | π * jk n jk = 0, ρn j ( Ẽn+ j -Ẽn j ) + ∆ t ∑ k∈N ( j) |∂ Ω jk | |Ω j | π * jk ũ * jk = 0, (6.21a) (6.21b) (6.21c)
where

       ũ * jk = n jk • ( ũ j + ũk ) 2 + 1 M i π j -πk 2 ãi , π * jk = π j + πk 2 -M i ãi ( ũk -ũ j ) 2 n jk . (6.22a) (6.22b)
The low Mach hypothesis (6.17) implies that that there exists two smooth functions A i and

B i of magnitude O(M 0 i ) such that we have t > 0 ũ(x k , t) • n = ũ(x j , t) • n + A i (x j , xk , t, n)|x k -x j |, p( xk , t) = p(x j , t) + M 2 i B i (x j , xk , t)|x k -x j | + O(M 2 i ∆ x2 ). (6.23a) (6.23b)
for any xk and x j in Di . Relations (6.22) and (6.23) yield that

       ũ * jk = ũ(x jk , t) • n jk - 1 2 ãi M i B i (x j , xk , t)|x k -x j | + O(M i ∆ x) + O(∆ x2 ), π * jk = p(x jk , t) - ãi 2 M i A i (x j , xk , t, n)|x k -x j | + O(M 2 i ∆ x) + O(M i ∆ x2 ). (6.24a) (6.24b)
Classic computations yield

∑ k∈ν( j) ũ(0) (x jk , t) • n jk = ∇ • ũ(0) (x j , t)∆ x + O(∆ x2 ), ∑ k∈ν( j) p(x jk , t) • n jk = ∇ p( x j , t)∆ x + O(M 2 i ∆ x2 ), (6.25a) (6.25b)
so that by injecting (6.24) and (6.25) into (6.21), we get

t > 0, x ∈ Di            ∂ t ρ(0) = O(∆ xM 0 i ) + O(∆ t), ∂ t ũ(0) + 1 M 2 i ∇ p = O(∆ xM 0 i ) + O ∆ x M i + O(∆ t), ∂ t Ẽ(0) = O(∆ xM 0 i ) + O(∆ t). (6.26a) (6.26b) (6.26c)

Low Mach behavior of the transport step in the bulk

The transport step in non-dimensional form for a cell Ω j ⊂ D i reads

Ω j ∩ Di ̸ = / 0, bn+1- j = bn+ j - ∆t | Ω j | ∑ k∈N j |∂ Ω jk | ũ * jk bn+ jk + bn+ j ∆ t | Ω j | ∑ k∈N j |∂ Ω jk | ũ * jk , (6.27)
where the flux terms bn+ jk are non-dimensional form of (5.23) and (5.24) with respect to the fluid i. By injecting injecting the low Mach solution into (6.27), standard computations show that it yields

t > 0, x ∈ Di , ∂ t b + ∇ • ( ũb ) -b( ∇ • ũ) = O(M i 0 ∆ x) + O(∆ t) + O(M i ∆ x), b ∈ { ρ, ρ ũ, ρ ṽ, ρ w, φ }. (6.28)

Low Mach regime accuracy analysis and fix for the bulk flow

Results presented in sections 6.2.1 and 6.2.2 are coherent with the literature: in the acoustic step, a truncation error term of magnitude O(∆ x/M i ) that appears in (6.26b) suggests that important loss of accuracy may occur when M i ≪ 1, while in the transport step, the truncation error remains uniform with respect to M i . Many works over the past years have proposed low Mach corrections for single-fluid flows [START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit[END_REF][START_REF] Li | An all-speed Roe-type scheme and its asymptotic analysis of low Mach number behaviour[END_REF][START_REF] Dellacherie | Construction of modified Godunov-type schemes accurate at any Mach number for the compressible Euler system[END_REF][START_REF] Bruel | A low Mach correction able to deal with low Mach acoustics[END_REF][START_REF] Chalons | An all-regime lagrange-projection like scheme for the gas dynamics equations on unstructured meshes[END_REF]. The fix usually amounts to use a centered pressure evaluation at the cell interfaces in the low Mach regime. In the case of a two-phase flow, let us recall two modifications that have been proposed in the literature.

Noting (M i ) jk a local evaluation of the Mach number M i at the interface ∂ Ω jk and setting θ jk = min((M 1 ) jk , (M 2 ) jk , 1), Peluchon et al. [START_REF] Peluchon | A robust implicit-explicit acoustic-transport splitting scheme for two-phase flows[END_REF] proposed to replace π * by altering terms related to velocity jump as follows:

π * ,θ ,AW jk = a k π j + a j (π k + p Γ jk ) a j + a k -θ jk a j a k a j + a k n jk (u k -u j + u Γ jk ). (6.29)
Chalons et al. [START_REF] Chalons | An all-regime Lagrange-Projection like scheme for 2D homogeneous models for two-phase flows on unstructured meshes[END_REF] studied a slightly different low Mach correction with a centered pressure which can be expressed as:

π * ,θ ,CP jk = (1 -θ jk ) π j + π k + p Γ jk 2 + θ jk π * jk . (6.30)
Let us underline that the two-phase models studied in both Peluchon et al. [START_REF] Peluchon | A robust implicit-explicit acoustic-transport splitting scheme for two-phase flows[END_REF] and Chalons et al. [START_REF] Chalons | An all-regime Lagrange-Projection like scheme for 2D homogeneous models for two-phase flows on unstructured meshes[END_REF] are different from our framework as they involve potential mixture regions where both fluids can simultaneously be present. Applying either (6.29) or (6.30) in our case succeeds in controlling the error term in the momentum update of the acoustic step in pure fluid i region. Indeed, either corrections will enable new truncation error estimate

t > 0, x ∈ Di ρ∂ t ũ + 1 M 2 i ∇ p = O(∆ xM 0 i ) + O θ ∆ x M i + O(∆ t), (6.31)
where θ is an upper bound for all θ jk .

Asymptotic behavior of existing low Mach schemes across the interface

We have seen in Section 6.2 that in the low Mach regime, our discretization strategy may suffer a loss of accuracy in each pure fluid region D i of the computational domain. This drawback is a classic pathology in the case of single fluid flows that can be improved by modifying the pressure discretization. We shall now investigate in the present section the behavior of the numerical scheme across the two-phase interface in the low Mach regime. We will see that different variables in the low Mach regime are triggered by the abrupt jump of the medium properties across the interface. Once again, we consider a smooth solution of (6.5) that fulfill the low Mach regime requirements (6.17). We consider a cell Ω j such that φ n j > 0 that lies in the vicinity of Γ so that N -( j) ̸ = / 0. Let V be a neighborhood of x j so that V ∩ D 1 ̸ = / 0 and V ∩ D 2 ̸ = / 0. Thanks to the smoothness of the velocity we know that ⌊ ũ⌋ i and its derivatives are of magnitude O(M 0 i ) within the entire region V for i = 1, 2. Thus there exist a smooth function A such that the variation of the normal velocity at the boundary Ω j reads

                   ũ(x j , t) • n jk 1 = ũ(x jk , t) • n jk + A(x jk , x j ,t, n k j ) |∆x jk | 2 + O(∆ x2 ) 1 for k ∈ N ( j), ũ(x k , t) • n jk 1 = ũ(x k j , t) • n jk + A(x k j , xk ,t, n jk ) |∆x k j | 2 + O(∆ x2 ) 1 for k ∈ N + ( j), ũ(x k , t) • n jk 2 = ũ(x k j , t) • n jk + A(x k j , xk ,t, n k j ) |∆x k j | 2 + O(∆ x2 ) 2 for k ∈ N -( j). (6.32a) (6.32b) (6.32c)
Similarly according to (6.17b), ∇ p i is of magnitude O(M 2 i ) within V ∩ D i so that there exist two smooth functions B 1 and B 2 respectively of magnitude O(M 0 1 ) and O(M 0 2 ) such that

                   p(x j , t) 1 = p(x jk , t) + M 2 1 B 1 (x jk , x j , t) |∆x jk | 2 + O(M 2 1 ∆ x2 ) 1 for k ∈ N ( j), ⌊ p(x k , t)⌋ 1 = p(x k j , t) + M 2 1 B 1 (x jk , xk , t) |∆x k j | 2 + O(M 2 1 ∆ x2 )) 1 for k ∈ N + ( j), ⌊ p(x k , t)⌋ 2 = p(x k j , t) + M 2 2 B 2 (x jk , xk , t) |∆x k j | 2 + O(M 2 2 ∆ x2 )) 2 for k ∈ N -( j). (6.33a) (6.33b) (6.33c)
The jump condition across the interface can be expressed as:

p(x jk , t) -p(x k j , t) 1 = p Γ jk 1 when k ∈ N -( j), (6.34) and ũ
(x jk , t) -ũ(x k j , t) 1 = ũ Γ jk 1 when k ∈ N -( j). (6.35)
It is important to note that ⌊ p(x k , t)⌋ 1 (resp. ⌊ ãk ⌋ 1 ) when k ∈ N -( j) is not of magnitude O(1) as p(x k ) (resp. a k ) is evaluated in D 2 but it is rescaled using the characteristic values of the fluid 1.

Nevertheless when k ∈ N -( j), one can exhibit express ⌊p(x k )⌋ 1 and ⌊ ãk ⌋ 1 using terms of magnitude O(1) as follows:

⌊ p(x k , t)⌋ 2 = ⌊ p(x k , t)⌋ 1 /(α 2 β 2 2 ), (6.36)
and

⌊ ãk ⌋ 1 = ⌊ ãk ⌋ 2 (α 2 β 2 ). (6.37)
Thus (6.33c) can now be expressed using characteristic values of the fluid 1, that is to say

⌊ p(x k , t)⌋ 2 = ⌊ p(x k , t)⌋ 1 (α 2 β 2 2 ) = p(x k j , t) + M 2 1 B 2 (x jk , xk , t) |∆x k j | 2 α 2 + O(M 2 1 ∆ x2 )α 2 1 . (6.38)
Following the same standard lines as in Section 6.2 for the bulk flow, we can now reinject the low Mach solution into the expression of the numerical scheme in order to express the truncation error. We examine the cell Ω j and consider that it is mostly occupied by the fluid 1 as we supposed that φ j > 0. We will study the behavior of two elements of the numerical scheme: the discretized pressure at ∂ Ω jk and the update of ρ(0) j when k ∈ N -( j) ̸ = / 0.

Low

= π j + πk + p Γ jk 2 1 + θ jk π * jk - π j + πk + p Γ jk 2 1 . (6.39)
By (6.32), (6.33), (6.34), (6.37), (6.38) and (6.30), tedious calculations show that π * ,θ,CP

jk 1 = p(x jk , t) 1 + O(M 2 1 ∆ x2 ) + θ jk M 1 ã j 1 ⌊ ãk ⌋ 2 α 2 β 2 ã j 1 + ⌊ ãk ⌋ 2 α 2 β 2 A(x j , xk , t, n jk )|x k -x j | -θ jk ⌊ ãk ⌋ 2 α 2 β 2 ã j 1 + ⌊ ãk ⌋ 2 α 2 β 2 B 1 (x j , x jk , t)|x jk -x j | M 2 1 +θ jk ã j 1 α 2 ã j 1 + ⌊ ãk ⌋ 2 α 2 β 2 B 2 (x jk , xk , t)|x k -x jk | M 2 1 + (1 -θ jk ) 2 B 1 (x j , x jk , t)|x jk -x j | -B 2 (x k , x jk , t)α 2 |x k -x jk | M 2 1 . (6.40)
Then we can see that in the fourth line of (6.40), when α 2 becomes large, it may produce important errors for the estimation of the pressure gradient in Ω j . Indeed, let us consider π * ,θ,CP

jk -π j 1 ∆ x = p( x jk ) -p( x j ) 1 ∆ x + O(α 2 M 2 1 ) + O(M 2 1 ) + O(θ jk M 1 ) + O(θ jk α 2 M 2 1 )
.

Using (6.33) we get π * ,θ,CP jk -π j 1 ∆ x = O(α 2 M 2 1 ) + O(M 2 1 ) + O(θ jk M 1 ) + O(θ jk α 2 M 2 1 ). (6.41)
Therefore, when α 2 = O(1/M 1 ), this suggests that the discrete gradient estimation in the vicinity of the interface may grow to reach a O(1/M 1 ) magnitude which violates the low Mach hypothesis (6.17b). Analysis given in Appendix. B shows that the AW type correction could involve a truncation error of magnitude O(α 2 ) to the discrete pressure gradient in the momentum update, which also disobeys he low Mach hypothesis (6.17b).

Low Mach flow across the interface: evolution of the density near the interface during the acoustic step

We now evaluate the discretized normal velocity ũ * jk 1

for k ∈ N -( j): by combining (6.32), (6.33), (6.34), (6.37), (6.38) and (5.18a) we obtain

ũ * jk 1 = ũ(x jk , t) 1 • n jk + ⌊ ãk ⌋ 2 α 2 β 2 ã j 1 + ⌊ ãk ⌋ 2 α 2 β 2 - 1 2 A(x j , xk , t, n jk )|∆x k j | - 1 ã j 1 + ⌊ ãk ⌋ 2 β 2 α 2 B 1 (x j , x jk , t) |∆x k j | 2 + B 2 (x k , x jk , t)α 2 |∆x k j | 2 M 1 +O(∆ xM 1 ) + O(∆ x2 ). (6.42)
We can now inject (6.42) and (6.24a) into (6.21a) and considering only M 0 1 terms we get

∂ t ρ(0) 1 = -∑ k∈N -( j) ⌊ ãk ⌋ 2 α 2 β 2 ã j 1 + ⌊ ãk ⌋ 2 α 2 β 2 - 1 2 A(x j , xk , t, n jk ) + O(∆ xM 0 1 ) + O(∆ t). (6.43)
Then one can see that due to the discrepancy of characteristic quantities that occurs when α 2 ≫ 1 and

β 2 ≫ 1, relation (6.43)
suggests that ρ(0) can no longer remain constant up to a terme of magnitude O(1).

Truncation error for the transport update in the vicinity of the interface

We now turn to the transport step and by using (6.42) we first remark the normal velocity estimation verifies ũ *

jk 1 = ũ(x jk , t) 1 • n jk + O(∆ xM 0 1 α 0 m β 0 m ) + O(M i ) 1 . (6.44)
Now, for b ∈ {ρ, ρu, ρv, ρw, ρE} and k ∈ N -( j), following (5.23) we know that b n+ jk is be obtained either by the upwind choice within the same fluid either by a linear Ghost fluid extrapolation. This means that in either case

bn+ jk = b(x jk , t) + O(∆ xM 0 1 α 0 m β 0 m ). (6.45)
Therefore, in the neighbourhood of Γ, we obtain a similar truncation error as in the bulk fluid case (6.28), that is to say if b ∈ { ρ, ρ ũ, ρ ṽ, ρ w, φ }, we have

∂ t b + ∇ • ( ũb ) -b( ∇ • ũ) = O(M i 0 ∆ xα 0 m β 0 m ) + O(∆ t) + O(M i ∆ xα 0 m β 0 m ). (6.46)
Consequently, (6.46) suggests that the transport step does not perturb the low Mach regime.

New low Mach correction

We have seen in Section 6.3 that within cells neighboring the interface Γ, even the classic low Mach corrected schemes failed to preserve the low Mach regime defined by (6.17). The cause of the error seems to originate from the jump of the characteristic values that are associated with the fluid across Γ. Indeed, if one considers the expression of the dimensional pressure in D i , we get from the asymptotic expansion with respect to M i that

p = pi p = p(0) ρi ĉ2 i + p(1) ρi ĉi û + p(2) ρi û2 + • • • in D i . (6.47)
In the low Mach regime defined by (6.17) we know that the first two terms in (6.47) are constant in D i so that the variations of p are dominated by the dynamic pressure with a magnitude ρi û2 . Across the interface Γ, the dynamic pressure with a magnitude undergoes a jump that is characterized by

ρ2 û2 ρ1 û2 = α 2 . (6.48)
This corroborates the fact that the error term of magnitude O 1 α m in (6.41) is issued from the discretized centered pressure term. For similar reasons, relation (6.43) shows that the the presence of α 2 and β 2 in (6.42) does not allow to recover an accurate estimate of ∇ • ( ũ(0) ) = 0 in the acoustic step. This results in variations of ρ(0) of magnitude O(1) that generate spurious oscillations in the thermodynamic pressure. In order to cure these problems, we propose to use a new numerical discretization for the pressure and the normal velocity terms at the cell interfaces. Our new numerical scheme should account for the variation of characteristic values associated with the dynamic pressure and also provide a reliable discretization of ∇ • ( ũ(0) ) in the vicinity of Γ. We introduce the following discretization of the pressure and the normal velocity at the interface ∂ Ω jk of a cell

Ω j for k ∈ N ( j) u * ,θ ,N jk = (1 -θ jk )n jk • u j + u k + u Γ jk 2 + θ jk n jk • a j u j + a k u k + u Γ jk a j + a k + π j -π k -p Γ jk a j + a k , (6.49a)
π * ,θ ,N jk = (1 -θ jk ) ρ k π j + ρ j (π k + p Γ jk ) ρ j + ρ k + θ jk π * jk . (6.49b)
We can now apply the lines presented in Section 6.3. Supposing that the flow verifies the low Mach hypotheses (6.17) and considering a cell

Ω j that is crossed by Γ, if k ∈ N -( j) we have that ũ * ,θ,N jk 1 = ũ(x jk , t) 1 • n jk + θ jk ⌊ ãk ⌋ 2 α 2 β 2 ã j 1 + ⌊ ãk ⌋ 2 α 2 β 2 - 1 2 A(x j , xk , t, n jk )|∆x k j | - 1 ã j 1 + ⌊ ãk ⌋ 2 β 2 α 2 B 1 (x j , x jk , t) |∆x k j | 2 + B 2 (x k , x jk , t)α 2 |∆x k j | 2 M 1 + O(∆ xM 1 ) + O(∆ x2 ), (6.50)
and π * ,θ,N

jk 1 = p(x jk , t) 1 + O(M 2 ∆ x2 ) + θ jk M 1 ã j 1 ⌊ ãk ⌋ 2 α 2 β 2 ã j 1 + ⌊ ãk ⌋ 2 α 2 β 2 A(x j , xk , t, n jk )|∆x k j | - α 2 (1 -θ jk ) 1 + α 2 B 1 (x j , x jk , t) |∆x k j | 2 -B 2 (x k , x jk , t) |∆x k j | 2 M 2 1 -θ jk ⌊ ãk ⌋ 2 α 2 β 2 ã j 1 + ⌊ ãk ⌋ 2 α 2 β 2 B 1 (x j , x jk , t) |∆x jk | 2 M 2 1 + θ jk ã j 1 α 2 ã j 1 + ⌊ ãk ⌋ 2 α 2 β 2 B 2 (x jk , xk , t) |∆x k j | 2 M 2 1 .
(6.51) (6.51) gives a discrete pressure gradient:

π * ,θ,N jk -π j 1 ∆ x = p( x jk ) -p( x j ) 1 ∆ x + O(M 2 1 ) + O(θ jk M 1 ) + O(θ jk α 2 M 2 1 ). (6.52)
The truncation error uniform with respect to the density ratio, sound speed ratio and Mach number.

Injecting the estimates (6.50) into (6.21) yields the following truncation error to the density evolution for the acoustic step:

∂ t ρ(0) j 1 = O(∆xM 0 1 α 0 m β 0 m ). (6.53)
As a conclusion, (6.52) and ( 6.53) provides a truncation error which is uniform respect to M 1 , α m and β m to the discrete gradient estimation and density evolution, the low Mach hypothesis (6.17

α m ).
As for the transport step, the dimensionless corrected intermediate velocity ũθ,N jk (6.51) can be written into the same form as (6.44). Therefore, we can get the same conclusion about the transport step: the discrete transport step with the new low Mach correction is capable of satisfying the necessary condition of AP property.

An important feature of the new low Mach corrected solver for the acoustic step is that it can be associated with an approximate Riemann solver

(ζ ; W l , W r ) → W N RP (ζ ; W l , W r ), defined by W N RP (ζ ; W l , W r ) =                W l , if ζ < -a l , W * ,N,θ l , if -a l < ζ < 0, W * ,N,θ r , if 0 < ζ < a r , W r , if a r < ζ , (6.54)
with intermediate states W * ,N,θ l and W N,θ , * r defined by

u * ,θ ,N lr = (1 -θ lr ) u l + u r + u Γ lr 2 + θ lr a l u l + a r u r + a r u Γ lr a r + a r + π # l -π # r -p Γ lr a l + a r , (6.55a)
u * ,θ ,N rl = (1 -θ lr ) u l + u r -u Γ lr 2 + θ lr a l u l + a r u r -a l u Γ lr a r + a r + π # l -π # r -p Γ lr a l + a r , (6.55b)
π * ,θ ,N l = (1 -θ lr ) ρ r π l + ρ l (π r + p Γ lr ) ρ l + ρ r + θ lr π * lr , (6.55c)
π * ,θ ,N r = (1 -θ lr ) ρ r (π l -p Γ lr ) + ρ l π r ρ l + ρ r + θ lr π * rl , (6.55d)
1 ρ * ,θ ,N l = 1 ρ l + a l +a r 2 (u r -u l + u Γ jk ) + π l -π r -p Γ jk a l (a l + a r ) , (6.55e)
1 ρ * ,θ ,N r = 1 ρ r + a l +a r 2 (u r -u l + u Γ jk ) + π r -π l + p Γ jk a r (a l + a r ) , (6.55f)
v * ,θ ,N l = v l , v * ,θ ,N r = v r , (6.55g)
w * ,θ ,N l = w l , w * ,θ ,N r = w r , (6.55h)
E * ,θ ,N l = E l - π * ,θ ,N l u * ,θ ,N -π l u l a l , (6.55i)
E * ,θ ,N r = E r + π * ,θ ,N r u * ,θ ,N -π r u r a r . (6.55j)
With analysis in C, the following constraints for the acoustic impedance a l and a r :

if π r -π l + p Γ jk ≥ 0,        a l = ρ l c l + γ l +1 2 ρ l 2 π r -π l + p Γ jk ρ r c r + ρ l c l + u l -u r -u Γ jk + , a r = ρ r c r + γ r +1 2 ρ r 2 π l -π r -p Γ jk a l + ρ r c r + u l -u r -u Γ jk + , (6.56)
if π r -π l + p Γ jk ≤ 0,        a r = ρ r c r + γ r +1 2 ρ r 2 π l -π r -p Γ jk c l ρ l + c r ρ r + u l -u r -u Γ jk + , a l = ρ l c l + γ l +1 2 ρ l 2 π r -π l + p Γ jk a r + ρ l c l + u l -u r -u Γ jk + , (6
.57) can satisfy the positivity of ρ * ,θ ,N l and ρ * ,θ ,N r . In low Mach regime, the second terms on the RHS of equation (6.56) and (6.57) could reach of magnitude of O(M) with respect to the characteristic value ρ ĉ. With the definition in (6.56) and (6.57), the acoustic impedance constraints in (5.10) is sufficient to preserve the density positiveness. To simplify the problem, we take the same acoustic impedance conditions as in (5.10).

As the same acoustic impedance is employed for low Mach schemes with AW, PC or the new correction, the same constraints on time step ∆t is adopted. These three low Mach corrections have almost no difference in CPU time, in the next section, we will check the performance in terms of accuracy.

Numerical validation of low Mach corrections

As presented in precedent sections, the new low Mach correction gives a uniform truncation error with respect to Mach number M, density ratio α m and sound speed ratio β m , while the acoustic impedance weighted type or centered pressure type corrections could not satisfy the necessary conditions for asymptotic-preserving property. To examine the numerical behavior for two-phase flow, we implement two classical test cases: the two-phase Gresho vortex and the static bubble.

Gresho vortex problem

The Gresho vortex [START_REF] Gresho | On the theory of semi-implicit projection methods for viscous incompressible flow and its implementation via a finite element method that also introduces a nearly consistent mass matrix. Part 2:Implementation[END_REF][START_REF] Miczek | New numerical solver for flows at various Mach numbers[END_REF] is a rotating flow with a time-independent solution of homogeneous Euler equations. Here we adapt this problem to two-phase flow. A quasi-incompressible fluid 1 is placed inside a quasi-incompressible fluid 2, with a vortex centered at (x, y) = (0.5, 0.5). By using polar coordinates (r, ϑ ), the density distribution can be expressed as:

ρ(r, ϑ ) =    ρ 0 1 r < 0.2, ρ 0 2 r > 0.2, (6.58)
We keep the same angular velocity distribution as proposed in [START_REF] Miczek | New numerical solver for flows at various Mach numbers[END_REF]:

(u r , u ϑ ) =          (0, 5r) 0 ≤ r < 0.2, (0, 2 -5r) 0.2 ≤ r < 0.4, (0, 0) r ≥ 0.4.
(6.59)

The velocity reaches its maximum value of u ϑ ,max = 1 at the interface between two different fluids (r = 0.2). As the centrifugal force is exactly balanced by pressure gradient, the pressure distribution can be given by: The reference pressure p 0 satisfies:

p(r, ϑ ) =          p 0 + 12.5ρ 1 r 2 0 ≤ r < 0.
p 0 = ρ 0 1 γ 1 M 2 1,max -π ∞ 1 = ρ 0 2 γ 2 M 2 2,max -π ∞ 2 , (6.61)
where π ∞ 1 and π ∞ 2 are two constant pressures related to the EOS of fluids 1 and 2. As presented in [START_REF] Miczek | New numerical solver for flows at various Mach numbers[END_REF] where two fluids with the same fluid property are considered. Problem #2 is dedicated to testing the ability of the low Mach correction in the case of large density ratio, while the initial sound speeds of different phases remain the same. Problem #3 is devoted to checking the ability of the low Mach correction that takes account of a large Mach variation, while the initial density remains the same at each phase. Problem #4 is a test case with both high density ratio and large Mach number variation. For the Problem #4, the low Mach scheme with AW or CP type correction provides the same results since the initial acoustic impedance ratio of two different fluids is set to 1. A global view of the vortex at t = 0.1 is presented in Fig. 6.5(b) or Fig. 6.5(c), we can find that the vortex is maintained. With the velocity magnitude profile on the horizontal line passing through the center of the vortex given in Fig. 6.6, we can see a more quantitative result: compared with the new correction, the velocity field of either AW or CP type correction is no more smooth, the vortex is distorted.

From the results of the Gresho vortex problem we conclude that, with a qualitative analysis, either AW or CP type low Mach scheme cannot preserve the two-phase vortex with large density ratio as a truncation error of magnitude of O( 1

α m
) is involved. This test case is not sensitive, even if the asymptotic-preserving conditions such as a constant thermodynamic pressure are not strictly respected, the low Mach scheme can also preserve the vortex.

Static bubble

To indicate the influence of existing low Mach corrections on scheme temporal stability with capillary effects, the "Static bubble" test case is certainly the first case of rudimentary simulation to implement. The simplicity of this case of simulation makes it possible to well isolate the capillary phenomenon, the theoretical solution being simply dictated by Laplace's law. A 2D bubble with a radius of R = 0.4 containing fluid 1 is placed in fluid 2, which are both incompressible. The initial pressure into fluid 2 is p 0 2 = 1. To satisfy the Laplace's law, the initial pressure in fluid 1 is given by: p 0 1 = p 0 2 + σ /R 0 . The surface tension coefficient of σ is set to 0.6144, and the Laplace number is defined as La = 2ρRσ /µ 2 . For problems with different densities in each phase, we keep the same dynamic viscosity and different Laplace numbers in each fluid. Problems #5 -#8 are designed with different fluid properties (see Table 6.2). With the same design as the Gresho vortex problem, Problem #5 considers two fluids with the same EOS and initial density. Problem #6 is devoted to checking the ability of low Mach correction that take into account large density ratio. While Problem #7 takes into account the variation of Mach number while the initial density of each phase remains the same. Problem #8 is a more complicated case with a high density ratio and a large Mach number variation.

The existence of discretization errors on the curvature act as an initial impulse that makes the bubble oscillating. With the effect of physical viscosity, the velocity amplitude of bubble oscillation decreases and the bubble should return to rest. In order to test the scheme numerical stability, configurations with high Laplace number should be considered, corresponding to low physical viscosity. In the present study, the Laplace is set to around 12000. As presented in [START_REF] Popinet | Numerical models of surface tension[END_REF], the redistancing step tends to perturb the curvature and prevents the system from reaching an exact balance. For this test case, this step is not activated.

Numerical results for Problems #5 -#8 are presented in Fig. 6.7. In this figure, we can find that the numerical results of some configurations are absent. From local zoom, we can find that these absences correspond to the divergence of the numerical result, the scheme stability is not ensured. With results presented in Fig. 6.7, we can get a global conclusion: numerical scheme without correction and with new low Mach correction converge and the velocity fluctuations are reduced to machine precision. Low-Mach scheme with AW or CP type correction can not always converge and maintain the scheme stable. For configurations in which the velocity fluctuations are reduced to machine precision, the low Mach scheme can always reduce the velocity fluctuations at a lower rate than without correction, proving that the numerical dissipation is clearly weakened.

For two fluids with the same EOS and initial density (Problem #5), all of the present low Mach corrections show good behavior of recovering the machine precision fluctuations. But the numerical results of these low Mach corrections are a little different as the pressure jump across the interface leads to different sound speeds.

For Problem #6 with large density ratio, the acoustic impedance weighted type and the new low Mach corrections can maintain a good convergence, as they provide the same truncation error for the momentum update which is independent of density ratio α m . While the CP type correction involves a truncation error of magnitude O( 1 α m ), the numerical result shows an immediate divergence.

When the velocity fluctuations amplitude is reduced to around machine precision, the low Mach scheme with AW type correction starts to diverge. Compared to the new low Mach correction, the AW type correction provides an extra truncation error of magnitude of O(1) on mass and energy update associated with u * , showing the importance of the correction on u * on preserving constant thermodynamic pressure. While for the Problem #7, we can observe similar results as Problem #6. At the beginning of the simulation, the numerical scheme with CP type and new low Mach correction can decrease velocity fluctuations. Compared to the CP type correction, the AW type correction provides an extra error of magnitude of O(1) related to the first term of RHS of π * ,θ ,AW in (6.29) as presented in the second line of (B.2), the numerical convergence is not observed. When the magnitude of dimensionless velocity fluctuations is decreased to 10 -8 , the phenomena of divergence is observed as the necessary condition of not satisfied according to the asymptotic analysis in (6.43).

As for Problem #8, both AW and CP types corrections provides the truncation error of magnitude

O( 1
α m ).
As expected, these two low Mach corrections cannot give satisfactory numerical results.

Conclusion

In this chapter, we studied the simulation of a compressible two-phase flows model with sharp interfaces in the low Mach regime. We proposed a definition of the low Mach regime that accounts for possible great discrepancies between the orders of magnitude that are associated with the fluid properties of each phase. We showed that classic low Mach fixes of the literature for Godunov-type solvers provide accurate results on Cartesian grids far from the interface. However, focusing on the numerical method proposed in [START_REF] Zou | An accurate sharp interface method for two-phase compressible flows at low-mach regime[END_REF], we analyzed the behavior of the scheme in the vicinity of the interface. We found out that the abrupt changes of order of magnitudes of the flow parameters across the interface dramatically impact the precision of the method. Indeed, one could see that for an initial condition that verifies the low Mach regime defined in this work, the computed solution will rapidly stop to satisfy the low Mach regime requirements in the vicinity of the interface. We then proposed two new numerical solvers that can cope with these important variations and that preserve low Mach regime approximated solutions. Numerical results were presented for the Gresho vortex test-case that showed good performance of the proposed low Mach correction. The new methods were able to successfully capture static bubble equilibrium with a high Laplace number.

Chapter 7

Numerical results

After previous verification of the correct implementation of the Level Set advection and analysis on low Mach correction, we now assess the performance of our two-phase flow solver on verification test cases. The approximate Riemann solver W N RP (6.55) with the new low Mach correction is employed for all test cases. To implement numerical experiments and validate our numerical model, we develop software for two-phase flows by using the C++ programming language. To enable portability, we implement the Kokkos library which is the C++ library for performance portability. As shown in Fig. 7.1, with the Kokkos library, the software could work on different architectures with or without accelerator. The software works on a high-performance machine "Poincaré" which is hosted by the Maison de la Simulation. There are 92 CPU nodes and 4 GPU nodes in Poincaré. Each CPU node is composed of 2 processors Sandy Bridge E5-2670 (2.60 GHz, 8 cores per processor, 32 GB of memory per node) while each GPU node is composed of 2 processors Sandy Bridge E5-2670 and 2 GPU Tesla K20.

Test case "Two-dimensional sloshing" could be applied to simulate maritime transport of liquefied natural gas. Test case "Rayleigh-Taylor instabilities" is employed to evaluate the performance of our model with large interface deformation. Test case "Static bubble in equilibrium" is utilized to assess the performance of our model with the presence of surface tension effects. Test case "Rising bubble" is frequently encountered in industrial processes. With this test case, we can study a complex flow with surface effects. Test case "1D non-isothermal problem" is set to check the good implementation of thermal effects. All test cases with phase change including "Stefan problem", "Sucking problem" and "Growing vapor bubble in a superheated liquid under zero gravity condition" are employed to validate our phase change model in 1D and 2D respectively.

Two-dimensional sloshing

In this configuration, two non-miscible inviscid fluids of different densities ρ 1 and ρ 2 with ρ 2 > ρ 1 are initially at rest in a rectangular tank, the lighter fluid being over the heavier one (see Fig. 7.2). Gravity g is acting in the vertical downward direction. Then the tank is subjected to a constant horizontal acceleration a 0 with a 0 /g = 0.01. Initial pressure distribution is hydrostatic (reference p 0 = 10 5 Pa). Gas is described by an ideal gas law with density ρ 1 = 1 kg m -3 and γ 1 = 1.4. Sound speed in the gas is around 370 m s -1 . Liquid is described by stiffened gas law, with density ρ 2 = 1000 kg m -3 , γ 2 = 7 and p ∞ 2 = 3.31 × 10 8 Pa. Sound speed in the liquid is around 1500 m/s. Gas and liquid heights are 1.25L and L respectively with L = 1m.

Since a 0 /g ratio is small, the interface position can be given analytically by a linear potential approach [START_REF] Chanteperdrix | Modélisation et Simulation Numérique des Ecoulements Diphasiques à Interfaces Libres. Application à l'étude des mouvements de liquides dans les réservoirs de véhicules spatiaux[END_REF]:

ξ (x,t) = a 0 g x - L 2 ∑ n≥0 4 Lk 2 2n+1 cos(ω 2n+1 t)cos(k 2n+1 x) ,
where ξ (x,t) is a function which give the theoretical interface position, k n = πn L is the wave number, and ω n can be given by:

ω 2 2n+1 = gk 2n+1 (ρ 2 -ρ 1 )
ρ 1 coth(k 2n+1 h 1 ) + ρ 2 coth(k 2n+1 h 2 ) . The evolution of computed interface positions on the left and right sides of the tank are compared to analytical solutions, at different spatial resolutions (Fig. 7.3 and Fig. 7.4). We observe that the present numerical scheme without low Mach correction is highly dissipative and cannot describe the interface. This is mainly due to the intrinsic dissipation of the scheme (first-order). A posteriori, Mach number of this flow is found to be around 2 × 10 -5 : in this low Mach regime, an adequate numerical scheme is mandatory. Proposed low Mach correction eliminates excessive numerical dissipation (even without any MUSCL reconstruction) and results coincide with the analytical solution even on a rather coarse mesh.

Rayleigh-Taylor instabilities

A test case involving more severe deformation of interface is investigated: "Rayleigh-Taylor instabilities". An accurate numerical description of the interface is needed to follow the complex evolution of interface within this flow. The computation domain is rectangular and filled by two non-miscible fluids of different densities ρ 1 = 1 and ρ 2 = 1.8. The two layers of fluid are initially superimposed with the heavier one over a lighter one. Initial interface is perturbed and located at y = 1sin(2πx). Reynolds . Both models converge but the way that they converge is quite different. Numerical results show that the incompressible method converges faster than the compressible method (which is expected as the compressible method is first-order while the incompressible method is second-order) but the main difference exists in the area of high gradients.

From this test case, we can find that in a low Mach regime, our numerical scheme can well predict the fluid behavior with large deformation.

Static bubble in equilibrium

This test case is the same configuration as the test case "Static bubble" presented in Section. 6.5.2. A 2D bubble with a radius R=0.4 containing an ideal gas is placed in a slightly compressible liquid. The gas is described by an ideal gas law with γ 1 = 1.4 and p ∞ 1 = 0 while the liquid is described by a stiffened gas law with γ 2 = 7.14 and p ∞ 2 = 300. The initial densities of these two fluids are equal with ρ 1,0 = ρ 2,0 = 1. The initial pressure in liquid is p 2,0 = 1. To confirm the Laplace's law, the initial pressure in the gas is given by: p 1,0 = p 2,0 + σ R 0 . The surface tension coefficient σ is set to be 1, and the Laplace number La = ρD 0 σ µ 2 is set to be 12000. show the spatial convergence of numerical scheme resolutions without and with low Mach correction respectively. From these two figures, we can find that the initial maximum amplitudes of spurious currents decrease as we refine mesh which shows a good convergence of the numerical scheme. Because of the total viscosity (including numerical viscosity and physical viscosity µ), both of these two models reduce the amplitude of spurious currents to machine precision and a good well-balanced property is observed.

In Fig. 7.8, without low Mach correction, the amplitudes of spurious currents reduce at different rates. As mesh refinement has an important influence on numerical viscosity which is proportional to the mesh size, these different rates at different resolutions show the importance of the numerical viscosity compared to the physical viscosity µ. In Fig 7 .9, with low Mach correction, the amplitudes of spurious currents reduce at the same rate. This phenomenon shows that the insignificance of numerical viscosity compared to the physical viscosity and the numerical viscosity is eliminated effectively.

Rising bubble

In this section, we want to study a more complex flow with surface tension effects, viscous effects, high-density ratio and large interface deformation. For this first problem, we simulate a configuration that has been studied with an incompressible Level Set method [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF]. It is an air bubble rising in a column of water at rest. The domain is shown in Fig. 7.10. The flow is characterized by the following dimensionless numbers: the Reynolds number Re = (2R) 3 gρ 2 /µ 2 = 1000 and the Bond number Bo = 4ρ 2 R 2 g/σ =200 (the importance of gravitational forces compared to surface tension force). The Numerical results are presented in Fig. 7.11, the evolution of the interface is shown and superimposed with reference results [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF] (incompressible Level Set method). We can find that these two models have similar behavior and the top of interfaces of these two models coincides very well. From t * = t/ g/R = 4.8, in the incompressible Level Set formulation, there are some small bubbles that detach and appear, while with the present model, this phenomenon is not predicted. This test case is constructed numerically, it is hard to reproduce experimentally. For the phenomenon such as the bubble coalescence and separation, more special treatments are required, such as the curvature estimation. As this problem is not the key issue of the present study, more attention should be paid when studying bubbly flow.

To show the performance of the present software in terms of efficiency, we give some numerical results of this test case. The software works on Poincare with one node (16 cores) with OpenMP parallelism. The whole simulation takes 23042 seconds and gives 8.189 million cells update per second.

1D non-isothermal problem

With this test case, we underline the advantage to simulate low Mach flows with a compressible solver. When considering such flows, compressibility effects could be important when considering thermal dilatation. We consider the test case "1D non-isothermal problem" as in [START_REF] Daru | A numerical method for the simulation of low mach number liquid-gas flows[END_REF]. A closed 1D tube of length 100 µm is consisting of a liquid layer of length 10 µm between of two layers of air as depicted in Fig. 7.12. The liquid layer is initially suited at the center of the tube. The initial pressure and Fig. 7.12 Sketch of the 1D non-isothermal tube temperature in the system are 101325 Pa and 293.15 K respectively. Both fluids are defined by the following characteristics: for gas, γ 1 = 1.4, p ∞ 1 = 0, specific heat at constant pressure C p,1 = 1004.5 J K -1 kg -1 , thermal conductivity K 1 = 0.0256 W m -1 K -1 , dynamic viscosity µ 1 = 1.82 × 10 -5 Pa • S. For liquid, a new equation of state is considered, which is barotropic (pressure depends only on density, for more details please refer in Appendix. A), with following properties: sound speed 1500 m s -1 , initial density ρ 0,2 = 1000 kg m -3 , specific heat C p,2 = 4184 J K -1 kg -1 , thermal conductivity K 2 = 0.6 W m -1 K -1 , dynamic viscosity µ 2 = 0.001 Pa • S. At initial time, the left wall is heated to T w = 373.15 K and the right wall is insulated.

The domain resolution consists of 100 grid points. The tube is then heated, and the liquid layer start moving due to dilatation of left gas layer. Fig. 7.13 presents the right interface trajectory in the early stages after the heating of the left wall. An oscillatory motion is observed. This oscillation is entirely related to the gas compressibility as presented in [START_REF] Daru | A numerical method for the simulation of low mach number liquid-gas flows[END_REF]. This phenomenon confirms the importance of compressibility for low Mach flows prediction. Fig. 7.14 shows the velocity profile at t = 1m s, while Fig. 7.15 and Fig. 7.16 present the temperature profile at t = 1 µs and t = 4 ms respectively. By comparing the interface trajectory, velocity profile and temperature profile with the reference [START_REF] Daru | A numerical method for the simulation of low mach number liquid-gas flows[END_REF], a good agreement could be found.

Test cases with phase change

To validate our phase change model, several test cases such as the "Stefan problem", "Sucking problem" and "Growing vapor bubble in a superheated liquid under zero gravity condition" should To simulate these incompressible phase change problems, both the liquid and vapor are described with a barotropic EOS, physical properties of liquid and vapor are given in Table . 7.1. To compare 7.1 Physical properties of phase change problems our numerical results with incompressible analytical solutions, we set that the fluid (vapor) at rest as less compressible (with higher sound speed) than the liquid.

(kg/m 3 ) (Pa • s) (J /kgK ) (W /mK ) (m/s) Vapor ρ 0,v µ v C p,v K v c 0,v 0.

Stefan problem

The Stefan problem is defined as in Fig. 7.17. Vapor and liquid are contained in a free-slip tube. The left of the tube is closed by a heated wall at a constant temperature T wall , which is higher than the saturation temperature T sat . The right side of the tube is opened, and treated as an outlet. Initially, the vapor is at rest, there is a linear temperature profile in the vapor, and a constant temperature profile in the liquid. The tube is then heated by the left wall, the liquid will boil at the interface, the interface will move toward the outlet. The interface position X(t) and the temperature profile T (x,t) are theoretically given by [START_REF] Sato | A sharp-interface phase change model for a mass-conservative interface tracking method[END_REF]:

X(t) = 2X K v ρ v C p,v t, (7.1)
T (x,t) = T wall + T sat -T wall erf(X) erf( x 2 K v /(ρ v C p,v )t ), (7.2)
where erf(x) is the error function, and X is the solution of the transcendental equation:

X exp X 2 erf(X) = C p,v (T wall -T sat ) √ ΠL heat . ( 7.3)
In the present test case, the wall temperature is set at 383.15K, the tube is heated under atmospheric pressure. The test is computed with four different meshes: ∆x = 2 × 10 -5 m, ∆x = 1 × 10 -5 m, ∆x = 5 × 10 -6 m and ∆x = 2.5 × 10 -6 m.

The interface evolution is shown in Fig 7 .18. We can find that the numerical results of the coarsest mesh is slightly different from the theoretical solution. The other numerical results coincide well with the analytical solution, a good convergence can be found between different resolutions. The temperature profile at t = 0.1 s is presented in Fig. 7.19, numerical results on temperature distribution show a good agreement with the analytical solution. Fig. 7.20 indicates the error in mass conservation in dimensionless form. At the end of the simulation, the error of the coarsest mesh is around 5 × 10 -4 . The mass conservation is the common issue of the Level Set method. The relative mass conservation error of the finest mesh is around 10 -5 , the mass conservation property is largely improved with mesh refinement.

Sucking problem

The sucking problem is schematically in Fig. 7.21 which is similar to the Stefan problem but with super-heated liquid. This test case is proposed by Welch and Wilson [START_REF] Welch | A volume of fluid based method for fluid flows with phase change[END_REF] to verify the phase change outlet. Similar to the Stefan problem, its analytical solution can be given as [START_REF] Fleau | Multifield approach and interface locating method for two-phase flows in nuclear power plant[END_REF][START_REF] Welch | A volume of fluid based method for fluid flows with phase change[END_REF]:

X(t) = P ′ ρ l ρ v 2 K l ρC p,l t, (7.4)
and

T = T sat + L heat C p,l P ′ exp x 2 0 Π 2 erf(x 0 , x), (7.5)
with x = x -X(t) + P ′ √ 2 and x = P ′ √ 2 , with P ′ can be given by the following equation:

P ′ exp P ′2 2 = 2 Π (T ∞ -T sat )C p,l L heat . (7.6)
This test is computed with three different resolutions: ∆x = 5 × 10 -5 m, ∆x = 2.5 × 10 -5 m and ∆x = 1.25 × 10 -2 m. Fig. 7.22 shows the interface position evolution while Fig. 7.23 gives the temperature distribution at t = 0.6 s. By comparing the numerical solution with the analytical solution, a good agreement can be found on both interface position and temperature distribution.

To observe the mass conservation property of the present model, by defining the total mass m total as: The total mass variation is shown in Fig. 7.24. From Fig. 7.24, we can find that the relative mass variation is around of magnitude 10 -4 , the mass conservation problem is a well-known problem of the Level Set method, with mesh refinement, the mass conservation property is largely improved.

m total = n cell ∑ n=1 ρ n V cell + t ∑ t ′ =0 ∑ Outlet ρuS∆t. ( 7

Growing vapor bubble in a superheated liquid under zero gravity condition

We consider a two-dimensional problem. Similar to the sucking problem, a bubble of vapor at saturation temperature is placed in a super-heated liquid as depicted in Fig. 7.25. The liquid will boil at the interface and the radius of the bubble increases. The analytical solution of the bubble radius R evolution is given as [START_REF] Sato | A sharp-interface phase change model for a mass-conservative interface tracking method[END_REF]:

R = 2X g K l C p,l ρ l t, (7.8)

Outlet

Outlet Outlet Outlet Fig. 7.25 Condition of simulation at a given time for growing vapor bubble in a superheated liquid under zero gravity condition where X g is the "growth constant" which can be obtained from:

1

M g = 2X 2 g 1 0 exp -X 2 g (1 -x) -2 -2(1 - ρ v ρ l x -1) dx, (7.9)
with M g = ρ v (L + (C p,l -C p,v )(T ∞-T sat )) ρ l C p,l (T ∞ -T sat ) . The analytical solution of the temperature distribution at position r and time t can be given as:

T =      T ∞ -2X g M g 1 1-R/r exp -X 2 g (1 -x) -2 -2(1 - ρ v ρ l x -1) for r > R,
T sat for r ≤ R.

(7.10a) (7.10b)

The degree of superheat is non-dimensionalized by the Jakob number:

Ja = ρ l C p,l (T ∞ -T sat ) ρ v L heat . (7.11)
The simulation domain is a box of 3 mm wide with boundary condition of outlet at each side. The simulations start at a time t 0 with a bubble radius equal to 1 mm, until a time 1.21 t 0 required for the final bubble radius to be 1.1 the initial radius (1.1 mm). Four Jakob numbers Ja = 12, 18, 24, 30 are tested which correspond to ∆T = T ∞ -T sat = 4, 6, 8, 10 K respectively. All tests have been implemented with three different meshes: 64 × 64 (coarse), 128 × 128 (medium), 256 × 256 (fine). Fig. 7.26 and Fig. 7.27 show the temperature field with the finest mesh at the initial and at the final time for a Jakob number equals to 12 and 30 respectively. For the finest mesh, the theoretical thermal boundary layer thickness for the initial bubble radius corresponds to 5∆x and 3∆x respectively. We can find that for the same bubble radius, thermal boundary layer thickness decreases as the Jakob number increases. Fig. 7.28 presents the velocity fields at t = 1.21t 0 with the finest mesh for Ja = 12 and Ja = 30 respectively. Radial expansion flow around bubble can also be visualized, the fluid velocity in the bubble is also null for both Ja = 12 and Ja = 30. Fig. 7.29a -Fig. 7.32a show the bubble radius evolution. By comparing the numerical results with analytical solutions, we can get a conclusion: the coarse mesh gives a result different from the theoretical solution because the boundary layer corresponds to only 1-2 ∆x in that mesh, but a good agreement is achieved for the medium and fine grids. It appears that the convergence rate is not monotonic because of the compressible effects. The acoustic phenomena always exist, especially, a simple velocity extension is used which can not strictly guarantee a divergence-free velocity field in 2D cases. This problem has been observed in [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF].

To visualize the temperature distribution in the thermal boundary layer, in Fig. 7.29b -Fig. 7.32b, the temperature distributions versus r at the final time for different Ja numbers are presented. For all numerical results, temperature distributions given by the coarse mesh are different from the analytical 2D "Growing vapor bubble in a superheated liquid under zero gravity condition" problem, a good agreement could be achieved with medium and coarse meshes. However, because of the presence of acoustic effects, the convergence rate is not monotonic. Although the convergence rate is not monotonic, the good agreement between numerical results with the medium and fine meshes can validate our phase change model in 2D cases.

Chapter 8

Conclusion and Perspective

The aim of this thesis is to develop a code to simulate compressible low Mach flows with phase change. First, a state of the art of the various existing methods for interface tracking and compressible flows was carried out. The advantages and disadvantages of interface description methods were presented. The Level Set method is chosen because of its advantages.

Considering the interface temporal evolution, there are two issues: Level Set advection and redistancing. For the first issue, a high-order coupled time-space approach called "One-Step" (OS) is proposed. Two classical separate time-space approaches WENO-RK and HOUC-RK are also implemented. The global stencil of a separate time-space scheme is much larger than that of a coupled time-space approach. Numerical tests on Level Set advection show the efficiency and accuracy of the OS scheme compared to classical schemes such as WENO-RK and HOUC-RK.

During the advection of the Level Set function, if the velocity field does not impose a rigid translation or rotation, the different contour lines are advected differently, the Level Set function will lose the distance property which will impact the accuracy of interface curvature estimation. We present two different high-order numerical procedures for redistancing: PDE (Partial differential equation) and Hopf-Lax formula based redistancing procedures. The PDE based method tries to reach a steady-state by solving a PDE equation, while the Hopf-Lax formula based method tries to find the minimum distance of a grid point to the interface by using the secant method. When the initial Level Set function is not convex, the secant method may converge to a wrong minimum distance to the interface, a procedure called "check procedure" is necessary. Numerical tests on redistancing show that without the check procedure, the Hopf-Lax based method has an advantage in terms of accuracy and efficiency. With the check procedure, the CPU time of Hopf-Lax based method largely increases, while its accuracy is always better than the PDE based method. Considering the PDE based method is more adapted to the Level Set advection procedure and easy to parallelize, as a reason, the retained scheme for Level Set redistancing is the PDE based method.

To take the capillary effects into account, the interface curvature is necessary. The classical Level Set (CLS) method and Height function method are tested. Numerical results show that the Height function method could provide better results with coarse meshes: better precision and able to estimate higher interface curvature. However, Height function methods are limited to the first-order accuracy while the classical Level Set method can reach theoretical second-order accuracy. When we refine the mesh, the classical Level Set method could provide results with better precision. Considering that the classical Level Set method is easy to implement and converges more quickly than Height function methods, the retained scheme for estimation is the classical Level Set method.

The physical model to describe two-phase flows verifies the compressible Navier-Stokes equations including effects of viscosity, external forces, heat transfer. Effects related to the surface tension and phase change are considered as jump conditions across the interface. For any system, to treat any operator with the most adapted method, it is important to use a technique called "Operator splitting". We give the principle of operator splitting and high-order methods. Then a splitting method called "Lagrange-projection" to treat the Euler system is presented, this splitting method has an advantage in treating interface jump condition. However, a complicated projection step is involved. To avoid this projection step, we adopt the acoustic-transport splitting, which splits the acoustic and transport phenomena separately. With the acoustic-transport splitting, the global model to describe two-phase compressible is decomposed into three subsystems: acoustic, transport and diffusion. The external force terms are added into the acoustic subsystem to get a well-balanced treatment.

Split subsystems are approximated with the finite volume method. An approximate Riemann solver is developed according to the acoustic subsystem. This approximate Riemann solver takes external forces as well as the velocity and pressure jumps associated with phase change and surface tension into account with a well-balanced discretization. The coupling across the interface is realized by the ghost-fluid method. As for the phase change model, to take the interface temperature as a boundary condition, in the literature [START_REF] Sato | A sharp-interface phase change model for a mass-conservative interface tracking method[END_REF][START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF], an implicit method is always applied for the stability issue. To avoid this implicit procedure, in this chapter we propose an explicit method to approximate the velocity jump. For any cell near the interface, we use another neighbor cell belonging to the same fluid to estimate the heat flux that takes the interface as a Dirichlet boundary condition.

Unfortunately, in the case of single-fluid flows, it is now well-established that standard Godunovtype methods may dramatically lose their accuracy in low Mach regime and when the computation domain is discretized over a quadrangle (resp. hexahedral) grid in 2D (resp. 3D). By implementing these low Mach corrections [START_REF] Chalons | An all-regime Lagrange-Projection like scheme for 2D homogeneous models for two-phase flows on unstructured meshes[END_REF][START_REF] Peluchon | A robust implicit-explicit acoustic-transport splitting scheme for two-phase flows[END_REF] into our model, we find that the numerical scheme is difficult to converge. With the analysis of flow property, we find that in low Mach regime, the pressure perturbation is dominated by the dynamic pressure. Moreover, the low Mach scheme should be able to preserve constant thermodynamics pressure in absence of thermal effects. With these arguments, we propose a new low Mach correction. With these features, we propose a new low Mach correction. With the "Gresho vortex" test case, we find that the new low Mach scheme could maintain the vortex in presence of high density ratio and high sound speed ratio. Numerical results of "Static bubble" show good property of convergence of the new low Mach scheme.

The present numerical model is validated with well-known test cases. These test cases are related to different important issues: large interface deformation, convergence of numerical scheme with capillary effects in presence of high Laplace number, dilatable effects and phase change. Numerical results provide a good agreement with analytical solutions or results in the literature, showing the good performance of the proposed model for two-phase low Mach flows.

Perspectives

Still, our model needs further improvements, in particular with the shortcoming of the Level Set method: mass conservation. This problem could be circumvented by coupling to the Volume of fluid method [START_REF] Ménard | Coupling level set/VOF/ghost fluid methods: Validation and application to 3D simulation of the primary break-up of a liquid jet[END_REF] or by using the cut-cell method [START_REF] Hu | A conservative interface method for compressible flows[END_REF]. The coupling with VOF method could involve the well-known problem such as interface reconstruction and the diffusion in the volume fraction transport while the cut-cell approach involves complex geometry and small cut cells. In this work, the precision of the acoustic and transport (except Level Set advection) subsystem estimation is limited to first-order accuracy. A higher-order scheme both in time and space could be considered, the reconstruction across the interface could be a big challenge as the interface is treated as a contact discontinuity.

For a low Mach flow simulated by a compressible solver, despite the excessive numerical diffusion, there is still another issue: stability condition requires very restrictive time steps for explicit schemes due to the importance of acoustic speed of sound compared to advective velocity. The semi-implicit could be interesting to implement. An implicit update for the acoustic system could significantly decrease the CPU time. When we refine the mesh, the time step could be limited by the diffusion subsystem, an implicit update on diffusion subsystem could be also investigated. The convergence of linear solver of two-phase flows with low-Mach correction could be far more difficult especially it could involve two different Mach numbers, high sound speed ratio and high density ratio.

Considering application into industrial processes, more complicated thermodynamic models should be studied. According to the fluid states close to the interface, the varying equilibrium state (saturation temperature) could be provided. Résumé: Une approche d'interface raide est présentée pour le calcul des écoulements diphasiques avec tension superficielle et changement de phase en régime à faible nombre de Mach. Pour développer un tel modèle, où de légers effets compressibles sont pris en compte ainsi que des fermetures thermodynamiques correctes, le liquide et le gaz sont considérés comme compressibles et décrits par un solveur compressible précis. Ce solveur compressible adopte une technique de décomposition appelée "décomposition du transport acoustique" qui décompose le système Euler en deux parties: acoustique et transport. Sur la base du sous-système acoustique, un solveur de Riemann approximatif qui tient compte des effets de tension superficielle et de changement de phase est développé. L'interface de l'écoulement diphasique est capturée par la méthode de Level Set et considérée comme raide. La problème de la capture d'interface de la méthode Level Set dans le cadre Eulérien est le point clé des simulations d'écoulement diphasique, et dans ce travail, nous proposons et adoptons des approches d'ordre élevé pour l'advection de l'interface, la redistanciation et l'estimation de la courbure. En régime à faible nombre de Mach, les solveurs compressibles conventionnels perdent en précision et une correction à faible Mach est alors nécessaire pour réduire la dissipation numérique. Pour une méthode d'interface raide, l'interface est traitée comme la discontinuité de contact via la méthode Ghost Fluid. Sans une région lisse à l'interface, une telle discontinuité existant à l'interface présente un énorme défi pour la conception d'un schéma numérique. La correction à faible Mach bien connue dans la littérature pourrait conduire à une erreur de troncature significative, en particulier pour les écoulements diphasiques avec de grands rapports de densité et de vitesse du son. Pour retrouver une bonne propriété de préservation asymptotique, nous proposons une nouvelle correction à faible Mach avec une analyse asymptotique rigoureuse. Plusieurs cas de test numériques ont été utilisés pour valider la présente approche numérique et montrer ses bonnes performances.

for computing two-phase flows with surface tension and phase change in low Mach regime. To develop such a model, where slight compressible effects are taken into account as well as correct thermodynamical closures, both the liquid and the gas are considered compressible and described by a precise compressible solver. This compressible solver adopt a splitting technique called "acoustic-transport splitting" which splits the Euler system into two parts: acoustic and transport. Based on the acoustic subsystem, an approximate Riemann solver that accounts for surface tension and phase change effects is developed. The interface between two-phase flows is captured by the Level Set method that is considered to be sharp. The interface capturing issue of the Level Set method within the Eulerian framework is the key point of the two-phase flow simulations, and in this work we propose and adopt high-order approaches for interface advection, redistancing and curvature estimation. In low Mach regime, conventional compressible solvers lose accuracy and a low Mach correction is then necessary to reduce the numerical dissipation. For a sharp interface method, the interface is treated as the shock-wave contact discontinuity via the Ghost Fluid method. Without a smooth region at the interface, such discontinuity existing at the interface presents a huge challenge to the design of a numerical scheme. The well-known low Mach fix in literature could lead to significant truncation error, especially for two-phase flows with large density and sound speed ratios. To recover a good asymptotic-preserving property, we propose a new low Mach correction with rigorous asymptotic analysis. Several numerical test cases have been employed to validate the present numerical approach and enlighten its good performance.

Université Paris-Saclay Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France

1 Introduction 1 . 1 1 . 1 5 . 4

 1111154 Industrial context . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Scientific context . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 The Interface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Sharp and diffuse interface methods . . . . . . . . . . . . . . . . . . . . . . 1.2.3 Lagrangian and Eulerian approaches . . . . . . . . . . . . . . . . . . . . . . 1.2.4 Simulation of low Mach two-phase flows . . . . . . . . . . . . . . . . . . . 1.3 Outline of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .2 Global Model 2.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Interface jump conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 Jump condition on mass . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.2 Jump condition on stress . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.3 Jump condition on energy . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Global model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Operator splitting 3.1 Basic idea . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 High-order splitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Euler system splitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Lagrange-Projection method . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Acoustic-transport splitting . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Decomposition of the global system . . . . . . . . . . . . . . . . . . . . . . . . . . Representation of a diffuse interface (left) and a sharp interface (right). . . . . . . . . 1.2 Representation of the interface Γ and two sub-domains D 1 ⊂ R d and D 2 ⊂ R d occupied by fluid 1 and fluid 2 respectively. . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Representation of Lagrangian approach with mobile mesh, dash line: the interface. . 1.4 Representation of Front-tracking mesh with fixed mesh and markers to represent the interface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Representation of Volume Of Fluid method. dash line: reconstructed interface . . . . 1.6 Representation of Level Set method. Blue dash line: interface; black dash lines: isocontours of φ = ±d. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Grids evolution of Lagrangian-Projection method. green area: Ω E j ∩ Ω L j , blue area: Ω E j ∩ Ω L j+1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Function sin 4 (πx). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Spatial Convergence of global (left) and local (right) L 1 errors for circle advection at t = 100. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Global (left) and local (right) L 1 errors versus computation time at t = 100 for circle advection (each symbol is a different spatial resolution). . . . . . . . . . . . . . . . 4.4 Zero-level contours of the advected circle for different numerical schemes with a mesh 32 × 32. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Results of the Zalesak disk after one full rotation with a resolution of 50 × 50. . . . .4.6 Global (left) and local (right) L 1 errors for the Zalesak disk rotation . . . . . . . . . 4.7 Zero-level contours at t = 3. Resolution 100 × 100 for different numerical schemes and exact solution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Temporal evolution of surface conservation for deformed circle. Different numerical schemes and resolutions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.9 Evolution of some iso-contours of Level Set function for the test case at different times: Severe interface deformation. red lines: iso-contours; black lines: interface; blue region: D 1 (t) = {x ∈ D | φ (x,t) > 0} . . . . . . . . . . . . . . . . . . . . . . . 4.10 Hopf-Lax Redistancing of a non-convex function φ with a mesh 32 × 32. . . . . . . . xv Representation of two different interface geometries. red point: active cell j; blue region: D 1 , green region: D 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Extrapolation of interface velocity u Γ and velocity jump. Colored cells: interface velocity u Γ and velocity jump initialized cells; white cells: undefined cells. . . . . . 6.1 1-D sketch of a cell j across the interface Γ, separating domains D i , with its full neighborhood V and partial neighborhoods N ± ( j). Distance function φ , positive at cell j. Discontinuous pressure field p across the interface with non-dimensional values at cells j and k. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Magnitude of the velocity field of two-phase Gresho vortex problem #1 with a resolution of 100 × 100. (a): without low Mach correction; (b) AW type correction; (c) CP type correction; (d) new low Mach correction . . . . . . . . . . . . . . . . . 6.3 Magnitude of the velocity field of two-phase Gresho vortex problem #2 with a resolution of 100 × 100. (a): without low Mach correction; (b) AW type correction; (c) CP type correction; (d) new low Mach correction. . . . . . . . . . . . . . . . . . 6.4 Magnitude of the velocity field of two-phase Gresho vortex problem #3 with a resolution of 100 × 100. (a): without low Mach correction; (b) AW type correction; (c) CP type correction; (d) new low Mach correction. . . . . . . . . . . . . . . . . . 6.5 Magnitude of the velocity field of two-phase Gresho vortex problem #4 with a resolution of 100 × 100. (a): without low Mach correction; (b) AW type correction; (c) CP type correction; (d) new low Mach correction. . . . . . . . . . . . . . . . . . 6.6 Velocity magnitude profile on the horizontal line passing through the center of the vortex. Red line: AW or CP correction; blue line: new correction. . . . . . . . . . . 6.7 Temporal evolution of dimensionless velocity fluctuations of a 2D bubble (La=12000) at resolution of R = 12.8∆x . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.1 Presentation of Kokkos. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Sketch of the sloshing tank with acceleration vectors. . . . . . . . . . . . . . . . . . 7.3 Temporal evolution of interface at x = 0 and at x = L with a uniform mesh 40 × 90. red line: present method with low Mach correction; blue line: present method without low Mach correction; black line: analytical solution. . . . . . . . . . . . . . . . . . . 7.4 Temporal evolution of interface at x = 0 and at x = L with a uniform mesh80 × 180. red line: present method with low Mach correction; blue line: present method without low Mach correction; black line: analytical solution. . . . . . . . . . . . . . . . . . . 7.5 Vorticity of Rayleigh-Taylor instabilities test case with a resolution of 320 × 640 at t * = 1 (left) and t * = 3 (right), black lines: the interface. . . . . . . . . . . . . . . . 7.6 Vorticity of Rayleigh-Taylor instabilities test case with a resolution of 320 × 640 at t * = 5. Left: present compressible Level Set formulation; Right: incompressible Level Set formulation [109], black lines: interface. . . . . . . . . . . . . . . . . . . List of tables 4.1 Level Set advection L 1 errors for test case #1 sin 4 (πx) by using the separate time-space approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Level Set advection L 1 errors for test case #1 sin 4 (πx) by using the coupled time-space approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Average global and local L 2 errors for test case #1 Circle . . . . . . . . . . . . . . . 4.4 Average global and local L 2 errors for test case #2 Square . . . . . . . . . . . . . . .4.5 Average global and local L 2 errors for test case #3 Ellipse . . . . . . . . . . . . . . . 6.1 Numerical parameters of Gresho vortex test problems . . . . . . . . . . . . . . . . . 6.2 Numerical parameters of static bubble test problems . . . . . . . . . . . . . . . . . . 7.1 Physical properties of phase change problems . . . . . . . . . . . . . . . . . . . . .

nK 1 θ 1 Introduction 1 . 1

 1111 Specific heat capacity at constant volume n jk Vector normal to the face ∂ Ω jk toward to cell k Vector Dimensionaless variable with respect to charateristic value of fluid i Dimensional density of fluid i with respect to characteristic density of fluid Constant value for acoustic impedance definition µ Dynamic viscosity N ( j) Neighbour cell of cell j N + ( j) Neighbour cell of cell j belong to the same fluid N -( j) Neighbour cell of cell j belong to another fluid • Dimensionaless variable ∂ Ω jk Cell face between cells j and k Dimensional sound speed of fluid i with respect to characteristic density of fluid Low Mach correction coefficient, approximation of Mach number Chapter Industrial context

Fig. 1 . 1

 11 Fig. 1.1 Representation of a diffuse interface (left) and a sharp interface (right).

D 1 D 2 ΓFig. 1 . 2

 1212 Fig. 1.2 Representation of the interface Γ and two sub-domains D 1 ⊂ R d and D 2 ⊂ R d occupied by fluid 1 and fluid 2 respectively.

Fig. 1 . 3

 13 Fig. 1.3 Representation of Lagrangian approach with mobile mesh, dash line: the interface.

1 D 2 Fig. 1 . 4

 1214 Fig. 1.4 Representation of Front-tracking mesh with fixed mesh and markers to represent the interface.

Fig. 1 . 5

 15 Fig. 1.5 Representation of Volume Of Fluid method. dash line: reconstructed interface

Fig. 1 . 6

 16 Fig. 1.6 Representation of Level Set method. Blue dash line: interface; black dash lines: isocontours of φ = ±d.

Fig. 4 . 1

 41 Fig. 4.1 Function sin 4 (πx).

Fig. 4 . 2 Fig. 4 . 3

 4243 Fig. 4.2 Spatial Convergence of global (left) and local (right) L 1 errors for circle advection at t = 100.

Fig. 4 . 4

 44 Fig. 4.4 Zero-level contours of the advected circle for different numerical schemes with a mesh 32 × 32.

Fig. 4 . 5

 45 Fig. 4.5 Results of the Zalesak disk after one full rotation with a resolution of 50 × 50.

Fig. 4 . 6

 46 Fig. 4.6 Global (left) and local (right) L 1 errors for the Zalesak disk rotation

Fig. 4 . 7

 47 Fig. 4.7 Zero-level contours at t = 3. Resolution 100 × 100 for different numerical schemes and exact solution.

Fig. 4 . 8

 48 Fig. 4.8 Temporal evolution of surface conservation for deformed circle. Different numerical schemes and resolutions.

  Fig. 4.10 Hopf-Lax Redistancing of a non-convex function φ with a mesh 32 × 32.

  Fig. 4.11 Hopf-Lax Redistancing of a function φ 0 involving d dt φ (x x x i i i ,t) = 0 with a mesh 32 × 32.

  secant method: end // projected gradient descent method: start y y y j k = argmin φ (y y y), ∀y y y ∈ B(x x x i i i ,t ′, k-1 ) ; while |y y y j kx x x i i i | < |t ′, k | do ȳ y y j+1 k = y y y j k -ε∆x ∇φ 0 (y y y j k ) |∇φ 0 (y y y j k )| sign(φ 0 (y y y j k )); y y y j+1 k = PROJ(ȳ y y j+1

2 ) 2 )

 22 Fig. 4.13 The initial function of test case #1 with a mesh 32 × 32. Left: 3D representation Right: isocontour of φ 0 black line: zero-level white line: other arbitrary isocontours

Fig. 4 .

 4 Fig. 4.14 Function after redistancing procedure of test case #1 circle with a mesh 32 × 32. Top: function after PDE based redistancing procedure; Bottom: Hopf-Lax based redistancing procedure; Left: 3D representation; Right: contour of φ 0 (In the right figures, black line: zero-level white line: other arbitrary isocontours)

Fig. 4 .

 4 Fig. 4.15 The initial function of test case #2 square with a mesh 32 × 32. Left: 3D representation Right: isocontour of φ 0 black line: zero-level white line: other arbitrary isocontours

  .

Fig. 4 .

 4 Fig. 4.17 The initial function of test case #3 ellipse with a mesh 32 × 32. Left: 3D representation Right: isocontour of φ 0 black line: zero level white line: other arbitrary isocontours

Fig. 4 .

 4 Fig. 4.18 Function after redistancing procedure of test case #3 ellipse with a mesh 32 × 32. Top: function after PDE based redistancing procedure; Bottom: Hopf-Lax based redistancing procedure; Left: 3D representation; Right: isocontour of φ 0 (In the right figures, black line: zero level white line: other arbitrary isocontours)

  Fig. 4.19 Zero-level contours at t = 3 with a resolution of 128 × 128, Red line: Without redistancing; Blue line: PDE based redistancing; Green line: Hopf-Lax based redistancing; Black line: theoretical solution.

Fig. 4 .

 4 Fig. 4.20 Temporal evolution of surface conservation for deformed circle with PDE based redistancing.

Fig. 4 .

 4 Fig. 4.21 Temporal evolution of surface conservation for deformed circle with Hopf-Lax based redistancing.

Fig. 4 .

 4 Fig. 4.24 Curvature estimation errors of Height function methods and classical Level Set method (CLS). Left: PND method; Right: RND method.

p

  for y ∈ Ω l ∩ Ω r means the jump of the pressure across the cell phase

Fig. 5 . 1

 51 Fig. 5.1 Mesh for a multi-dimensional problem with an intersection point x jk at the cell face.

cell 1 Fig. 5 . 2

 152 Fig. 5.2 Ghost cell representation

Fig. 5 .

 5 Fig. 5.3 representation of temperature profile of two-pase problems. T Γ : interface temperature

3 Fig. 5 . 4

 354 Fig. 5.4 Representation of two different interface geometries. red point: active cell j; blue region: D 1 , green region: D 2 .

Fig. 5 . 5

 55 Fig. 5.5 Extrapolation of interface velocity u Γ and velocity jump. Colored cells: interface velocity u Γ and velocity jump initialized cells; white cells: undefined cells.

( 2 p

 2 H1) the fluid velocity normal to the interface (u + u Γ jk ) • n is continuous across the interface, (H2) The interface pressure jump: = p Γ (y,t), for y ∈ Γ(t), (6.1)

. 16 )

 16 Before going any further, let us sump what shall be referred to as a flow in the low Mach regime for our two-phase model : such flow shall verify

Fig. 6 .

 6 Fig. 6.1 1-D sketch of a cell j across the interface Γ, separating domains D i , with its full neighborhood V and partial neighborhoods N ± ( j). Distance function φ , positive at cell j. Discontinuous pressure field p across the interface with non-dimensional values at cells j and k.

  Mach flow across the interface: behavior of the pressure near the interface during the acoustic step Let us consider π * ,θ,CP jk the non-dimensional form of corrected flux defined by (6.30), for k ∈ N -( j) ̸ = / 0 we have that π * ,θ,CP jk 1

1 )

 1 ) is generally maintained without truncation error of large magnitude such as O( 1 M and O( 1

Fig. 6 . 2 Fig. 6 . 3

 6263 Fig. 6.2 Magnitude of the velocity field of two-phase Gresho vortex problem #1 with a resolution of 100 × 100. (a): without low Mach correction; (b) AW type correction; (c) CP type correction; (d) new low Mach correction

Fig. 6 .

 6 Fig. 6.2 shows snapshots of Problem #1 at t = 0.1 for the numerical scheme without and with low Mach correction respectively. Without low Mach correction, the effect of dissipation is obvious and the vortex has completely disappeared. All of the presented low Mach corrections give the same numerical results for problem #1 since the properties of the two fluids are the same. With a low Mach correction, the vortex is not degraded, showing the efficiency of low Mach correction. Fig. 6.3 -Fig. 6.5 present numerical results of the Gresho test case with different fluid properties. From these figures, we can get a global conclusion: numerical scheme without low Mach fix can not preserve the vortex while the new low Mach correction always provides a satisfactory prediction. Fig. 6.3 presents numerical results of Problem #2 with large density ratio. For the numerical result related to the CP type correction, there are some instabilities as soon as t = 10 -3 . As presented in the asymptotic analysis, the CP type correction provides an error of magnitude of 1 α m for two fluids with the same sound speed. The acoustic impedance weighted correction can well preserve the vortex as well as new low Mach correction since the new low Mach fix gives the same truncation error on the momentum equation as the AW low Mach correction. The correction on u * does not seem to play an important role in this test case.

Fig. 6 . 4

 64 Fig. 6.4 Magnitude of the velocity field of two-phase Gresho vortex problem #3 with a resolution of 100 × 100. (a): without low Mach correction; (b) AW type correction; (c) CP type correction; (d) new low Mach correction.

Fig. 6 . 5

 65 Fig. 6.5 Magnitude of the velocity field of two-phase Gresho vortex problem #4 with a resolution of 100 × 100. (a): without low Mach correction; (b) AW type correction; (c) CP type correction; (d) new low Mach correction.

Fig. 6 .Fig. 6 . 6

 666 Fig. 6.4 shows the numerical results of Problem #3. All these low Mach corrections can well preserve the vortex. For two fluids of the same initial density and large sound speed ratio, all these low Mach corrections provide with errors of magnitude O(1), uniform with respect with M, α m and β m , but the asymptotic-preserving conditions are not strictly respected by AW or CP type correction.

1 0Fig. 6 . 7

 167 Fig. 6.7 Temporal evolution of dimensionless velocity fluctuations of a 2D bubble (La=12000) at resolution of R = 12.8∆x

Fig. 7 . 1

 71 Fig. 7.1 Presentation of Kokkos.

Fig. 7 . 2

 72 Fig. 7.2 Sketch of the sloshing tank with acceleration vectors.

Fig. 7 . 3

 73 Fig. 7.3 Temporal evolution of interface at x = 0 and at x = L with a uniform mesh 40 × 90. red line: present method with low Mach correction; blue line: present method without low Mach correction; black line: analytical solution.

Fig. 7 . 4

 74 Fig. 7.4 Temporal evolution of interface at x = 0 and at x = L with a uniform mesh 80 × 180. red line: present method with low Mach correction; blue line: present method without low Mach correction; black line: analytical solution.

Fig. 7 . 5 Fig. 7 . 6 Fig. 7 . 7

 757677 Fig. 7.5 Vorticity of Rayleigh-Taylor instabilities test case with a resolution of 320 × 640 at t * = 1 (left) and t * = 3 (right), black lines: the interface.

Fig. 7 . 8 Fig. 7 . 9

 7879 Fig. 7.8 Temporal evolution of velocity fluctuations of a 2D bubble (La=12000) at different resolutions: without low Mach correction.

Fig. 7 .

 7 Fig. 7.8 and Fig 7.[START_REF] Bouchut | Nonlinear stability of finite Volume Methods for hyperbolic conservation laws: And Well-Balanced schemes for sources[END_REF] show the spatial convergence of numerical scheme resolutions without and with low Mach correction respectively. From these two figures, we can find that the initial maximum amplitudes of spurious currents decrease as we refine mesh which shows a good convergence of the numerical scheme. Because of the total viscosity (including numerical viscosity and physical viscosity

9

 9 Fig. 7.8 and Fig 7.[START_REF] Bouchut | Nonlinear stability of finite Volume Methods for hyperbolic conservation laws: And Well-Balanced schemes for sources[END_REF] show the spatial convergence of numerical scheme resolutions without and with low Mach correction respectively. From these two figures, we can find that the initial maximum amplitudes of spurious currents decrease as we refine mesh which shows a good convergence of the numerical scheme. Because of the total viscosity (including numerical viscosity and physical viscosity

Fig. 7 .

 7 Fig. 7.10 Sketch of the domain for test case "rising bubble".

Fig. 7 .

 7 Fig. 7.11 Interface at different non-dimensional time t * = t R/g of present method (black line) and reference (blue circles)[96].

Fig. 7 .

 7 Fig. 7.13 Temporal evolution of left liquid interface position.

Fig. 7 .

 7 Fig. 7.14 Velocity profile at t = 1 ms in the tube.

Fig. 7 .

 7 Fig. 7.15 Temperature profile at t = 1 µs in the tube.

Fig. 7 .

 7 Fig. 7.16 Temperature profile at t = 4 ms in the tube.

Fig. 7 .

 7 Fig. 7.17 Condition of simulation at a given time for the Stefan problem

Fig. 7 .Fig. 7 . 6 Fig. 7 .

 7767 Fig. 7.18 Temporal evolution of the interface position.

Fig. 7 .Fig. 7 .

 77 Fig. 7.21 Condition of simulation at a given time for the Sucking problem

Fig. 7 . 5 Fig. 7 .

 757 Fig. 7.23 Temperature profile of the sucking problem at t = 0.6 s.

Fig. 7 .Fig. 7 .

 77 Fig. 7.26 Temperature distribution for ∆T = 4 K at t = t 0 (a) and t = 1.21 t 0 (b)

Fig. 7 .profile at t = 1 .21 t 0 Fig. 7 .

 7107 Fig. 7.28 Representation of velocity field at t = 1.21 t 0 , black lines: the interface.

Fig. 7 .Fig. 7 .Fig. 7 .

 777 Fig. 7.30 Numerical results for Ja = 18 (∆T = 6K)

Titre:

  Une méthode d'interface raide pour écoulement diphasique à faible nombre de Mach avec le changement de phase Mots clés: interface raide, écoulement diphasique compressible, correction à faible nombre de Mach, changement de phase

Table 4 .

 4 1 Level Set advection L 1 errors for test case #1 sin 4 (πx) by using the separate time-space approaches

		10e-2	2.91	6.41e-3	2.84	5.33e-3	2.74
	80 × 80	9.00e-4	3.61	7.14e-4	3.17	6.69e-4	2.99
	160 × 160 8.49e-5	3.41	8.56e-5	3.06	8.41e-5	2.99
	320 × 320 1.01e-5	3.07	1.05e-5	3.03	1.05e-5	3.00

Table 4 .

 4 2 Level Set advection L 1 errors for test case #1 sin

		.1.
		OS5		OS7
	Mesh	L 1 error L 1 order L 1 error L 1 order
	20 × 20	9.24e-3	-	2.91e-3	-
	40 × 40	3.44e-4	4.75	2.87e-5	6.66
	80 × 80	1.12e-5	4.94	2.39e-7	6.91
	160 × 160 3.53e-7	4.99	1.89e-9	6.9
	320 × 320 1.11e-8	4.99	1.48e-11	7.00

Table 4

 4

.3. From this table, we can find both of these two redistancing methods converge. According to the global errors, the Hopf-Lax based method is more accurate with errors of one order of magnitude smaller

Table 4 .

 4 -7 with a mesh of 256 × 256 and its local errors are always two orders magnitude smaller than the PDE based method. According this test case, in terms of both accuracy and efficiency, the Hopf-Lax based method without check procedure has a better performance than the PDE based method. 3 Average global and local L 2 errors for test case #1 Circle

		PDE based method	Hopf-Lax based method
	Mesh	L G 2	L L 2	t cpu	L G 2	L L 2	t cpu
	32 × 32	1.42e-3	1.32e-3	1.48e-2	1.18e-4	6.71e-5	6.85e-3
	64 × 64	3.84e-4	3.17e-4	1.16e-1	2.52e-5	7.91e-6	2.95e-2
	128 × 128 1.22e-4	1.49e-4	9.13e-1	5.69e-6	9.70e-7	1.42e-1
	256 × 256 6.01e-5	8.59e-5	7.57e0	1.33e-6	1.22e-7	7.85e-1

Table 4 .

 4 4. Comparing the errors and the CPU time, we can get the same conclusion as the test case "Circle" redistancing, in the aspect

		PDE based method	Hopf-Lax based method
	Mesh	L G 2	L L 2	t cpu	L G 2	L L 2	t cpu
	32 × 32	1.70e-3	1.03e-3	1.91e-2	2.65e-6	1.77e-6	3.66e-3
	64 × 64	7.97e-4	3.57e-4	1.25e-1	4.93e-7	9.28e-7	1.70e-2
	128 × 128 3.82e-4	1.25e-4	9.37e-1	1.46e-7	3.26e-7	8.88e-2
	256 × 256 1.87e-4	4.42e-5	6.50e0	4.58e-8	1.15e-7	5.48e-1

Table 4 .

 4 [START_REF] Ansari | Numerical simulation of compressible two-phase flow using a diffuse interface method[END_REF] Average global and local L 2 errors for test case #2 Square of both accuracy and efficiency, the Hopf-Lax based method has a better performance than the PDE based method.

Table 4 .

 4 5 Average global and local L 2 errors for test case #3 Ellipse

	-7 with a

Table 6 .

 6 Table. 6.1, we define Problem #1 -#4 with different fluid properties. 1 Numerical parameters of Gresho vortex test problems Problem #1 is original test case from

				Fluid 1			Fluid 2
		#	ρ 1	M 1,max γ 1	ρ 2	M 2,max γ 2	α m	β m
		1	1	10 -3	1.666	1	10 -3	1.666	1	1
	Gresho	2	1	10 -3	1.666	100	10 -3	7.14	0.01	1
	vortex	3	1	10 -3	1.666	1	10 -5	7.14	1	0.01
		4	1	10 -5	1.666	100	10 -3	7.14	0.01	0.01

1/2 , (4.54)

Acknowledgements

, thanks a lot to my parents, my sister, my girlfriend Yuxue Shang and all my friends for supporting me spiritually throughout my life. I would also like to express my gratitude to the China Scholarship Council (CSC), Maison de la Simulation, LIMSI and

Appendix A Numerical resolution of system with a barotropic EOS This work concentrates on low Mach two-phase flows, to assess the performance of our numerical approach on verification test cases of incompressible flows where fluid density variation is independent of thermal effects, we adapt a barotropic EOS: the pressure depends only on fluid density:

where p 0 , ρ 0 , c 0 are three constant physical quantities representing the reference pressure, density and sound speed. The fluid temperature is independent of the fluid density and pressure, equations of mass conservation and momentum balance are identical as (2.1a) and (2.1b), while the energy conservation equation can be simplified as:

2) is a simplification of (2.1c) by filtering the acoustic phenomenon. In this thesis, we will first concentrate on the numerical resolution of the global system equipped with a Mie-Gruneisen EOS, the numerical method to resolve (A.2) will be given latter. The system equipped with a barotropic EOS can be expressed as:

Numerical resolution of system with a barotropic EOS where p i = p EOS i (ρ i ). With the acoustic-transport splitting technique, the global system can be divided into: an acoustic subsystem

the transport subsystem:

the diffusion subsystem:

Numerical resolution of the acoustic subsystem

Subsystem (A.4) shares the same wave structure with W RP , thus the acoustic subsystem can be discretized as:

with the definition of u * , π * and {ρ∇Ψ} j given (5.17) and (5.18).

Numerical resolution of the transport subsystem

We write the transport subsystem with conservative variables which causes conservative terms to appear:

This system is approximated with an upwind scheme using the values of the normal velocities u * jk computed at each cell edge ∂ Ω jk during the acoustic step. The discretized equations then read

where b

′ n+ k,ghost is a ghost value that associated with the so-called ghost-cells Ω k when k ∈ N -( j) as represented in figure 5.2. Here, we want to mention that {ρ, (ρu), (ρv), (ρw)} ′ n+ k,ghost are define in the same way as in Section 5.2.2. While the temperature in ghost cells are approximated by a second-order extrapolation, more detailed about the temperature extrapolation method can be found in [START_REF] Abadie | On the combined effects of surface tension force calculation and interface advection on spurious currents within volume of fluid and level set frameworks[END_REF][START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF].

Numerical resolution of the diffusion subsystem

Terms related to diffusion effects in (A.6) are treated in the same way as the approximation of diffusion subsystem equipped with a Mie-Gruneisen EOS. For more details, please refer to Section 5.2.3.

Appendix B

Behavior of the pressure near the interface with AW type correction Let us consider π * ,θ,AW jk the non-dimensional form of corrected flux defined by (6.30), for k ∈ N -( j) ̸ = / 0 we have that π * ,θ,AW

By (6.32), (6.33), (6.34), (6.37), (6.38) and (6.29), tedious calculations show that π * ,θ,AW

Then we can see that when α 2 becomes large, it may produce important errors for the estimate of the pressure gradient in Ω j . Indeed, let us consider π * ,θ,AW

Using (6.33) we get π * ,θ,AW

Behavior of the pressure near the interface with AW type correction the error term of magnitude O(

) is issued from the first part of RHS of π * ,θ,AW jk in (6.29) and it may become significantly large. For example when α 2 β 2 ≤ O(1) then the error becomes of magnitude O(α 2 ).

Appendix C

Subcharacteristic conditions for density positivity

We follow the work presented in [START_REF] Bouchut | Nonlinear stability of finite Volume Methods for hyperbolic conservation laws: And Well-Balanced schemes for sources[END_REF]. With intermediate densities presented in (6.55), the acoustic impedances a l and a r should satisfy the following conditions:

then the approximate Riemann solver could preserve positiveness of densities ρ l , ρ r , ρ * ,θ ,N l , ρ * ,θ ,N r . With the second law of thermodynamics, we have

where s, c p and τ are the entropy, specific heat capacity at constant pressure and specific volume respectively. According to the equation of state, c p can be expressed as:

Combining (C.2) and (C.3), we can obtain that:

For a polytropic fluid, we have:

Subcharacteristic conditions for density positivity

where

. With (C.5), we can make the following assumptions:

An inverse function ζ can be defined that:

With assumptions in (C.6), we can derive the following relationship:

According to the monotonicity of ζ and (6.55), the conditions (C.1) can be expressed as:

ζ (a r ) .

(C.10)

These conditions (C.10) include the positivity of ρ * ,θ ,N l and ρ * ,θ ,N r and can be satisfied with the following acoustic impedances:

Proof:

Here we give the proof for the constraint on acoustic impedance a l . The definition of a l always satisfies the first condition at the first line of (C.10). To derive required constraint on the second condition, we define a variable X

If X ≤ 0, the second condition at the first line of (C.10) is always satisfied. Assume that X ≥ 0, then a l can be rewritten as: a l = ρ l ( p ′ (ρ l ) + ϒX). Multiplying by ρ l , the second condition at first line of (C.10) reads:

1 -1 2

By defining a variable ϖ that:

with (C.9), we can always have:

where 0 < ϖ ≤ 1 and ϒ ≥ 1. Thus the definition of a l can always satisfy (C.14), the positivity of ρ * ,θ ,N l is then guaranteed. With the same analysis, we can justify that constraints on a r presented in (C.11) and (C.12) can preserve the positivity of ρ * ,θ ,N r .