Mireille Haddad
email: elhaddad@ljll.math.upmc.fr

AND Toni Sayah
email: toni.sayah@usj.edu.lb.

INTERFACE TRANSPORT SCHEME OF A TWO-PHASE FLOW BY THE METHOD OF CHARACTERISTICS

Keywords: Two-phase flow, Level-Set function, finite element method, characteristics method, boundary condition under pressure effect, projection method

In this paper, we study an interface transport scheme of a two-phase flow of an incompressible viscous immiscible fluid. The problem is discretized by the characteristics method in time and finite elements method in space. The interface is captured by the Level-Set function. Appropriate boundary conditions for the problem of mould filling are investigated, a new natural boundary condition under pressure effect for the transport equation is proposed and an algorithm for computing the solution is presented. Finally, numerical experiments show and validate the effectiveness of the proposed scheme.

Introduction.

This paper is devoted to the study of an interface transport scheme separating two immiscible incompressible viscous fluids. This problem involves a wide range of real-life physical phenomena having a major importance in several industrial applications; within which we are especially interested in modeling mould filling in iron foundry. In many such fluid flows, the physical time scale and length scales are so small that credible experiments are very expensive. Thus numerical analysis appears to be the only way to understand and solve the problem. The development of a reliable computational strategy of such problems requires the accurate discretization and tracking of the free surface.

There are, however, major challenges in the context of multiphase flow modeling. First, we have to take into account the evolution of the interface and its topological changes. Second, we have to deal with the non-linearity for the convection of the flow and the interface. Third, we must assign appropriate boundary conditions to the flow and transport equations. In addition, care must be taken in treating the geometrical and topological singularities across the interface. We also have to maintain a sharp interface resolution, including the cases of interface folding, breaking and merging. Finally, we should respect the physical properties such as the mass conservation for any incompressible fluid flows.

Since the seminal work of Harlow and Welch [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF], several methods have been developed to solve the problem of interfacial flows numerically, we cite these studies [START_REF] Anderson | Diffuse-interface methods in fluid mechanics[END_REF], [START_REF] Cuvelier | Some numerical methods for the computation of capillary free boundaries governed by the Navier-Stokes equations[END_REF], [START_REF] Shyy | Computational Fluid Dynamics with moving Boundaries[END_REF] and the references therein as examples. The most popular way is to divide these methods into two main categories according to the mesh:

(1) In the Lagrangian methods, a set of equidistant markers is used to track the interface motion.

Each computational cell carries always the same fluid portion and the mesh moves with the fluid. As time evolves, marker points have to be relocated along the new interface. Furthermore, this method requires the transfer of information between the interface and the fixed mesh once the interface has been moved. In addition, in this approach it is hard to deal with the evolution of interface markers when the interface becomes severely stretched or deformed. More details about these methods and their implementation can be found in the references [START_REF] Unverdi | A front-tracking method for viscous incompressible multi-fluid flows[END_REF], [START_REF] Esmaeeli | Direct numerical simulations of bubbly flows. Part 1. Low Reynolds number arrays[END_REF], [START_REF] Esmaeeli | Direct numerical simulations of bubbly flows. Part 2. Moderate Reynolds number arrays[END_REF], [START_REF] Hirt | An arbitrary Lagrangian-Eulerian computing method for all flow speeds[END_REF] and [START_REF] Donea | Finite Element Methods for Flow Problem[END_REF]. To overcome the drawbacks of the Lagrangian methods, researchers attempted hybrid approaches; one of the most common method proposed was the arbitrary Lagrangian-Eulerian (ALE) method [START_REF] Torres | The point-set method: front-tracking without connectivity[END_REF], [START_REF] Shin | Modeling three-dimensional multiphase flow using a level contour reconstruction method for front tracking without connectivity[END_REF], [START_REF] Baensch | Finite element discretization of the Navier-Stokes equations with a free capillary surface[END_REF], [START_REF] Behr | Stabilized space-time finite element formulations for free-surface flows[END_REF] and [START_REF] Ganesan | A coupled arbitrary Lagrangian-Eulerian and Lagrangian method for computation of free surface flows with insoluble surfactants[END_REF]. They consist in solving the problem using a mesh and then this is moved according to the flow velocity field. The ALE schemes lead to satisfying results but the intricacy of implementation seems to be expensive, especially in three-dimensional problems.

(2) The Eulerian methods, on the other hand, are characterized by a fixed mesh and a traveling fluid from one grid cell to another. The most commonly used approaches among the Eulerian methods are those of the interface capturing methods. The location of the interface is determined by the advection of either a characteristic function or a zero-isocontour of a continuous function by the fluid velocity. The obtained transport equation expresses that the interface is a material line propagating with the fluid. The most popular Eulerian methods approaches are:

