A Cartesian Method with Second-Order Pressure Resolution for Incompressible Flows with Large Density Ratios

: An Eulerian method to numerically solve incompressible biﬂuid problems with high density ratio is presented. This method can be considered as an improvement of the Ghost Fluid method, with the speciﬁcity of a sharp second-order numerical scheme for the spatial resolution of the discontinuous elliptic problem for the pressure. The Navier–Stokes equations are integrated in time with a fractional step method based on the Chorin scheme and discretized in space on a Cartesian mesh. The biﬂuid interface is implicitly represented using a level-set function. The advantage of this method is its simplicity to implement in a standard monoﬂuid Navier–Stokes solver while being more accurate and conservative than other simple classical biﬂuid methods. The numerical tests highlight the improvements obtained with this sharp method compared to the reference standard ﬁrst-order methods.


Introduction
Bifluid problems are ubiquitous in nature and in many industrial applications such as combustion in engines, water waves energy converters, and jet printers to cite only a few. In such applications, the density ratio between the two fluids can be large, for instance the ratio is equal to 1000 between water and air. Accurate numerical modeling and numerical simulations of these kind of phenomena are then necessary, in particular to optimize such devices.
In this paper we are thus concerned with the numerical modeling of incompressible bifluid flows with large density ratios, like air and water, and by the accurate description of the phenomena occurring at their interface. We present a sharp Cartesian method for the simulation of incompressible flows with high density and viscosity ratios. This method is an extension of the second-order Cartesian method for elliptic problems with immersed interfaces developed in [1].
Cartesian grids are an attractive alternative to body fitted meshes. Indeed, they avoid complex mesh generation as well as mesh adaptation when unsteady interfaces are considered. Moreover, the numerical resolution of the governing equations can be simplified with an easy parallelization and the use of standard linear algebra libraries. Generally speaking, numerical schemes are easy to implement on a Cartesian mesh because a dimensional splitting is often possible; however, some numerical modeling is necessary near a complex interface that does not fit the background Cartesian grid. This is the case for fluid structure interface and moreover for bifluid interface where the properties of the flow are discontinuous. Indeed, applying naively a numerical scheme originally devised for a flow with constant or continuously varying density will lead to a non-consistent treatment of the interface. Most of the time, it will result in severe stability issues if the density ratio is large as highlighted in [2] and references therein; therefore, as already mentioned, one has to devise specific numerical schemes at the vicinity of the interface. This region is called narrow band and is the set of numerical points that have at least one neighbor on the other side of the interface.
Conservative or non-conservative approaches can both be used to face this issue. Among the non-conservative approaches, one solution is to regularize the properties of the fluids in the vicinity of the interface, so that the density, viscosity, and their derivatives are continuous in the whole computational domain. This idea leads to the well-known "Continuous Surface Force" (CSF) method [3], where the discontinuous quantities are smoothed near the interface, and in case of a fluid with surface tension, this surface tension is taken into account as a smooth volume force instead of a surface force. This method is widely used (see for instance [4,5]) because it offers a straightforward way to implement the presence of two fluids in an already existing monofluid Navier-Stokes code. However, the exact way that the regularization should be performed is not always clear, and spurious oscillations at the bifluid interface can appear due to errors in the pressure gradient computations. Another non-conservative method introduced by Kang, Fedkiw, and Liu [6] after the CSF is the Ghost Fluid Method (GFM). It is based on a firstorder method developed in [7] to solve an immersed interface elliptic problem, with a dimensional splitting making the method easy to implement. The resulting linear system is symmetric and has the same structure as the usual matrix to discretize a Poisson equation with variable coefficient on a Cartesian grid. This method has been used successfully in a lot of works, for instance [8,9]. One drawback is that the method is only first-order accurate near the interface [2] and a loss of conservativity of the momentum of each fluid near the interface can occur leading to erroneous velocities.
Non-conservative methods are often associated with a level-set representation of the interface [4] because the level-set method is itself intrinsically non-conservative at the discrete level, and convenient to use on a Cartesian grid.
The other family of methods is based on the conservative form of the Navier-Stokes equations, where mass and momentum fluxes of each fluid are explicitly computed, see for instance [2,[10][11][12][13]]. An explicit interface representation is necessary even if the interface do not coincide with grid points. Conservative methods are generally more stable than non-conservative methods. The price for this increased stability is an additional amount of numerical developments due to the interface reconstruction, which can be performed from information carried by Lagrangian markers or by cell quantities such as volume fractions.
Another approach has been developed recently [14,15] to deal with large density ratio. In this approach, a fully second-order method is obtained at the interface (for both velocity and pressure) with several physical boundary treatments, including velocity and traction boundary conditions.
In this paper we aim to preserve as much as possible the simplicity of the Ghost Fluid method of [6], avoiding an explicit identification of the volume fractions near the interface, while improving the accuracy and stability of the pressure computation. We thus propose a method, mainly based on the improvement of the Ghost Fluid method, where the discontinuities across the interface are taken into account in a sharp way with a secondorder scheme inspired from [1]. This second-order treatment improves the conservativity of the method, as it will be proved numerically in the section devoted to numerical validations.
After having described the governing equations for the incompressible bifluid flows that we consider (Section 2), the discretization of these equations in each fluid and at the interface are presented (Sections 3 and 4). The second-order numerical resolution of the elliptic problem arising from the computation of the pressure is introduced (Section 5), and the overall is validated on several two-dimensional numerical test cases (Section 6).