• The discontinuous approach was introduced by Hirt and Nichols [START_REF] Hirt | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF] by the so-called volumeof-fluid (VOF) method. This is the first interface capturing method using the concept of transport of a scalar field. The main idea of the VOF method is to use a scalar field to locate the two fluids. It is a very popular method for modeling free boundaries in hydrodynamics problems. The advantages of this tool are numerous: it can manage topology changes in the interface such as breaks and reconnections, it naturally conserves mass, and it can easily be extended to 3D space (see [START_REF] Hirt | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF] and [START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF]). There are however major drawbacks for this method, one is the necessity to advect a discontinuous function, which requires a specific numerical treatment of the transport equation, another is represented by the difficulties in determining the precise location of the interface as well as its geometric characteristics. Moreover, even if some reconstruction algorithms are effective in improving the precision, they are complex and expensive to implement in the three dimensional space. Finally, the quality of this method will depend on both the reconstruction method of interface and of the numerical scheme for solving the transport equation. • The first algorithm on the continuous interface capturing approach was suggested by Dervieux and Thomasset [START_REF] Dervieux | A finite element method for the simulation of a Rayleigh-Taylor instablity[END_REF]. Later, the development of this approach evolved based on the Level-Set method, we may cite [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF], [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flows[END_REF], [START_REF] Chang | A level set formulation of Eulerian interface capturing methods for incompressible fluid flows[END_REF], [START_REF] Zhang | A curvilinear level set formulation for highly deformable free surface problems with application to solidication[END_REF], [START_REF] Zhao | Capturing the behavior of bubbles and drops using the variational level set approach[END_REF] and [START_REF] Sussman | An adaptive level set approach for incompressible two-phase flows[END_REF]. The concept of this method is to define a regular scalar function through the interface, a distance function, for which the zero level set is the interface that we are looking to describe. Solving the transport equation by advecting this distance function makes it possible to predict the evolution of the interface propagating with the fluid velocity field. The Level-Set method takes into account the topological changes naturally. This method considerably simplifies the interface convection problem, by transforming the discontinuous function into a continuous one. It also makes it easier to compute the geometric characteristics from the distance function and creates a smoother extension to three-dimensional problems. A number of disadvantages can however challenge the performance of the method. During the computing process after several time iterations of the computing process, the distance property and the mass are not preserved. A re-initialization algorithm is usually used to keep the Level-Set function a signed distance function; in addition a mass conservation process is required in order to respect the physical properties.

Based on this state of art, we are motivated to select the context of the Level-Set method. Also, we decouple the purely convective part of the flow and of the transport equation from the advection part; the treatment of the nonlinear convection term is thus reduced to a problem of research of characteristic curves. Furthermore, we select the finite element scheme for the discretization in space and we maintain the convenience and effectiveness of Eulerian grid. In addition, we choose the Galerkin finite element for the discretization in space as it assumes the minimal regularity for the existence and the uniqueness of a solution.

The outline of the paper is as follows:

• In section 2, we present the problem of mould filling.

• In section 3, we analyze the corresponding continuous problem.

• In section 4, we introduce the discrete problem and recall its main properties.

• In section 5, we study the problem by the projection method.

• In section 6, we correct the Level set by re-initialisation and mass conservation.

• In section 7, we show numerical results of validation for gas-water modelling.

• In section 8, we show numerical results of validation for gas-liquid metal modelling.

Problem of mould filling

To formulate the appropriate mathematical model to our problem, we start by the physical hypothesis that we have assumed. We consider an unsteady and laminar flow of two immiscible fluids. In this flow, the two fluids are supposed to be viscous, Newtonian and of large density ratio. In addition, the fluids are considered incompressible and isothermal, thus neglecting the variations of density and viscosity due to changes in pressure and temperature. Furthermore, by assuming that both fluids are homogenous, we believe that the viscosity and density are constants in each fluid. The two fluids are immiscible and the separation zone between the fluids is a sharp interface of zero thickness wherein the physical properties of the two fluids change abruptly. To treat the variations of the physical properties across the interface, it is necessary to implement jump conditions. In our study, we neglect the surface tension between the two fluids. We assume the interface is impermeable, thus the mass transfer across the interface is neglected (See [START_REF] Batchelor | An Introduction to Fluid Dynamics[END_REF] and [START_REF] Lamb | Hydrodynamics[END_REF]).

We consider an interval [0, T ] ⊂ IR, where T is a positive real number, and an arbitrary time t ∈ [0, T ].

Let Ω be a bounded simply connected open domain in IR d , d = 2, 3 , with a Lipschitz-continuous connected boundary ∂Ω. We denote by n the outward unit normal vector to the interface ∂Ω and (e 1 , e 2 ) the canonical base of R 2 (respectively (e 1 , e 2 , e 3 ) the canonical base of R 3 ).

We suppose that Ω represents a mould containing two fluids, thus, at each time t ∈ [0, T ], it is divided into two open sub-domains Ω 1 and Ω 2 evolving in time and separated by the interface Γ such that

Ω = Ω 2 (t) ∪ Ω 1 (t) and Ω 1 (t) ∩ Ω 2 (t) = ∅.
We denote by ∂Ω i the boundary Ω i , i = 1, 2, which is divided into four parts such that We denote by u i and p i the velocity and the pressure in the domain Ω i , i = 1, 2 respectively, and by µ i and ρ i the constant dynamic viscosities and densities of the fluid in Ω i respectively. The bi-fluid flow motion is described in each subdomain and at each time t ∈]0, T [ by the following incompressible Navier-Stokes equations.

∂Ω 1 ∩ ∂Ω 2 = Γ(t), ∂Ω = Γ 0 ∪ Γ 1 ∪ Γ 2
       ρ i ∂u i ∂t + (u i .∇)u i -div(2µ i Du i ) + ∇p i = f i in Ω i (i = 1, 2), div u i = 0 in Ω i (i = 1, 2), (2.1)
where

Du i = 1 2 (∇u i + t ∇u i )
is the deformation rate tensor and f i represents a density of body forces in

Ω i , i = 1, 2.
The density and the viscosity of the fluid can be written in Ω as

ρ(x, t) = ρ 1 1 x∈Ω1(t) + ρ 2 1 x∈Ω2(t) (respectively µ(x, t) = µ 1 1 x∈Ω1(t) + µ 2 1 x∈Ω2(t) ),
where 1 x∈Ωi is the characteristic function of the subdomain Ω i , i = 1, 2. We denote by f the data, u the velocity and p the pressure of the fluid in Ω such that f = f i , u = u i and p = p i in Ω i , i = 1, 2. The problem (2.1) can be rewritten as: at each time t ∈]0, T [

       ρ ∂u ∂t + (u.∇)u -div(2µ Du) + ∇p = f in Ω 1 ∪ Ω 2 , div u = 0 in Ω 1 ∪ Ω 2 .
(

This system must be endowed with adequate boundary conditions and initial conditions, thus, we will consider on the inlet Γ 0 a non-homogeneous boundary condition of Dirichlet type: u = U in , on the free surface Γ 2 a do-nothing boundary condition: (2µDu -pI).n = 0, and on the wall Γ 1 we will compare two different boundary conditions:

(BCU)       
Dirichlet boundary conditionss: u = 0 or Navier boundary conditions:

u.n = 0 and αu.τ + t n(2µDu -pI).τ = 0, where τ is the tangential unit vector and α is the friction coefficient. In fact, many researchers consider boundary conditions of Dirichlet type. However, as noted by Serrin [START_REF] Serrin | Mathematical principles of classical fluid mechanics[END_REF], they are not always realistic and in general lead to boundary layers phenomena next to the walls (as we will see in section 7, Numerical Results). Navier [START_REF] Navier | Sur les lois de l'équilibre et du mouvement des corps élastiques[END_REF] has proposed a so-called slip boundary provided with friction, at the wall that allows taking into account the slip of fluid next to the boundaries and measuring the friction effect. Besides, the system (2.2) is also completed with interface conditions imposing the continuity of the velocity and the balance of the normal stress with the surface tension across the interface Γ(t), namely: In this paper, we neglect the surface tension.

We intend to work with the following problem: at each time t ∈]0, T [,

                                                       ρ ∂u ∂t + (u.∇)u -div(2µDu) + ∇p = f in Ω 1 ∪ Ω 2 , div u = 0 in Ω, u = U in on Γ 0 , (BCU) on Γ 1 , (2µDu -pI).n = 0 on Γ 2 , [2µDu -pI].n = 0 on Γ, [u] = 0 on Γ, u(x, 0) = 0 in Ω, (2.3)
where f = -ρge d is the gravitational force vector, I is the identity matrix and U in designates the velocity of the flux at the inlet.

For the interface transport, the main challenge is to handle geometrical and topological changes. Thus we solve the problem using the Level-Set function on a fixed uniform mesh. In particular, we follow here Osher and Sethian (see [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF]), we introduce the signed distance function to the interface Γ(t):

φ(x, t) = ± min y∈Γ(t) |x -y|, (2.4)
where the function φ is set to be negative in the domain Ω 1 (t) and positive in the domain Ω 2 (t). Hence, at each time step, the fluid interface corresponds to the zero isocontours of the continuous function φ:

Γ = {x ∈ R d /φ(x, t) = 0, ∀t ≥ 0}.
The density ρ and the viscosity µ can be rewritten in Ω as:

ρ(x, t) = ρ 1 (φ ≤ 0) + ρ 2 (φ > 0) and µ(x, t) = µ 1 (φ ≤ 0) + µ 2 (φ > 0).
The interface is then captured, at each time step, by the advection of the Level-Set function by the fluid velocity. It can be described by the following transport equation:

∂φ ∂t + u.∇φ = 0. (2.5)
After a very small time, discontinuities appear over the interface next to the boundaries because there is no uniqueness of the solution for a general continuous velocity field u in this strong formulation. To avoid them we may find the solution in the sense of viscosity. This method was introduced by P.L.Lions et M.G. Crandall [START_REF] Crandall | Viscosity solutions of Hamilton-Jacobi equations[END_REF] and selects the weak discontinuous physically significant solution by adding an artificial viscosity -ε∆φ that vanishes as ε → 0 [START_REF] Barles | Solutions de viscosité et équations elliptiques du deuxième ordre[END_REF].

We denote by h the grid size of the mesh, we choose the parameter ε to be proportional to h as in [START_REF] Olsson | A conservative level set method for two phase flow[END_REF].

The choice of ε is very delicate, a small ε gives better conservation of the area (volume) bounded by the zero contour of φ since the volume error increases proportionally to ε. There are however numerical restrictions on how small we can choose this parameter [START_REF] Olsson | A conservative level set method for two phase flow II[END_REF]. In our application, numerical tests show that below a critical value of epsilon (for ε ≤ h/50) discontinuities occur over the interface, for h/50 ≤ ε < h/3 oscillations occur over the interface and for ε ≥ h/3 we obtain neither discontinuities nor oscillations but the volume error is remarkably large (as we will show in the numerical results).

Appropriate boundary and initial conditions must also be assigned to φ. There exists in the literature variant of boundary conditions that can be assigned to φ, for example Dirichlet boundary conditions, homogeneous Neumann boundary conditions and contact line boundary conditions, which depends on the wettability property of a solid surface by a liquid via the Young equation [START_REF] Qian | Molecular Hydrodynamics of the Moving Contact Line in Two-phase Immiscible Flows[END_REF].

In this work, we will consider homogeneous Neumann boundary conditions, which are used in a lot of applications (See for example [START_REF] Olsson | A conservative level set method for two phase flow[END_REF]) and we will propose a new natural boundary condition, which will be called non-homogeneous boundary condition under pressure effect. This last one calculates the slope of the angle between the free surface and the wall boundaries in a natural way that makes sense from a physical point of view. The non-homogeneous boundary condition on Γ 1 under pressure effect is a boundary condition on the pressure:

We multiply the first equation of (2.3) by n.

ρ∂ t u.n| Γ1 + ρ (u.∇)u .n| Γ1 -div(2µDu).n| Γ1 + ∇p.n| Γ1 = -ρge d .n| Γ1 (2.6)
The first term (∂ t u).n = ∂ t (u.n) = 0 since u.n = 0.

And we get the following boundary condition on Γ 1 :

∂ n p| Γ1 = -ρge d .n| Γ1 -ρ (u.∇)u .n| Γ1 -2µDu.n| Γ1 .
In our application, the viscosity is a very small number. Thus we neglect the second term in the right hand side of the last equation, and therefore the boundary condition on the pressure on Γ 1 can be written as:

∂ n p| Γ1 = -ρge d .n| Γ1 -ρ (u.∇)u .n| Γ1 .
We neglect the pressure of gravity in Ω 2 that contains a fluid with very small viscosity. We impose on Γ 0 an average velocity U in since generally we know the quantity of fluid that enters during the experiment. Otherwise, we impose on the interface Γ a do nothing boundary condition : [2µDu -pI].n = 0. Since the viscosity is a very small number, we will assume that the normal forces acting on the interface are negligible, namely 2µDu.n is almost 0 on the interface (See figure [START_REF] Esmaeeli | Direct numerical simulations of bubbly flows. Part 1. Low Reynolds number arrays[END_REF] and [START_REF] Hirt | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF] in the numerical results for a comparison between the forces acting on the interface and the forces exerted on the entry that validate the hypothesis in 2D and 3D problems). As a result, the total pressure vanishes on the interface Γ, which is the zero isocontour. Hence the pressure and the Level-Set function have the same isovalues φ = p = 0 on the free surface that implies that φ can be considered as equal to p up to a multiplicative function c (φ = -cp) in a neighborhood of the boundary Γ.

The non-homogeneous boundary condition under pressure effect can be written on a neighberhood of Γ ∩ Γ 1 as:

∂ n φ = -c ∂ n p -p ∂ n c = -c∂ n p ( since p is almost 0). (2.7)
Also, since ∇φ = 1, we normalize the boundary condition and we get:

∂ n φ = ρge d + ρ (u.∇)u .n ρge d + ρ (u.∇)u .n . (2.8)
We denote by G(u) the right hand side of (2.8). We impose two different boundary conditions on the transport equation on (2.9) We will show in the numerical results that these boundary conditions under pressure effects give the physical slope of the angle between the free surface and the interface as well as decrease the volume error remarkably. The equation (2.5) endowed with boundary conditions will be written in the following form:

       ∂ t φ + u.∇φ -ε∆φ = 0 in Ω, (BC φ) on Γ 1 , ∂ n φ = 0 on Γ 0 ∪ Γ 2 , φ(x, 0) = φ 0 in Ω, (2.10)
where φ 0 is the initial position of the interface.

Our system of equations will be (2.3) and (2.10). Well-posedness results for the general weak formulation of the Navier-Stokes problem for two-phase flows including the interface jump condition have been analyzed only for special cases. The case of a bounded domain Ω for arbitrary time intervals [0, T ], T > 0 was treated in [START_REF] Tanaka | Global existence of two-phase non-homogeneous viscous incompressible fluid flow[END_REF]; it provided a well-posedness result for the Navier-Stokes problem in a weak formulation.

Description of the algorithm

In order to write the variational formulation of the previous problem, we introduce the following Sobolev spaces (m and p ∈ N):

W m,p (Ω) = {v ∈ L p (Ω), ∂ α v ∈ L p (Ω), ∀|α| ≤ m}, H m (Ω) = W m,2 (Ω),
equipped with the following semi-norm and norm:

|v| m,p,Ω = { |α|=m Ω |∂ α v(x)| p dx} 1 p and v m,p,Ω = { K≤m |v| p K,p,Ω } 1 p .
We denote by X u the sub-space of H 1 (Ω) defined by:

X u = {u ∈ H 1 (Ω) / u = U in on Γ 0 and u = 0 on Γ 1 } : in the case of Dirichlet , {u ∈ H 1 (Ω) / u = U in on Γ 0 } : in the case of Navier . (3.1)
X v the sub-space of H 1 (Ω) defined by:

X v = {v ∈ H 1 (Ω) / v = 0 on Γ 0 ∪ Γ 1 } : In the case of Dirichlet , {v ∈ H 1 (Ω) / v = 0 on Γ 0 } : In the case of Navier . (3.2)
We denote by

M = L 2 (Ω), Y = H 1 (Ω) and we suppose that f ∈ L 2 (Ω) d .
The weak formulation of the problem with the Dirichlet boundary conditions can be written as:

Find (u, p) ∈ X u × M , φ ∈ Y such that:            ∂φ ∂t + u.∇φ, r +ε(∇φ, ∇r) -ε Γ1 G(u)r = 0 ∀r ∈ Y, ρ ∂u ∂t + ρ(u.∇)u, v +2(µDu, Dv) -(p, div v) = -(ρge d , v) ∀v ∈ X v , (q, div u) = 0 ∀q ∈ M. (3.3)
Let us now write the weak formulation with the Navier boundary conditions (BCU). To implement them, we couple the velocity componants in one equation and we use the penalty method as follows:

β -1 u.n + t n(2µDu -pI).n = 0 and αu.τ + t n(2µDu -pI).τ = 0. (3.4)
where β is a penalty coefficient, which is a small number [START_REF] Babuska | The Finite Element Method with Penalty[END_REF].

Thanks to the condition of incompressibility, it follows the relation ∆u = 2 div Du. The symmetry of the deformation tensor yields:

(Du, Dv) = (Du, ∇v 2 ) + (Du, t ∇v 2 ) = (Du, ∇v 2 ) + ( t Du, t ∇v 2 )
= (Du, ∇v)

and

div 2Du = div ∇u + div( t ∇u) = div ∇u + ∇(div u) = div(∇u).
Than the weak formulation of the problem with Navier boundary conditions can then be written as:

Find (u, p) ∈ X u × M , φ ∈ Y such that:                    ∂φ ∂t + u.∇φ, r +ε(∇φ, ∇r) -ε Γ1 G(u)r = 0 ∀r ∈ Y, ρ ∂u ∂t + ρ(u.∇)u, v +(µ∇u, ∇v) -(p, div v) - Γ1 2µ Du -pI n .v ds = -(ρge d , v) ∀v ∈ X v , (q, div u) = 0 ∀q ∈ M. (3.5)
The integral over the boundary Γ 1 in (3.5) can be rewritten by decomposing the test function v in the following way:

v = (v.n)n + (v.τ )τ .
That implies by using the definition of Navier boundary conditions(3.4):

Γ1 (2µDu -pI)n .v ds= Γ1 t n(2µDu -pI)n v .n ds + Γ1 t n(2µDu -pI)τ v .τ ds = Γ1 β -1 (u.n)(v.n)ds + Γ1 α(u.τ )(v.τ )ds = Γ1 β -1 t u(n t n)vds + Γ1 α t u(τ t τ )vds, (3.6)
Then, the variational formulation can be written as:

Find (u, p) ∈ X u × M , φ ∈ Y such that for all (v, q) ∈ X v × M and ∀r ∈ Y                    ∂φ ∂t + u.∇φ, r +ε(∇φ, ∇r) -ε Γ1 G(u)r = 0. ρ ∂u ∂t + ρ(u.∇)u, v +(µ∇u, ∇v) -(p, div v) - Γ1 β -1 t u(n t n)vds - Γ1 α t u(τ t τ )vds = (-ρge d , v).
(div u, q) = 0.

(3.7)

The discrete problem

In this section, we present the numerical strategy we have designed to resolve the continuous coupled system obtained in the previous section. Our approach is based on the characteristics method combined with a finite element method.

Discretization in time.

We propose a time discretization of (3.5) and (3.7 ) by the method of characteristics. This method, also known as the Lagrange-Garlekin method was introduced by Benqué [START_REF] Benqué | A finite element method for the Navier-Stokes equations[END_REF] and analyzed in [START_REF] Pironneau | On the transport-diffusion algorithm and its applications to the Navier-Stokes equations[END_REF]. The main idea behind this method is that the convection operator (the non linear term) can be turned into a total derivative by using a Lagrangian formulation. Thus, the treatment of the nonlinear convection term is reduced to a problem of searching the characteristic foot X(x; s; t), i.e the position of the particle at the previous time. This approach allows us to avoid theoretically the constraint CFL (Courant-Friedrichs-Levy) on the time step and it has been shown that it has very good stability properties, we may cite [START_REF] Pironneau | On the transport-diffusion algorithm and its applications to the Navier-Stokes equations[END_REF], [START_REF] Morton | Stability of the Lagrange-Galerkin method with non-exact integration[END_REF] and [START_REF] Fourestey | Stabilité des méthodes de Lagrange-Galerkin du premier et du second ordre[END_REF]. Furthermore, only the right-hand side has to be updated at each iteration during the resolution. Thanks to this formulation, it is theoretically possible to follow the particles over time along their trajectory by solving, for each particle, an ordinary differential equation called characteristics equation:

   ∂X ∂t (x, s; t) = u(X(x, s; t), t), X(x, s; s) = x, (4.1)
where the characteristics curve X(x, s; t) denotes the position at time t of a fluid particle located at position x at the time s.

We introduce a partition of the interval [0, T ] into N subintervals [t n , t n+1 ], such that ∆t = T N , the points t n = n∆t, for n = 0, ...., N , and denote by

u n (x) = u(x, t n ), p n (x) = p(x, t n ),ρ n (x) = ρ(x, t n ), µ n (x) = µ(x, t n ) and φ n (x) = φ(x, t n ).
Using the following approximation of the total derivative along the characteristic curves, we approximate Du Dt at the time t = t n+1 by:

Du Dt n+1 ≈ u(x, t n+1 ) -u(X n (x), t n ) ∆t , (4.2)
where X n (x) is the approximation of X(x, t n+1 : t n ).

Same for

Dφ Dt we approximate it at the time t = t n+1 by:

Dφ Dt n+1 ≈ φ(x, t n+1 ) -φ(X n (x), t n ) ∆t . (4.3)
Then along the characteristic curves, the variational formulation with Dirichlet boundary conditionss becomes:

               φ n+1 -φ n oX n ∆t , r +(ε∇φ n+1 , ∇r) -ε Γ1 G(u n )r = 0 ∀r ∈ Y, ρ n u n+1 -u n oX n ∆t , v + (µ n ∇u n+1 , ∇v) -(p n+1 , div v) = (-ρ n ge d , v) ∀v ∈ X v , (q, div u n+1 ) = 0 ∀q ∈ M, (4.4)
and the variational formulation with Navier boundary conditions becomes: We introduce the discrete spaces X The velocity is discretized with the Mini-Element:

                       φ n+1 -φ n oX n ∆t , r +(ε∇φ n+1 , ∇r) -ε Γ1 G(u n )r = 0 ∀r ∈ Y, ρ n u n+1 -u n oX n ∆t , v + (µ n ∇u n+1 , ∇v) -(p n+1 , div v) - Γ1 α t u n+1 (τ t τ )vds - Γ1 β -1 t u n+1 (n t n)vds = (-ρ n ge d , v) ∀v ∈ X v , (q, div u n+1 ) = 0 ∀q ∈ M.
u,h ⊂ X u , X v,h ⊂ X v , M h ⊂ M , Y h ⊂ Y
X u,h = {u h ∈ X u ; ∀k ∈ τ h ; u h|k ∈ P b (k) d } and X v,h = {v h ∈ X v ; ∀k ∈ τ h ; v h|k ∈ P b (k) d },
where the space P b (k) is spanned by functions in P 1 (k) and the bubble function on k (for each element k, the bubble function is equal to the product of the barycentric coordinates associated with the vertices of k).

The pressure is discretized with classical continuous finite element of order one:

M h = {q h ∈ M ∩ C 0 (Ω); ∀k ∈ τ h , q h|k ∈ P 1 (k)}.
The Level-Set function is also discretized with classical continuous finite element of order one:

Y h = {r h ∈ Y ∩ C 0 (Ω); ∀k ∈ τ h , r h|k ∈ P 1 (k)}.
The discrete system corresponding to the variational formulation with Dirichlet boundary conditionss can be written in the following form:

Find (u n+1 h , p n+1 h ) ∈ X u,h × M h and φ h ∈ Y h such that                φ n+1 h -φ n h oX n ∆t , r h +(ε∇φ n+1 h , ∇r h ) -ε Γ1 G(u n h )r h = 0 ∀r h ∈ Y h , ρ n h u n+1 h -u n h oX n ∆t , v h + (µ n h ∇u n+1 h , ∇v h ) -(p n+1 h , div v h ) = (f n , v h ) ∀v h ∈ X v,h , (q h , div u n+1 h ) = 0 ∀q h ∈ M h , (4.6)
where ρ n h et µ n h are the corresponding discrete densities and viscosities.

The discrete system corresponding to the variational formulation with Navier boundary conditions can be written in the following form:

Find (u n+1 h , p n+1 h ) ∈ X u,h × M h and φ h ∈ Y h such that                        φ n+1 h -φ n h oX n ∆t , r h +(ε∇φ n+1 h , ∇r h ) -ε Γ1 G(u n h )r h = 0 ∀r h ∈ Y h , ρ n h u n+1 h -u n h oX n ∆t , v h + (µ n h ∇u n+1 h , ∇v h ) -(p n+1 h , div v h ) - Γ1 α t u n+1 h (τ t τ )v h ds - Γ1 β -1 t u n+1 h (n t n)v h ds = (-ρ n h ge d , v h ) ∀v ∈ X v , (q h , div u n+1 h ) = 0 ∀q h ∈ M h . (4.7)
In the following, we call the schemes (4.6) and (4.7) by the "classical method" for corresponding Dirichlet and Navier boundary conditions.

Projection method

In this section, in order to reduce the CPU time and the used memory for the simulation of the problem, we use the projection method to solve the Navier-Stokes problem. This method was introduced by Chorin [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF], [START_REF] Chorin | On the convergence of discrete approximations to the Navier-Stokes equations[END_REF] and Temam [START_REF] Temam | Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires (I)[END_REF], [START_REF] Temam | Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires (II)[END_REF]. The problem is discretized by the characteristics method in time and the pair P b -P 1 of finite elements in space.

The algorithm of the projection method is based on the decomposition of the velocity vector field into a vector of divergence free and another irrotational. Typically, the algorithm is decomposed at each time step into three steps: the first step computes an intermediate velocity that does not satisfy the incompressibility condition; the second step projects this intermediate velocity on the set of divergence free functions to get the value of the pressure solution of the problem; The third step updates the final velocity from the obtained results. The algorithm is summarized as follows : We start with u 0 h = 0 and p 0 h = 0. Given u n , find (u n+1 , p n+1 ) such that :

Step1-Computation of the intermediate velocity u * :

                         ρ n u * -u n oX n ∆t -div(µ n ∇u * ) + ∇p n = ρ n ge d in Ω, [ µ n ∂u * ∂n ] = 0 in Γ, u * = U in on Γ 0 , (BCU) Dirichlet or Navier on Γ 1 , µ n ∂u * ∂n = 0 on Γ 2 .
(5.1)

Step 2-Computation of the pressure p n+1 :

                 div u * -div ∆t ρ n ∇(p n+1 -p n ) = 0 in Ω, [ p n+1 -p n ] = 0 on Γ, ∂ ∂n (p n+1 -p n ) = 0 on Γ 0 ∪ Γ 1 , p n+1 -p n = 0 on Γ 2 .
(5.2)

Step 3-Computation of the final velocity u n+1 :

• For Dirichlet boundary conditions:

u n+1 = u * - ∆t ρ n ∇(p n+1 -p n ) in Ω.
(5.3)

• For Navier boundary conditions, we solve the following problem:

                   -ε ∆u n+1 + u n+1 = u * - ∆t ρ n ∇(p n+1 -p n ) in Ω. u n+1 = U in on Γ 0 , 1 β u n+1 .n + t n(ε ∇u n+1 ).τ = 0 on Γ 1 , ∂u ∂n n+1 = 0 on Γ 2 , (5.4)
where ε is the parameter of penalization, which is a small number.

The computed pressure is in M = {q ∈ H 1 (Ω)/q| Γ2 = 0} and not in L 2 (Ω) and the final velocity belongs to the space H 0 div (Ω) = {v ∈ L 2 (Ω)/ div v = 0 in Ω} in the case of Dirichlet boundary conditionss and in H 1 (Ω) in the case of the Navier boundary conditions.(see Bell and Marcus [START_REF] Bell | A second-order projection method for variable-density flows[END_REF]). For the discretization in space, we introduce the following discrete spaces:

X uh = {u h ∈ X u ; ∀k ∈ τ h ; u h|k ∈ P 1 (k) d }, X vh = {v h ∈ X v ; ∀k ∈ τ h ; v h|k ∈ P 1 (k) d }, M h = {q h ∈ M ∩ C 0 (Ω); ∀k ∈ τ h , q h|k ∈ P 1 (k)} L h = {u h ∈ L 2 (Ω); ∀k ∈ τ h ; u h|k ∈ P 1 (k) d }.
The discrete variational formulation can be written in the following form:

1-Find u * h ∈ X uh such that for all v h ∈ X vh , we have:

• For Dirichlet boundary conditions:

(ρ n h u * h -u n h oX n ∆t , v h ) + (µ n h ∇u * h ∇, v h ) + (∇p n h , v h ) = Ω ρ n h ge d .v h .
(5.5)

• For Navier boundary conditions:

(ρ n h u * h -u n h oX n ∆t , v h ) + (µ n h ∇u * h ∇, v h ) + (∇p n h , v h ) - Γ1 α t u n+1 h (τ t τ )v h ds - Γ1 β -1 t u n+1 h (n t n)v h ds = Ω ρ h n ge d .v h .
(5.6)

2-Find p n+1 h

∈ M h such that for all q h ∈ M h , we have:

(div u * h , q h ) + ( ∆t ρ n h ∇(p n+1 h -p n h ), ∇q h ) = 0.
(5.7)

3-Find u n+1 h such that,
• For Dirichlet boundary conditions: Find u n+1 h ∈ L h such that for all v h ∈ L h , we have:

(u n+1 h , v h ) = (u * h , v h ) -( ∆t ρ n h ∇(p n+1 h -p n h ), v h ).
(5.8)

• For Navier boundary conditions: Find u n+1 h ∈ X uh such that for all v h ∈ X vh , we have:

(ε ∇u n+1 h , ∇v h ) - Γ1 β -1 t u n+1 h (n t n)v h ds + (u n+1 h , v h ) = (u * h , v h ) - ∆t ρ n h ∇(p n+1 h -p n h ), v h .
(5.9)

Level set correction

It is well known that numerous errors affect the numerical algorithm and perturb the mass conservation as time evolves in two phase flow modeling. In this section, we introduce several corrections and ameliorations of the algorithm in order to get satisfying results.

6.1. Algorithm of reinitialization. At the initial time, all levels lines are calculated using the definition of the Level-Set method. As time evolves, the advection of the Level-Set by a velocity field causes the contour lines to become very tight (steep Level-Set) in some areas and spaced (flat Level-Set) in others, thus the method becomes imprecise and algebraic distance property ∇φ = 1 is lost. This may cause numerical errors that affect interface shape, its geometric characteristics and, moreover, the mass conservation.

To overcome this, we use the algorithm of reinitialization proposed by Sussman et al. [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flows[END_REF], which is based on the following equation

∂Φ ∂τ = sign(φ)(1 -∇Φ ) in Ω × (0, τ ), (6.1) with Φ(x, t, τ = 0) = φ(x, t), (6.2)
where τ is an imaginary time. We solve (6.1) iteratively until it reaches a steady state, we obtain the distance property ||∇Φ|| = 1.

In order to discreticize the equation (6.1), we rewrite it in the following form:

∂Φ ∂τ + w∇Φ = sign(φ) with w = sign(φ) ∇Φ ∇Φ , (6.3)
where ∆τ is the time step corresponding to the imaginary time τ . By using the characteristics method, the discrete variational formulation can be written as:

Find Φ n+1 ∈ Y such that ∀r ∈ Y , Φ n+1 -Φ n oX n (w n ,∆τ ) ∆τ , r -(sign(φ) n , r) = 0. (6.4)
The space discretization follows exactly the space discretization of the Level-Set transport equation. The function sign(φ) n , is approximated numerically by a smoothed function. This smoothness is important to obtain better properties of conservation and to insure stability (See [START_REF] Peng | A PDE-based fast local level set method[END_REF]):

sign(φ) n h = φ n h (φ n h ) 2 + h 2 ∇φ n h 2 . ( 6.5)
The discrete variational formulation can be written as:

Find Φ n+1 h ∈ Y h such that ∀r h ∈ Y h , Φ n+1 h -Φ n h oX n (w h n ,∆τ ) ∆τ , r h -(sign(φ) n h , r h ) = 0.
(6.6) 6.2. Mass conservation. The resolution of the transport equation of the Level-Set function causes the diffusion of a small amount of mass at each time step. It can be either an increase or a decrease of the error according to the topological changes of the interface. As time evolves these errors will typically accumulate. But the flow we considered is incompressible this implies that the volume occupied by any of the fluids should be preserved as well. There exists in the litterature many approaches that can be used in order to preserve the mass. We may cite for example Chang et al. [START_REF] Chang | A level set formulation of Eulerian interface capturing methods for incompressible fluid flows[END_REF], Sussman and Fatemi [START_REF] Sussman | An efficient interface-preserving level set redistancing algorithm and its application to interfacial incompressible fluid flow[END_REF].

In this paper, we follow the method proposed by Smolianski [START_REF] Smalonski | Numerical Modeling of Two-Fluid Interfacial Flows[END_REF], which seems simple, cheap and very efficient in our case. The simplicity of the method comes from the fact that the mass conservation can be enforced by adding a three lines algorithmic step. The key observation is that the error in mass balance should be very small within one time-step, usually this is done by using a sufficiently accurate scheme for the convection of level-set function. In our case we were able to reduce the volume error remarkably by using the new proposed boundary condition (as we will see in the numerical results), which makes the computational strategy remarkably cheap and efficient. The concept of the method is to vary the zero isocontour at each time step by moving the level-set function, i.e. by adding to Φ some signed constant c Φ , where |c Φ | is the distance between the old and new zero-level sets such that the new level-set function Φ new reduces the error of the corresponding mass and defines a new domain

Ω new 2 = {x ∈ Ω : Φ new > 0}.
The expression of c Φ is given by the formula

S exact -S(Ω 2 ) = Ω new 2 dx - Ω2 dx = Γ (c Φ n).ndΓ + O(c 2 Φ ) = c Φ Γ dΓ + O(c 2 Φ ), (6.7)
where S exact is the exact area (or volume in 3D) of the region occupied by the second fluid, S(Ω 2 ) is the numerical area of Ω In this experiment, a rectangular water column, initially at rest on a rigid horizontal plane is confined by the wall and a gate. The gate is suddenly removed at time t = 0 and the water starts to collapse due to gravity, the fluid spreads out and the height of the column falls. The initial height of water column H is 0.05715 m. The density of water is 1000 Kg.m -3 , its kinematic viscosity is 10 -6 m 2 .s -1 and its dynamic viscosity is 10 -3 Kg.m -1 .s -1 ; whereas the density of gas is 1 Kg.m -3 , its kinematic viscosity is 10 -5 m 2 .s -1 and its dynamic viscosity is 10 -5 Kg.m -1 .s -1 , the gravitational acceleration is g = 9.8m.s -2 .

2D experiment of Martin and Moyce

The computational domain is set to be 14H and 1.4H in length and height, respectively.

The space mesh sizes of two simulations are h = 1/50 and h = 1/100, respectively. The time mesh sizes of two simulations are ∆t = 2/50 and ∆t = 2/100, respectively.

In 2 dimensional problem, convergence was obtained for the regularity parameter ε = h/5. Throughout the simulation, the top surface of the computational domain is set to be outflow boundary condition. The simulated non-dimensional position of wave front l (l = l/H) and height of residual water column h (h = h/H) varying with non-dimensional time t = t g/H are compared with experimental results available in [START_REF] Martin | Some gravity wave problems in the motion of perfect liquids[END_REF]. 4), ( 5) and ( 6) illustrate a comparison between the experiment and the simulated interface between the gas phase and the liquid phase at different times, which proves that the present numerical method can be used to simulate large-density ratio immiscible two phase flow problems. Also, it can be noticed that the simulated interface shape and its motion are very similar to the experimental results. Figure [START_REF] Bell | A second-order projection method for variable-density flows[END_REF] shows the comparison of non-dimensional wave front between experiments and the two simulation results. From Figure [START_REF] Bell | A second-order projection method for variable-density flows[END_REF], it can be learned that the wave fronts of two simulations both agree well with the experiments. Figure [START_REF] Benqué | A finite element method for the Navier-Stokes equations[END_REF] gives the non-dimensional residual water column height varying with time. Simulation results from figures [START_REF] Bell | A second-order projection method for variable-density flows[END_REF] and [START_REF] Benqué | A finite element method for the Navier-Stokes equations[END_REF] indicate that the convergence study of the present numerical method is satisfied. It can be seen that the finer mesh h = 1/100 gives results closer to the experiment than the mesh h = 1/50 which is expected.

Figures (3), (

Dimensionless time (s)

Dimensionless time (s)

3D experiment of Martin and Moyce

The computational domain is set to be 8H,1H and 1H in length, width and height, respectively. The space mesh sizes of two simulations are h = 1/50 and h = 1/60, respectively. (For memory and cpu limitations we consider h = 1/60 for the fine mesh.)

The time mesh sizes of two simulations are ∆t = 2/50 and ∆t = 2/60, respectively. In 3 dimensional problems, convergence was obtained for the regularity parameter ε = h/5. The simulated position of wave front L and height of residual water column H varying with time t are compared with 3D experimental results of Martin and Moyce available in [START_REF] Saez | Etude numérique du remplissage 3D en fonderie[END_REF]. Figures [START_REF] Chang | A level set formulation of Eulerian interface capturing methods for incompressible fluid flows[END_REF], ( 10) and ( 11) illustrate the interface shape between the gas phase and the liquid phase at different times, which proves that the present numerical method can be used to simulate large-density ratio immiscible two phase flow problems for 3 dimensional problems. Furthermore, the simulated interface shape and its motion are very similar to the results presented in [START_REF] Saez | Etude numérique du remplissage 3D en fonderie[END_REF]. Figure [START_REF] Crandall | Viscosity solutions of Hamilton-Jacobi equations[END_REF] shows the comparison of wave front between experiments and the two simulation results.

From figure [START_REF] Crandall | Viscosity solutions of Hamilton-Jacobi equations[END_REF], it can be learned that the wave fronts of two simulations both agree well with the experiments.

Filling test case (2D).

In this case, we consider the initial domain

Ω ∪ Γ =]0, a[×]0, b[ ∪ ]c, c + d[×] -e, 0] ,
which is composed of a two-dimensional rectangular mould where a = 2 dm, b = 1 dm, d = 0.4 dm and e = 0.2 dm. We suppose that this mould contains a small amount of fluid defined by the Γ = [c, c + d] (as shown in the figure [START_REF] Dervieux | A finite element method for the simulation of a Rayleigh-Taylor instablity[END_REF] ).

The boundary of Ω is decomposed as

∂Ω = Γ ∪ Γ 0 ∪ Γ 1 ∪ Γ 2 where Γ 0 = [c, c + d] × {0} , Γ 1 = [0, c] × {0} ∪ {c, c + d} × [-e, 0] ∪ [c + d, a] × {0} ∪ {0, a} × [0, b] and Γ 2 = [0, a] × {b},
as shown in the figure [START_REF] Dervieux | A finite element method for the simulation of a Rayleigh-Taylor instablity[END_REF]. The density of water is 1 Kg.dm -3 , its kinematic viscosity is 10 -4 dm 2 .s -1 and its dynamic viscosity is 10 -4 Kg.dm -1 .s -1 ; whereas the density of gas is 10 -3 Kg.dm -3 , its kinematic viscosity is 10 -3 dm 2 .s -1 and its dynamic viscosity is 10 -6 Kg.dm -1 .s -1 .The gravitational acceleration is g = 0.98 dm.s -2 . Furthermore, we take U in = 0.4 dm.s -1 , h = 0.02, ∆t = 2h and ε = h/5. The non-dimensionalized incompressible Navier stokes equation can be written as: Figures [START_REF] Fourestey | Stabilité des méthodes de Lagrange-Galerkin du premier et du second ordre[END_REF], ( 17) and ( 18) show a comparison of the numerical results for the 4 considered cases at t = 0.168 s, t = 0.5 s and t = 1.5 s. Figures [START_REF] Fourestey | Stabilité des méthodes de Lagrange-Galerkin du premier et du second ordre[END_REF], ( 17) and [START_REF] Ganesan | A coupled arbitrary Lagrangian-Eulerian and Lagrangian method for computation of free surface flows with insoluble surfactants[END_REF] show that in case i, boundary layers phenomena appear between the free surface and the wall boundaries, which are corrected in case ii by using Navier-Boundary conditions. We still have a non-physical phenomenon concerning the right angle showed between the free surface and the wall boundaries. In the case iii, we impose natural boundary conditions under pressure effect that gives a natural angle between the free surface and the wall boundaries. We notice an increase of the volume of fluids among the cases, which will be justified later. In case iv, we present numerical results after the mass conservation correction step.