Flow Equations
We consider a rectangular domain Ω filled with two viscous incompressible fluids with different densities and viscosities. The subdomains Ω − and Ω + corresponding to the two fluids are separated by an interface Γ as depicted in Figure 1.
In this work, these domains are implicitly defined with a scalar function φ, usually called the level-set function, see Section 2.2, that takes different values in each subdomain with a fixed value on the interface. For instance we chose φ = 0 on Γ, φ > 0 in Ω + and φ < 0 in Ω − . The unit normal to the interface is denoted n and the unit tangent vector is denoted η. The density is denoted and the viscosity is denoted where H is the Heavyside function, i.e., The flow is modeled in the whole domain with the incompressible Navier-Stokes equations: with g the gravitational acceleration vector, τ the viscous stress tensor: and the term σκ∇H accounting for the surface tension effects, with σ the surface tension itself and κ the local curvature of the interface between the fluids. This formulation of the bifluid incompressible Navier-Stokes equations contains a singular term which is not trivial to handle. Alternatively, the flow can also be modeled in each subdomain with the incompressible Navier-Stokes equations: The above equations are completed by jump conditions at the interface Γ between the two fluids. In what follows, jumps are defined by the notation [ψ] = ψ + − ψ − . The first jump conditions describe the balance between the normal stresses at the interface and the surface tension σ, with κ the local curvature of the interface Γ, [µ(∇u · n, ∇v · n) · η + (∇u · η, ∇v · η) · n] = 0.
Others jump conditions can be derived from continuity properties across the interface. For instance, for a viscous fluid, the velocity field is continuous across the interface [v] = 0.
Since the material derivative of (6) and (7) is zero, one can write which leads to ∇p The jump condition for the pressure p can be simplified. We differentiate the jump on the velocity in the tangential direction: Moreover, because the velocity is divergent-free on each side of the interface, 0 = [∇ · u] = [(∇u · n, ∇v · n) · n + (∇u · η, ∇v · η) · η].
Finally, we will use Equations (9) and (12) to compute the pressure jump at the interface.

Interface Description
The interface between the two fluid subdomains is implicitly defined by a scalar function φ. Local geometrical information on the function φ is needed to obtain an accurate discretization of the interface. To this purpose we use the level set method, introduced by Osher and Sethian [16] and described in [17][18][19]. A common choice for the level-set function φ is the signed distance function to the interface: The zero isoline of φ thus represents implicitly the interface Γ immersed in the computational domain.
We assume that the interface is smooth enough, so that the derivatives of the level-set function in the vicinity of the interface are well-defined. A useful property of the level-set function is a straightforward computation of its normal with the formula In the same way, the curvature of the interface can be computed with the formula We impose a curvature threshold 1/h, with h the grid spacing, corresponding to the minimum size of a bubble for a given spatial discretization.
For a moving interface as it is usually the case for bifluid problems, the flow density and viscosity are updated with φ tracking the interface thanks to the transport equation where the velocity fieldsû coincides with the flow velocity field u on the interface Γ. Different choices for the value ofû in Ω + and Ω − can be a priori used. A natural choice we have considered isû = u in the whole domain: Another possible choice is the extension velocity introduced in [20].