ρ ∂u ∂t + (u.∇)u - 1 Re div(2µ Du) + ∇p = -ρ 1 1 F r 2 e d , ( 7
In this section, we will show the effectiveness of the non-homogeneous Neumann boundary condition under pressure effect for the Level-Set function in reducing the volume dissipation at each time step. In order to do this, we analyze the volume error evolving with time:

err V = |V e -V n | |V n | and V e = V initial + |Γ 0 | U in t,
where V e is the exact volume, V n is the numerical volume and V initial is the initial volume. Figure [START_REF] Gerlach | Comparison of volume-of-fluid methods for surface tensiondominant two-phase flows[END_REF] shows a comparison between the volume error by using the classical method and the projec- To the left: classical method; To the right : projection method tion method for ε = h/5. The classical method and the projection method give very close results. The volume error in the classical method is better than that in the projection method. However, the CPU time and memory in the projection method are remarkably lower. Thus we will continue our work in 3D using the projection method. As for the comparison between the 4 considered cases, we may infer that the volume error reaches more than 17% in the case i, it decreases in the case ii with Navier boundary conditions wich allows the fluid to slip, it decreases remarkably in the case iii with the new proposed boundary condition under pressure effect. And because the error became small enough at each time step, we were able to apply the mentioned mass conservation algorithm in the case iv where the volume is conserved in the whole computational process as shown in figure [START_REF] Gerlach | Comparison of volume-of-fluid methods for surface tensiondominant two-phase flows[END_REF].

Figure [START_REF] Esmaeeli | Direct numerical simulations of bubbly flows. Part 1. Low Reynolds number arrays[END_REF] shows a comparison between the variation with time of the force acting on the interface and the one acting on the entry. The force acting on the interface is negligible and is very small with respect to the entry force. From figure [START_REF] Esmaeeli | Direct numerical simulations of bubbly flows. Part 1. Low Reynolds number arrays[END_REF], we show that the forces acting on the interface are negligible with respect to those acting on the entry, thus the pressure vanishes on the interface. Figure [START_REF] Hirt | An arbitrary Lagrangian-Eulerian computing method for all flow speeds[END_REF] shows the evolution of the error err V during time for the cases i, ii, iii and iv. It can be seen that the volume error follows the same way of variations as in 2D but it reaches a maximum of around 27% for the case i, it decreases with the case ii for a maximum about 14% and with the case iii for about a maximum of 9% and we are able to apply the simple mass conservation algorithm to preserve the mass property where the volume error is less than 0.1% during the whole computation in 3D. From figure [START_REF] Hirt | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF], we show that the forces acting on the interface are negligible with respect to those acting on the entry in 3 dimensional problem, thus the pressure vanishes on the interface. 30), two sequences are shown and are compared with the simulated interface. The figures present video excerpts at 0.25s intervals during the filling of two separate castings to ascertain the degree of reproducibility of the pour. The researchers were careful to ensure that the experimental conditions in each case were reproduced as exactly as possible. The shape and the size of the mould cavity as well as the experimental gating system are presented in figure [START_REF] Lamb | Hydrodynamics[END_REF].

In the experiment, the density of air is 1Kg.m -3 , the kinematic viscosity is set to be 1e -5 m 2 .s -1 and the dynamic viscosity is set to be 1e -5 Kg.m -1 .s -1 . The density of the aluminium liquid is 2385Kg.m -3 , the kinematic viscosity is set to be 0.55e -6 m 2 .s -1 and the dynamic viscosity is set to be 1e -3 Kg.m -1 .s -1 . The pouring temperature is 720 o C.

It should be noted that slight differences in behavior can be seen in experiments, even though they were designed to repeat identically. This is to be expected and is in the nature of turbulent phenomena.

It is worth noting that illustrations of modelling the experiment in three extreme cases: a) highly viscous laminar flow, b) low-viscosity laminar flow, and (c) turbulent flow are presented in [START_REF] Sirrell | The benchmark test[END_REF]. The nearest approximation to the experimentally observed results was consistent with high Reynolds numbers (40,000 or more). Most models predict the filling time correctly at around 2.0 seconds. In our simulation, our algorithm is not designed to simulate turbulent flow however an artificial viscosity for the air and metal liquid are used while keeping the same ratio of viscosity (18.18:1), since this ratio has a great influence on the shape of the interface. The density of air is 1Kg.m -3 , the kinematic viscosity is set to be 1.818e -3 m 2 .s -1 and the dynamic viscosity is set to be 1.818e -3 Kg.m -1 .s -1 ; whereas the density of the aluminium liquid is 2385Kg.m -3 , the kinematic viscosity is set to be 1e -4 m 2 .s -1 and the dynamic viscosity is set to be 2.385e -5 Kg.m -1 .s -1 . It follows that Reynold number is Re = 10000. The simulated filling time is 1.98 s which is very close to the filling time of the experiment. It is admitted that since the aluminum liquid is directly poured into the cavity during experiment, consequently, the inlet velocity should be actually varying with time. In current study, an average inlet velocity of 0.5 m.s -1 is used as in [START_REF] Saez | Etude numérique du remplissage 3D en fonderie[END_REF]. The averaged velocity is estimated from the gating system and numerical experiences. The mesh size is considered as h = 0.002, ∆t = 2h. For ease of mesh generation, the gating system is slightly modified, as shown in figure [START_REF] Lamb | Hydrodynamics[END_REF] to the right. We will see that the influence of this small modification of the gating system is negligible, and the simulated interface shape still agrees with the experimental results.