Navier-Stokes Monofluid Solver and Numerical Method for Interface Evolution
The computational domain is discretized on a two-dimensional uniform Cartesian grid with a grid spacing ∆x = ∆y = h. However, the following approach stands for nonuniform Cartesian meshes. The points on the Cartesian grid are defined as M i,j = (x i , y j ). Similarly, we denote by u ij the approximation of u at the point (x i , y j ). In what follows, all the unknowns are collocated in space on Cartesian meshes. An odd-even coupling can sometimes be observed when one uses collocated unknowns on a Cartesian grid. This problem can be fixed using some corrections such as [21,22]. However, in all the applications considered in this paper we have not observed any odd-even coupling.
In this section we present the numerical solver devised for monofluid incompressible Navier-Stokes equations that is used in each fluid when the interface is not taken into account. We also provide information about the numerical methods used to compute the evolution of the level-set function.

Flow Computation
We use a classical projection method [23,24] to solve the Navier-Stokes equations in each fluid outside the narrow band. In what follows, a non-incremental projection method is used, i.e., the guess value for the pressure in the prediction step is zero.
We thus compute successively: The convective terms are computed with a fifth-order WENO scheme, and the viscous terms with an explicit second-order centered finite-difference scheme. The time integration is performed with a first-order explicit Euler scheme, which is consistant with the use of a non-incremental version of the projection scheme.
The pressure appearing in Equation (19) is computed through the resolution of a Poisson equation in order to enforce the divergence-free condition. At each point in one subdomain, the following relationship is satisfied: We will provide additional details about the jump conditions that have to be satisfied across the interface for this problem in Section 4. On the exterior boundary of the domain Ω, Neumann boundary conditions are satisfied: where u b is the value of the velocity to be imposed on the external boundaries. We compute at each iteration an adaptive time step taking into account the restrictions due to convection, viscosity, and surface tension. The convective time step restriction is given by with |u| max and |v| max the maximum magnitudes of the horizontal and vertical velocities. The viscous time step restriction is given by We also apply a time step restriction associated with the surface tension evaluated only in the narrow band. This time step restriction is similar to the one in [6] and in [8], and is in this context usually the most restrictive: Finally, at each time step, the overall algorithm is the following: • Prediction: evaluate convective and diffusive fluxes and compute u * , • Interface evolution: convect the level-set with velocity u and re-initialize if necessary, • Construction and resolution of the linear system for the pressure, • Correction step: update velocity with pressure gradient.

Numerical Method for the Level-Set Evolution
The computation of the level-set function should be performed very accurately when one deals with moving interfaces. Indeed, as the level-set method is not intrinsically conservative, a lack of accuracy in the computation of the level set evolution results often in a substantial loss of mass for one of the fluids. It can also increase the problem of transfer of momentum between both fluids and generate spurious velocity oscillations. Moreover, if one wants to compute the curvature of the interface from (15), the level-set function needs to be accurate enough (at least third-order) so that the finite difference formulas used to discretize (15) are consistent.
Unfortunately, the property of the signed distance function is usually lost when the interface evolves with the flow velocity. The norm of the level set gradient can be far from unity. These gradients variations of the level-set are harmful to the accuracy of the numerical evaluation of the normal to the interface and the curvature. To circumvent the problem, Sussman et al. [4] introduced a reinitialization algorithm to recover the signed distance function through the resolution of the eikonal equation Several methods have been developed over the years to perform this reinitialization step, either by using a relaxation method and searching a stationary solution to a time dependent Hamilton-Jacobi equation [4] or using a Gauss-Seidel-based method as in fast marching methods [25,26] or fast sweeping methods [27].
The usual level set strategy for an evolving interface is thus the following: • Use a transport equation to update φ. • From time to time, reinitialize φ with the signed distance function.
One of the most widespread option in the literature is to combine a fifth-order WENO scheme [28] with a RK3 scheme, for the transport and for the reinitialization through a relaxation equation. It provides a high-order yet stable resolution; however, although the reinitialization procedure, performed with such a numerical scheme, may improve mass conservation, it also introduces some error by slightly moving the interface as shown in [29]. For this reason, Russo and Smereka [29] introduced a subcell fix taking into account the interface location in the reinitialization procedure. This technique was extended to a higher order accuracy in [30] through the use of third-order ENO schemes near the interface.
Usually, this reinitialization steps are performed uniformly every each n iterations. A recent study [31] proposes a strategy to sample the reinitialization steps based on interface deformation criteria. This should also be coupled with the high-order decentered reinitialization scheme of [30] near the interface.
In what follows, we use the classical option of the fifth-order WENO scheme for the spatial discretization. Since it is the most commonly used technique in the literature, it will allow us to distinguish the effects of the new scheme for the pressure computation from the effects of the reinitialization technique.

Notations
Let us introduce some definitions and notations. A grid point is defined to be irregular if at least one of its neighbors is on the other side of the interface, i.e., if the sign of φ changes between this point and at least one of its neighbors, see Figure 2. The set of irregular points is called the narrow band. All the other points are called regular grid points.
We define the interface point I i,j,E = (x i,j,E , y j ) as the intersection of the interface Γ and the East segment [M ij M i+1j ], if it exists. Similarly, the interface points I i,j,W = (x i,j,W , y j ), respectively. With this notation the same interface point can be described in two different ways The set of interface points is denoted Γ h , see Figure 2 for an illustration.

Modeling Choices for the Discontinuities across the Interface
In this study, the values of the viscosity and the densities are discontinuous across the interface; therefore, if the numerical scheme for solving incompressible Navier-Stokes equations described previously in Section 3.1 was applied on the irregular grid points of the narrow band, the approximations of the following terms would not be consistent:

•
Prediction step: viscous terms. • Correction step: divergence of the predicted velocity, elliptic operator, and gradient of the pressure.
This lack of consistency could eventually lead to stability problems. The computation of the convective terms is not mentioned in the above enumeration because it is performed with a fifth-order WENO scheme, which automatically provides spatial adaptivity. Therefore, we assume that the gradients computed with the WENO scheme are decentered near the interface, and consequently, consistent. Moreover, the levelset function is classically evolved with such a scheme, because it is crucial to have a good accuracy in the computation of the interface evolution; thus, it seems coherent to have the same numerical scheme for the convection of the interface and the convection of the fluids.
In this work, we use two different strategies to handle the lack of consistency near the interface, one for the viscous terms and another for the pressure computation. For the viscous term in the prediction step, we follow a continuous approach and regularize the quantities used for the computation of the viscous terms. It has been proven in [32,33] that this continuous approach provides correct accuracy for high Reynolds number flows. It has also been used successfully in [2,3]. A sharp approach for the viscous terms could probably improve the accuracy of the simulations; however, the complexity of the computations would be increased due to the treatment of the jump conditions for the viscous terms (5) implying derivatives of the velocity components in both normal and tangential directions. Moreover, if one needs to use an implicit treatment of the viscous terms, such a sharp treatment would become more complex to handle.
The discretization of the prediction step becomes: In practice the viscosity and the inverse of the density were regularized by a discrete convolution [2]: Then, we discretize the viscous terms with a classical second-order centered scheme as in Section 3.1.
In the correction step, according to the jump conditions (4)- (12) presented in Section 2, the pressure satisfies an elliptic problem with discontinuous values of the solution and its derivative across the interface: Because the viscous terms in the prediction step are handled with a regularization approach, we have [µ] = 0 and ∇·τ ρ = 0. Therefore, the pressure computed for the correction step satisfies rather: [p] = σ κ on Γ, The details of the resolution of this elliptic problem will be provided in Section 5. Let us remark however that as the values of p are discontinuous across the interface, two values of p will be created at each point of the interface considered in the numerical scheme.