a lateral inlet is simulated. This experiment was originally developed by the corresponding author of these papers [START_REF] Chorin | On the convergence of discrete approximations to the Navier-Stokes equations[END_REF][START_REF] Crandall | Viscosity solutions of Hamilton-Jacobi equations[END_REF]. The experimental apparatus are mainly made of transparent Perspex. The schematic setup, the physical size of apparatus, and the boundary conditions are shown in Fig. 7. The length, width, and depth of the mould cavity are 300, 300, and 30 mm, respectively. The inlet velocity is designed to be 1.0 m/s. During the experiment, the liquid (water mixed with blue ink) is used for better visualization; and the filling process is recorded by high speed CCD video camera. During the simulation, the water density and the gas density are set to be 1000 and 1 kg m À3 , respectively; and the kinematic viscosity of water and gas are given to be 1eÀ6 and 1eÀ5 m 2 s, respectively. The grid number used in the current example is 100 Â 100 Â 10 in each axis direction, respectively. Fig. 8 shows the interface comparison between the experimental results and the simulation results. Figures (a), (c), (e), and (g) in Fig. 8 are high speed CCD camera captured interface at 0.35, 0.56, 0.67, and 1.35 s, respectively; whereas other sub-figures are the corresponding simulation results. The color bar in Fig. 8 represents the velocity component in Y-axis direction, namely, the inflow direction. And the color of vector plot represents different phases; specifically, the blue 1 color refers to the gas phase while the red color of vector arrow represents the water phase. Since this experiment is a real three dimensional flow problem, half of mould cavity in the width direction is blanked for better illustration. As shown in Fig. 8, at 0.56 s, there is a large gas bubble entrapped in the water phase, and at 0.67 s this large bubble severely deforms and another large bubble will be entrapped; these phenomena are clearly presented both in numerical simulation and experiment which demonstrates a good agreement; and at other times the phase interfaces between simulation and experiment also agree well. These facts demonstrate that the present sharp interface numerical method can be used to predicate the gas entrapment phenomenon during mould filling.

To demonstrate the capability of robustly handling topology changes of small gas bubbles, the motion of the entrapped bubble in Fig. 8(d) is investigated. As illustrated in Fig. 9, at 0.82 s the shape of the entrapped bubble is roughly similar to a sphere; at 0.87 s, due to the interaction of two phase flows, the bubble becomes more slender and the gas velocity inside the bubble also becomes larger; at 0.90 s the bubble is separated into two bubbles and the separated larger bubble moves upwards and finally starts to escape at 0.97 s. Aforementioned facts are reasonable and thus it can be concluded that the present particle level set method can conserve mass and be suitable for topological changing problems. Figures ( 27)-( 28)-( 29)-( 30)-( 31) and [START_REF] Olsson | A conservative level set method for two phase flow II[END_REF] show that the simulated interface location agrees with the experimental results at t = 0.74 s, t = 1 s, t = 1.24 s, t = 1.5 s, t = 1.74 s and t = 2 s. In this section, the domain Ω is a mould for an industry with a shape given in figure [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF]. In this simulation, The density of air is 1Kg.m -3 , the kinematic viscosity is set to be 1.818e -3 m 2 .s -1 and the dynamic viscosity is set to be 1.818e -3 Kg.m -1 .s -1 ; whereas the density of the aluminium liquid is 2385Kg.m -3 , the kinematic viscosity is set to be 1e -4 m 2 .s -1 and the dynamic viscosity is set to be 2.385e -5 Kg.m -1 .s -1 . The entry surface is 0.045 m 2 . The Reynolds number is Re = 10000. Furthermore, we take U in = 0.5 m/s, h = 0.02, ∆t = 2h and ε = h/5. The considered mesh contains 59364 vertices and 294544 tetrahedrons. Figures [START_REF] Peng | A PDE-based fast local level set method[END_REF], [START_REF] Pironneau | On the transport-diffusion algorithm and its applications to the Navier-Stokes equations[END_REF] and [START_REF] Pironneau | Stability and convergence of a Galerkin-characteristics finite element scheme of lumped mass type[END_REF] show the evolution of the interface at the initial time, t = 0.4 s, t = 0.8 s, t = 1.04 s, t = 2 s, t = 3 s, t = 4 s, t = 7.2 s and t = 8.4 s. Figure [START_REF] Qian | Molecular Hydrodynamics of the Moving Contact Line in Two-phase Immiscible Flows[END_REF] show the evolution of the error during the time for the case iv. One more time, the figure [START_REF] Qian | Molecular Hydrodynamics of the Moving Contact Line in Two-phase Immiscible Flows[END_REF] shows the efficiency of the algorithm even for practical cases.

Conclusion

In this paper, a three-dimensional algorithm for the modelling of interface transport scheme in two-phase flow in the cases of gas-water phases and gas-liquid metal phases was proposed, which simulates industrial cases. The capability of the algorithm to model two-phase problems with large-density and large-viscosity ratio is tested. The broken dam problem and Campbell Benchmark test are validated, it was shown that the simulated interface agrees well with the experimental results. In the future, we would like to perform many ameliorations for the algorithm. We are currently working on solving the 3D code on parallel computer using an iterative solver which reduces remarkably the time of the computation. Also we would like to improve the algorithm in order to be able to model turbulent flow which models the best gas-liquid metal interface shape. In addition, we will be working on introducing the temperature variation between the two phases and the mould in addition to the solidification process that follows the filling in modelling moulds casting for iron foundries.

Figure 1 .

 1 Figure 1. An arbitrary domain Ω

  [u]| Γ = 0 and [2µDu -pI]| Γ .n = 0, where [.]| Γ denotes the jump of quality across Γ in the normal direction of Ω 1 , i.e. [.]| Γ = .| Ω1 -.| Ω2 .



  Homogeneous Neumann boundary conditions ∂ n φ = 0 and Non-homogeneous Neumann boundary conditions under pressure effect ∂ n φ = G(u).

(4. 5 ) 4 . 2 .

 542 Discretization in space. Let τ h be a regular family of triangulations of Ω by triangles of tetrahedron k, of parameter h.

  and we denote by u n+1 h , p n+1 h et φ n+1 h respectively the discrete velocity, pressure and Level-Set function.

Figure 32 :

 32 Figure32: plan expérimental[Martin et Moyce,1952]

Figure 2 .

 2 Figure 2. Diagram of typical apparatus.

Figure 3 .

 3 Figure 3. On the top: experiment; On the bottom: approximated interface at t = 0s. IsoValue 0.025975 0.075925 0.125875 0.175825 0.225775 0.275725 0.325675 0.375625 0.425575 0.475525 0.525475 0.575425 0.625375 0.675325 0.725275 0.775225 0.825175 0.875125 0.925075 0.975025 Vec Value 0 0.082147 0.164294 0.246441 0.328588 0.410735 0.492882 0.575029 0.657176 0.739323 0.82147 0.903617 0.985764 1.06791 1.15006 1.2322 1.31435 1.3965 1.47865 1.56079

Figure 6 .

 6 Figure 6. On the top: experiment; On the bottom: approximated interface at t = 0.2s .

Figure 7 .

 7 Figure 7. Comparison of experimental and approximated wave front variation with collapsing time.

Figure 8 .

 8 Figure 8. Comparison of experimental and approximated height of residual water variation with collapsing time.

Figure 9 .Figure 10 .

 910 Figure 9. Approximated interface in 3 dimensions; To the left : t = 0.0317029s; To the right: t = 0.0634057s.

Figure 11 .

 11 Figure 11. Approximated interface in 3 dimensions; To the left : t = 0.253623s; To the right: t = 0.375905s.

Figure 13 .

 13 Figure 13. Comparison of experimental and approximated height of residual water variation with collapsing time in 3 dimensions.

Figure 14 . 1 Figure 15 .

 14115 Figure 14. The 2D mould .

. 1 )•

 1 where Re = u ref l ref /ν ref is the Reynolds number, F r = u ref / l ref g the Froude number, the nondimensional time t ref = t g/l ref with ρ ref the density of the first fluid ρ 1 , µ ref the viscosity of the first fluid µ 1 , l ref is the diameter of the inlet and u ref is the average velocity of the fluid. Then ρ ref = 1, µ ref = 10 -4 , l ref = 0.4, u ref = 0.4. It follows that the non-dimensional parameters are Re = 1600 and F r = 0.64. For the numerical tests, we consider four cases: Case i: Solution in the sense of viscosity with boundary conditions of Dirichlet type for the velocity and homogeneous Neumann type for the Level-Set function before mass conservation. • Case ii: Solution in the sense of viscosity with boundary conditions of Navier type for the velocity and homogeneous Neumann type for the Level-Set function before mass conservation. • Case iii: Solution in the sense of viscosity with boundary conditions of Navier type for the velocity and non-homogeneous Neumann type for the Level-Set function before mass conservation. • Case iv: Solution in the sense of viscosity with boundary conditions of Navier type for the velocity and non-homogeneous Neumann type for the Level-Set function after mass conservation.

Figure 17 .

 17 Figure 17. At t = 0.5 s, On the top: From left to right: case i, case ii. On the bottom: From left to right: case iii, case iv.

Figure 18 .

 18 Figure 18. At t = 1.3 s, On the top: From left to right: case i, case ii. On the bottom: From left to right: case iii, case iv.

Figure 19 .

 19 Figure[START_REF] Gerlach | Comparison of volume-of-fluid methods for surface tensiondominant two-phase flows[END_REF]. Comparison of the volume error for the cases i, ii, iii and iv for ε = h/5. To the left: classical method; To the right : projection method

Figure 20 .

 20 Figure 20. Variation of the force acting on the interface and the force acting on the entry with time.

Figure 21 .Figure 22 .Figure 23 .

 212223 Figure 21. To the left: The cube mould. To the right : The cube mesh.

Fig. 12 .Figure 26 .

 1226 Fig. 12. The schematic setup of the benchmark test.

Figure 27 .Figure 28 .Figure 29 .Figure 30 .Figure 31 .Figure 32 . 8 . 2 .

 27282930313282 Figure 27. Interface shape of the benchmark test at t = 0.74 s

Figure 33 .

 33 Figure 33. To the left: The mould shape; To the right: the mould mesh.

Figure 34 .

 34 Figure 34. The approximated interface at the initial time, t = 0.4 s and t = 0.8 s.