Gradient and Divergence for Correction Step
The predicted velocity u * obtained after the prediction step (18) is defined only on grid points. We need to compute the divergence of this predicted velocity in order to solve the elliptic Equation (20). However, since the two fluids have different properties across the interface and the derivatives of the velocity are not necessarily continuous, we need to use a decentered stencil on each side of the interface. Consequently, we have to compute two values for u * on each interface point, one for each side of the interface. In practice, as jump conditions for u * are not available, we perform simply linear extrapolations from the grid values on the interface points. Then, to compute the divergence of u * on an irregular grid point M i,j , we use a standard five point stencil, see Figure 3. Formally this is equivalent to a standard first-order decentered scheme.
More precisely, we denote u * S the value of the predicted velocity u * on the nearest point to M i,j in the south direction (possibly an interface point), with coordinates (x S , y S ). Similarly, we define u * N , u * W and u * E and the associated coordinates (x N , y N ), (x W , y W ) and (x E , y E ). The discretization reads Similarly, in order to keep a consistent discretization, the gradient of the pressure p appearing in Equation (19) is also computed with an adapted decentered stencil near the interface, see Figure 3. More precisely, with the same notations S, N, W and E, as before, the discretization reads If one of the discretization point is an interface point, we consider the value of the pressure on this point corresponding to the same subdomain than point M i,j . Indeed we recall that, since the pressure is discontinuous across the interface, two pressure unknowns (one for each subdomain) are computed at an interface point. If no interface point is involved, the numerical scheme reduces to the classical second-order central finite differences scheme. Figure 3. Example of geometrical configuration, with points involved in the discretization of the divergence of the predicted velocity and pressure gradient at grid points M i,j and M i+2,j+2 in black.

Numerical Resolution of Elliptic Problems with Immersed Interfaces
The elliptic problem with discontinuous values across an interface (28)-(30) is solved with the second-order method developed in [1]. The interface points are used to impose the jump conditions across the interface in the numerical scheme. Since the pressure is discontinuous across the interface, two unknowns are created, one for each side of the interface. The accuracy of this method is based on the use of unknowns located at the interface. The size of this linear system is thus augmented with two unknowns for each interface point. These interface unknowns are used to discretize the flux jump conditions and the elliptic operator accurately enough to get a second-order convergence in maximum norm. Actually, to this purpose, near the interface the elliptic operator needs to be discretized with a first-order truncation error, and the fluxes with a second-order truncation error. For a visual explanation of the discretization we refer to Figure 4. The advantage of using this method, compared to the reference work of [6] is that the jump conditions in the correction step are solved with second-order accuracy instead of firstorder. The drawback is that the linear system is not symmetric anymore and it is solved with the preconditioned GMRES method.

Discrete Elliptic Operator
We use a standard five-point stencil including the grid point M i,j and its nearest neighbors in each direction: interface or grid points. More precisely, we denote p S the value of the solution on the nearest point in the south direction, with coordinates (x S , y S ).
Similarly, we define p N , p W , and p E and the associated coordinates (x N , y N ), (x W , y W ), and (x E , y E ). Since the density is piecewise constant, the discretization reads where ρ ± stands for ρ + or ρ − .

Discrete Flux Transmission Conditions
As written before, at each interface point we create two additional unknowns, called interface unknowns. We denote them by p ± i,j,γ with γ = E, W, N or S. The interface unknowns carry the values of pressure on each side of the interface.
Contrarily to [1], we do not have a jump condition on the normal derivative, but on the whole gradient, as expressed in formula (30). To obtain the same number of equations and unknowns we have to chose in which direction we want to project this gradient equality. As we use a Cartesian grid, it is easier to discretize the x− and y−derivatives than derivatives in other directions. Consequently, we discretize the following jump conditions at each interface point I i,j,γ , with γ = N, S, W, E.
We want the truncation error of the discretization of flux equality (34) and (35) to be second-order accurate in order to solve the problem with a second-order accuracy. A possible configuration of the interface is illustrated in Figure 4. In the x-direction, it is straightforward to compute a second-order approximation of the x-derivative with three a priori nonequidistant points. For example we approximate the flux on the left side of interface point I i,j,E , if it exists, with the values of p on the points M i−1,j , M i,j and I i,j,E with the formula: The right x-derivative is approximated in the same way.
The same discretization holds for the y-derivative. The formulas (36) and (37) are consistent if both grid points involved in the formula, for instance M i−1,j and M i,j , belong to the same domain. If on one side of the interface the two closest grid points aligned with the intersection point do not belong to the same subdomain, then the second-order discretization is not possible anymore. In this case, we use a first-order discretization involving only two points: the interface point and the closest grid point on the same side of the interface. Such a case is illustrated on Figure 4. In fact, this first-order discretization is equivalent to the Ghost Fluid method [6].
Let us notice that, because we use a dimensional splitting for the jump conditions across the interface, it is quite straightforward to eliminate the interface unknowns from the linear system. We simply inject expressions (36) and (37) in the jump condition (34), and use the resulting equality to express p ± i,j,E as a function of p i−1,j , p i,j , p i+1,j and p i+2,j . This expression for p ± i,j,E can then be used in the discretization of the elliptic operator (32). The local curvature κ i,j,γ at the interface point I i,j,γ is computed in the following way. We first compute on all irregular grid points the value with centered second-order finite-difference formulas. Then we perform a one-dimensional linear interpolation of these values on the interface points.

Numerical Results and Validations
This section is devoted to the numerical simulations and validations. After having studied the well-balanced behavior of the proposed approach, i.e., the capacity of the method to preserve equilibria state for a bubble, in Section 6.1, we perform a quantitative validation for the well-known dam break problem from Section 6.2. We then present other numerical simulations for the rising of bubbles with different sizes that can qualitatively been compared to some reference results in Section 6.3.

Equilibria Preservation for a Bubble: the Parasitic Oscillations
This first test case aims to assess the influence of the interface curvature error on the stability of the numerical scheme. A bubble is located at the center of the computational domain in Sections 6.1.1 and 6.1.2. Due to the Laplace law and the concavity of the interface, the pressure inside the bubble is larger than the pressure outside. If the curvature of the interface is computed numerically, the errors due to the numerical approximation in the right-hand side of Equation (33) cause small errors in the resolution of the pressure equation and the system is thus not well balanced. These errors create artificial values of the velocity near the interface which should theoretically be zero. These artificial velocities are often called parasitic currents. The amplitude of the parasitic currents is an indication of the stability and the accuracy of the numerical method, and especially of the pressure computational step. Indeed, they are the only source of numerical errors.
In what follows several comparisons with other references methods, the Ghost Fluid and the CSF methods in Section 6.1.1 and the Volume of Fluid method in Section 6.1.2, are presented on slightly different test cases.

Comparison with the Ghost Fluid and the CSF Methods
We use the same parameters as in [8], where a Ghost Fluid and a CSF method were implemented. The amplitude of the parasitic currents, compared to the results in [8], are reported in Table 1. The L ∞ and L 2 norms are computed over the whole domain Ω. The initial configuration is described in Figure 5. As it can be observed, the amplitude of the parasitic currents generated by our method is several orders of magnitude smaller than those of the CSF method, and significantly lower than those of the Ghost Fluid method when the grid is refined.  We now compare the behavior of our method to the Volume of Fluid method developed in [12]. The density and viscosity ratio are both chosen to be one for this test case.
The coefficient σ is chosen so as to obtain an Ohnesorge number Oh = µ σρD satisfying Oh 2 = 1 12000 . The maximum velocity is computed for varying grids at non-dimensional A comparison between our method and the Volume of Fluid method is presented in Table 2. The new method provides a better accuracy than the Volume of Fluid method for the coarsest grid. As expected, the Volume of Fluid method outperforms our new approach for finer meshes due to more sophisticated schemes near the interface. Nonetheless, Table 2 show a second-order accuracy for our new method.

Collapse of a Water Column: the Dam Break Problem
This test case is studied in [2,34], and based on experiments conducted in [35]. The initial configuration is a water column at rest in air. The initial height and width of the column are both 5.715 cm. The domain size is 40 cm×10 cm. The physical constants are the same than for the rising bubble (Section 6.1.1). For more details, we refer the reader to [2]. We present in Figure 6 the interface evolution at non-dimensional times T = t g/H = 0, 1, 2, 3, 4, with H the initial height of the water column. The computations are performed with 256 × 64 points. Figure 7 present the temporal evolution of the water front, compared to the experimental results [35], to results with the Ghost Fluid method used for pressure resolution, and to the conservative method of Raessi and Pitsch [2]. We observe that the front propagation is in agreement with the experimental results and the results of the conservative method [2]. This means that, though the method is not strictly conservative, the numerical errors due to momentum transfer across the interface are not large enough to slow down the propagation of the front. It is not the case for instance for the Ghost Fluid method, as it can be noticed in Figure 7 and has been reported in [2].