Figure 35 .Figure 36 .Figure 37 .

 353637 Figure 35. The approximated interface at t = 1.04 s, t = 2 s and t = 3 s.

  Introduction et description de l' expérience s résultats de cette expérience sont très souvent utilisés pour valider les modèles d' étude des oulements instationnaires très peu visqueux à surface libre. La géométrie du test permet la validation s modèles en 2D, mais étant réellement 3D, elle nous permettra de tester les deux cas de figure proche 2D et 3D. Nous ne citerons pas ici les nombreuses références aux résultats de cette périence, en choisissant de nous reporter directement aux données source[START_REF] Martin | Some gravity wave problems in the motion of perfect liquids[END_REF].considère une colonne d' eau au repos, retenue par une membrane de papier ciré (imperméable) trêmement fin, dans un réservoir parallélépipédique en matière plastique transparente (Figure32) tte membrane est maintenue par un film de cire à une bande métallique du réservoir. Elle est libérée moyen d' un courant électrique de haute intensité : lorsque le courant circule, la cire fond, et la lonne s' écroule dans un canal d' Altuglas. Le milieu ambiant est de l' air.

2 . By denoting L = Γ dΓ, we approximate c Φ by c Φ S exact -S(Ω 2 ) L(Γ) . (6.8) Then the corrected level set function becomes Φ new Φ + c Φ ||∇Φ||. (6.9) 6.2.1.1

  Figure 4. On the top: experiment; On the bottom: approximated interface at t = 0.07s.

	IsoValue 0.025975 0.075925 0.125875 0.175825 0.225775 0.275725 0.325675 0.375625 0.425575 0.475525 0.525475 0.575425 0.625375 0.675325 0.725275 0.775225 0.825175 0.875125 0.925075 0.975025
	Vec Value 0 0.0980932 0.196186 0.29428 0.392373 0.490466 0.588559 0.686652 0.784745 0.882839 0.980932 1.07902 1.17712 1.27521 1.3733 1.4714 1.56949 1.66758 1.76568 1.86377	IsoValue 0.025975 0.075925 0.125875 0.175825 0.225775 0.275725 0.325675 0.375625 0.425575 0.475525 0.525475 0.575425 0.625375 0.675325 0.725275 0.775225 0.825175 0.875125 0.925075 0.975025
		Vec Value 0 0.0980932 0.196186 0.29428 0.392373 0.490466 0.588559 0.686652 0.784745 0.882839 0.980932 1.07902 1.17712 1.27521 1.3733 1.4714 1.56949 1.66758 1.76568 1.86377
	082147 0.164294 0.246441 0.328588 0.410735 0.492882 0.575029 0.657176 0.739323 0.82147 0.903617 0.985764 1.06791 1.15006 1.2322 1.31435 1.3965 1.47865 1.56079	IsoValue 0.025975 0.075925 0.125875 0.175825 0.225775 0.275725 0.325675 0.375625 0.425575 0.475525 0.525475 0.575425 0.625375 0.675325 0.725275 0.775225 0.825175 0.875125 0.925075 0.975025
	IsoValue 0.025975 0.075925 0.125875 0.175825 0.225775 0.275725 0.325675 0.375625 0.425575 0.475525 0.525475 0.575425 0.625375 0.675325 0.725275 0.775225 0.825175 0.875125 0.925075 0.975025	Vec Value 0 0.082147 0.164294 0.246441 0.328588 0.410735 0.492882 0.575029 0.657176 0.739323 0.82147 0.903617 0.985764 1.06791 1.15006 1.2322 1.31435 1.3965 1.47865 1.56079
	Vec Value 0 0.0952957 0.190591 0.285887 0.381183 0.476478 0.571774 0.66707 0.762365 0.857661 0.952957 1.04825 1.14355 1.23884 1.33414 1.42943 1.52473 1.62003 1.71532 1.81062	IsoValue 0.025975 0.075925 0.125875 0.175825 0.225775 0.275725 0.325675 0.375625 0.425575 0.475525 0.525475 0.575425 0.625375 0.675325 0.725275 0.775225 0.825175 0.875125 0.925075 0.975025
		Vec Value 0 0.0952957 0.190591 0.285887 0.381183 0.476478 0.571774 0.66707 0.762365 0.857661 0.952957 1.04825 1.14355 1.23884 1.33414 1.42943 1.52473 1.62003 1.71532 1.81062

Figure 5. On the top: experiment; On the bottom: approximated interface at t = 0.1s .

  Figure 12. Comparison of experimental and approximated wave front variation with collapsing time in 3 dimensions.Figure[START_REF] Cuvelier | Some numerical methods for the computation of capillary free boundaries governed by the Navier-Stokes equations[END_REF] gives the dimensional residual water column height varying with time. Simulation results from figures (12) and (13) indicate that the convergence study of the present numerical method is satisfied in 3 dimensional problems.

		0.5
			experiment
			simulation h=1/50, dt=2/50
		0.45	simulation h=1/60, dt=2/60
		0.4
		0.35
	Wave front (L)	0.25 0.3
		0.2
		0.15
		0.1
		0.05
		0	0.05	0.1	0.15	0.2	0.25	0.3	0.35	0.4
						Time (s)
		0.06
									experiment
									simulation h=1/50, dt=2/50
									simulation h=1/60, dt=2/60
		0.05
		0.04
	Residual height (H)	0.03
		0.02
		0.01
		0
		0	0.05	0.1	0.15	0.2	0.25	0.3	0.35	0.4
						Time (s)

  Variation of the force acting on the interface and the force acting on the entry with time.

		0.25
										Forces acting on the interface
										Forces acting on the entry
		0.2
	Dimensionless Force	0.1 0.15
		0.05
		0	0	1	2	3	4	5	6	7	8	9	10
							Dimensionless Time
	Figure 25.

8. Numerical Results: gas-liquid metal modelling 8.1. Campbell benchmark. Since numerical tests in Section 1 are actually gas-water two-phase flow problems, to demonstrate the capability of the present numerical method for gas-liquid metal two phase flows during real casting mould filling process ( with a density ratio (2385:1) and viscosity ratio (18.18:1)). A benchmark test of filling process designed by Professor Campbell

[START_REF] Sirrell | The benchmark test[END_REF]

is simulated. The casting in the benchmark test is a plate sized 200mm × 100mm × 10mm in X, Y, and Z directions, respectively. The mould cavity is made of resin bonded sand. The pouring liquid is 99.999% pure aluminum. The filling process is recorded with in situ X-ray imaging technology. The instant when the stopper was lifted was designated as time zero. The experiment was repeated several times. In Figures

[START_REF] Lamb | Hydrodynamics[END_REF]

-(

Acknowledgements

Authors acknowledge support provided by the joint funding program from the Lebanese National Council for Scientific Research CNRS-L and the Université Saint-Joseph USJ (Lebanon) and the funding provided by Laboratoire Jacques-Louis Lions-Université Pierre et Marie-Curie (France).

The formula (6.8) is accurate up to O(c 2 Φ ). First, it is noteworthy that if S(Ω 2 ) > S exact , we have c Φ < 0 and the level-set function Φ will be moved downward. if S(Ω 2 ) < S exact , we have c Φ > 0 and the level-set function Φ will be moved upward. Remark 6.1. By using the reinitialization algorithm, which gives ∇Φ = 1, we can approximate Φ new by Φ new = Φ + c Φ . But numerically, we have numerical errors and we never reach the relation ∇Φ = 1 and it is better to use the previous expression (6.9).

Briefly, the proposed algorithm is described as the following: We start with the initial conditions. At each time step:

(1) We solve the Navier-Stokes equation.

(2) We advect the Level-Set function by the fluid velocity.

(3) We update the new density and viscosity. ( 4) We reinitialize the level-Set function.

(5) We apply the mass conservation algorithm.

Numerical results : gas-water modelling

In the following sections, we perform numerical simulations using the FreeFem ++ software [START_REF] Hecht | New development in FreeFem++[END_REF]. The Navier-stokes code used for the flow was verified in Freefem++ with Abboud [START_REF] Abboud | Schémas à deux-grilles pour la résolution du problème de Navier-Stokes instationnaire incompressible[END_REF] and was validated by the Turek benchmark in Freefem++ examples [START_REF] Schäfer | Benchmark Computations of Laminar Flow Around a Cylinder[END_REF]. The transport equation code used for the Level-set was validated and verified in Freefem++ with Pironneau [START_REF] Pironneau | Stability and convergence of a Galerkin-characteristics finite element scheme of lumped mass type[END_REF]. In all the numerical results, we consider the imaginary time step for the reinitialization ∆τ = ∆t/10. For Navier boundary conditions, we choose α = 0 and β = 1e -6, which allow slip without friction. 7.1. The broken dam problem. In order to validate the effectiveness of the proposed algorithm in modelling interface transport in two phase-flow problems especially it is capability to handle large density ratio (1000:1) and large-viscosity ratio (10:1), we consider the broken dam experiment done by J.C. Martin and W. J. Moyce [START_REF] Martin | Some gravity wave problems in the motion of perfect liquids[END_REF]. The schematic setup of the experiment is illustrated in figure [START_REF] Babuska | The Finite Element Method with Penalty[END_REF].

apitre II

Résolution du problème de