Rising of Air Bubble in Water
We study the evolution of fluid bubbles rising in an heavier fluid, and compare our results to several methods in the literature. The initial configuration is described in Figure 8.

Comparison with the Ghost Fluid Method
We consider air bubbles rising in water, as in the test case proposed for the Ghost Fluid method in [6]. The value of the physical parameters are ρ water = 1000 kg/m 3 , µ water = 1.137 × 10 −3 kg/ms, ρ air = 1.226 kg/m 3 , µ air = 1.78 × 10 −5 kg/ms, σ = 0.0728 kg/s 2 , g = −9.8 m/s 2 . (41) We consider two cases: a small bubble with R = 1/300 m and a large one R = 1/3 m. In the first case, the surface tension plays an important role in the evolution of the interface because of the high bubble curvature. In the second case, the surface tension has less influence, and larger deformations occur. The interface of the small bubble and the vorticity values are plotted at times t= 0, 0.02, 0.035, 0.05 in Figure 9. The interface of the large bubble and the vorticity values are plotted at times t= 0, 0.2, 0.35, 0.5 in Figure 10. Our numerical results are in good agreement with [6]. 6 R 9 R in [6] or 10 R in [36] 2 R in [36] 3 R in [6] gaz liquid R Figure 8. Initial fluid domain for the test case of the rising bubble in water in the References [6,36].

Comparison with SPH and the Level-Set Method
This test case is taken from [36], and inspired from a test case presented in [4]. It gives us the opportunity to compare our method to another class of methods, based on the SPH formulation. The initial configuration is described on Figure 8. The values of the physical parameters are .025 m, ρ water = 1000 kg/m 3 , µ water = 1.137 × 10 −3 kg/ms, ρ air = 1.226 kg/m 3 , µ air = 1.78 × 10 −5 kg/ms, σ = 0.0728 kg/s 2 , g = −9.8m/s 2 .
The evolution of the interface and the vorticity values are plotted on Figure 11, for 120 × 200 grid points. We observe that the interface deforms in a way similar to the results in [36].

Two Air Bubbles in Water
Two air bubble are initially at rest in water, see Figure 12. At the same time that they are rising, their interaction produces larger deformations than for a single bubble. This test case is meant to assessing the conservativity of the new method, as the increase in pression resolution is meant to increase this conservativity.
The value of the physical parameters are In Figure 13, the interface evolution and the vorticity values for the Ghost Fluid method (used for the pressure computation) and the new method are plotted. One observes small oscillations of the vorticity near the interface between two fluids for the Ghost Fluid method, but not for the new method. The conservation of the partial mass of air and momentum of air in each dimension for the new method, and the Ghost Fluid method are depicted in Figures 14-16 respectively. As air is the less dense fluid, it is more prone to big oscillations if erroneous mass and momentum transfers between the two fluids occur. One can observe that with the new method and its second-order pressure resolution, the oscillations of these quantities are notably less present than for the Ghost Fluid. The partial mass is better conserved too.

Conclusions
We have developed a new method on Cartesian grids for the simulation of incompressible flows with large density ratios. This method relies on a sharp resolution of the pressure term across the interface defined by a level-set function. The advantage of the proposed approach is its simplicity to implement in an existing Cartesian monofluid Navier-Stokes solver as for instance the Ghost Fluid or CSF methods. It is only necessary to modify the stencil for the pressure equation at some irregular grid points by adding one additional point. This Cartesian scheme uses thus additional unknowns located on the interface to discretize with second-order accuracy the jump conditions across the interface. The viscous term is treated with a regularizing approach that allows us to eliminate terms in the jump conditions. Indeed, it has been shown [32,33] that the regularization has no significant impact on the accuracy of the results. Numerical results show that this method leads to more accurate and stable results than some reference methods as for instance the well known Ghost Fluid or CSF methods. The conservation of volume of each phase is also better conserved. We thus take advantage of both the regularity of the interface defined by the level-set function (curvature is numerically consistent), and the mass conservation of each phase. The CLSVOF [11] can also provide good mass conservation and precision of the interface but with a more complex numerical implementation.
Future works include an extension of the method to three-dimensional problems with interactions with solids including the numerical simulation of wave energy converters [37,38]. We also aim to study the sensibility of the reinitialization procedure of the level-set function defining the bifluid interface, as for instance the one developed in [31].