Benjamin Lalanne

Lucia Rueda Villegas

Sébastien Tanguy
email: tanguy@imft.fr

Frédéric Risso

On the computation of viscous terms for incompressible two-phase flows with Level Set/Ghost Fluid Method

Keywords:

In this paper, we present a detailed analysis of the computation of the viscous terms for the simulation of incompressible two-phase flows in the framework of Level Set/Ghost Fluid Method when viscosity is discontinuous across the interface. Two pioneering papers on the topic, Kang et al. [10] and Sussman et al. [26], proposed two different approaches to deal with viscous terms. However, a definitive assessment of their respective efficiency is currently not available. In this paper, we demonstrate from theoretical arguments and confirm from numerical simulations that these two approaches are equivalent from a continuous point of view and we compare their accuracies in relevant test-cases. We also propose a new intermediate method which uses the properties of the two previous methods. This new method enables a simple implementation for an implicit temporal discretization of the viscous terms. In addition, the efficiency of the Delta Function method [24] is also assessed and compared to the three previous ones, which allow us to propose a general overview of the accuracy of all available methods. The selected test-cases involve configurations wherein viscosity plays a major role and for which either theoretical results or experimental data are available as reference solutions: simulations of spherical rising bubbles, shape-oscillating bubbles and deformed rising bubbles at low Reynolds numbers.

Introduction

In pioneering works on simulation of incompressible two-phase flows with interface tracking or interface capturing algorithms, discontinuities of pressure field, density field or viscosity field were tackled by smoothing the singular terms across the interface on a finite layer of about 2 or 3 grid cells thick. This kind of approaches has been proposed in the framework of Front-Tracking algorithms [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF][START_REF] Unverdi | A front-tracking method for viscous, incompressible multi-fluid flows[END_REF], for Volume-of-Fluid methods [START_REF] Gueyffier | Volume-of-fluid interface tracking with smoothed surface stress methods for three-dimensional flows[END_REF][START_REF] Lafaurie | Modelling merging and fragmentation in multiphase flows with SURFER[END_REF], and for Level-Set methods [START_REF] Sussman | A Level Set approach for computing solutions to incompressible two-phase flow[END_REF]. In the literature, they are referred to as Continuum Surface Force methods or Delta Function methods. Such methods still provide satisfactory results in many situations. Nevertheless, they suffer from some drawbacks such as the presence of intense spurious currents at high capillary-numbers, which can induce strong numerical instabilities. Also, they have to face the tricky problem of describing a surface as a volume of small thickness. A few years ago, a new class of method appeared in the literature, where singular terms across an interface are treated as jump conditions, which prevents from interface thickening as well as it reduces spurious currents. Significant contributions to this new approach, well known as Sharp Interface method, can be found in [START_REF] Popinet | A front-tracking algorithm for accurate representation of surface tension[END_REF] for Front-Tracking methods, in [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF][START_REF] Liu | A boundary condition capturing method for Poisson's equation on irregular domains[END_REF][START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF] for Level-Set methods, and [START_REF] François | A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework[END_REF] Volume-of-Fluid methods. In [START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF], a Level-Set method was coupled with a Volume-of-Fluid method, but the way to impose the jump conditions depended only on the Level-Set function -the volume fraction being just used to improve interface advection and to compute interface curvature. Most of these papers emphasize on improvements resulting from a sharp description of the surface tension, but they are not conclusive about the effectiveness of Sharp Interface methods for the discretization of viscous terms. However, Kang et al. [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF] proposed a new discretization of the viscous terms based on the general principle of the Ghost Fluid Method [START_REF] Fedkiw | A non-oscillatory Eulerian approach to interfaces in multimaterial flows (The Ghost Fluid Method)[END_REF][START_REF] Liu | A boundary condition capturing method for Poisson's equation on irregular domains[END_REF]. The main feature of this new approach consists in computing the contraction of the divergence of the viscous stress tensor in a Laplacian operator of the velocity components µ V rather than its full expression (with µ the dynamic viscosity and V the velocity field). This contraction is only possible for incompressible flow when viscosity is uniform; it is hence apparently not suited to compute the viscous terms on grid cells that are crossed by an interface between two different fluids. To solve this issue for two-phase flows, once µ V has been computed, the authors propose to add explicitly the jump of the normal viscous stresses at the interface, resulting from the jump of viscosity, to the pressure jump while solving the Poisson equation for the pressure. Another approach to compute the viscous terms in the framework of Sharp Interface methods, inspired of a previous work [START_REF] Liu | A boundary condition capturing method for Poisson's equation on irregular domains[END_REF], was proposed by Sussman et al. [START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF]. In a similar way of what it is usually done in Continuum Surface Force methods, the full expression of the divergence of the viscous stress tensor ∇.(2µD) is discretized, the viscous stress tensor being defined by 2µD for a Newtonian fluid, with D the rate-of-deformation tensor. In this approach, contrary to previous work [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF], it seems that the contribution of the normal viscous stresses has not to be added explicitly to the pressure jump condition at the interface. Valuable numerical validations of this method are proposed in [START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF]. However, whether these two approaches are formally equivalent or not remained an open question. In particular, the fact that jump of the normal viscous stresses has to be explicitly added to the pressure jump in [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF] and not in [START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF] still required to be justified.

In this paper, we focus on the computation of viscous terms for incompressible two-phase flows using Level-Set methods.

In the framework of Sharp-Interface methods, we shall demonstrate that the two approaches introduced in [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF] and [START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF] are equivalent from a theoretical point of view, even if their numerical implementation is quite different.

From this original demonstration, a new intermediate method is proposed, which allows to make more easily the implicit temporal discretization of the viscous terms as the two methods proposed in [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF] and [START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF].

To complete this survey of numerical methods, we also discuss the consistency of these three Sharp-Interface methods with the Delta-Function method that is based on a Continuum Surface Force approach.

Section 2 describes the formalism of this set of four methods and proves their equivalence. Then, Section 3 proposes three numerical benchmarks that are particularly relevant to assess the reliability of these methods and compare their accuracy: spherical rising bubbles, damped shape oscillations of a bubble, and deformed rising bubbles at low Reynolds number.

Equations & numerical methods

Equation of motion for two-phase flows, jump condition and singular source term

Let us consider the Navier-Stokes equations for a single-phase incompressible flow; the following two equations express respectively momentum and mass conservation:

ρ D V Dt = -∇ p + ∇ • (2µD) + ρ g, (1) ∇ • V = 0, ( 2
)
where V is the velocity field, D Dt is the Lagrangian derivative, ρ is the fluid density, µ is the fluid viscosity, p is the pressure field, g is the gravity acceleration, and D is the rate-of-deformation tensor defined as:

D = ∇ V + ∇ V T 2 . ( 3
)
For incompressible flows with constant viscosity, it can be easily shown that:

∇ • (2µD) = µ V , (4)
which leads to the following simplified expression for Navier-Stokes equation:

ρ D V Dt = -∇ p + µ V + ρ g. ( 5
)
Let us now consider an incompressible two-phase flow, the whole-domain-formulation [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF][START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF] leads to the following expression that includes the capillary forces:

ρ D V Dt = -∇ p + ∇ • (2µD) + σ κ nδ Γ + ρ g, ( 6
)
where σ is the surface tension, κ is the interface curvature, n is the normal vector at the interface and δ Γ is a multidimen- sional Dirac distribution localized at the interface. If no phase change occurs, eq. ( 2) stays valid. As density and viscosity can vary the across interface, they can be defined as:

ρ = ρ gas + (ρ liq -ρ gas )H Γ = ρ gas + [ρ] Γ H Γ , (7)
µ = µ gas + (µ liq -µ gas )H Γ = µ gas + [µ] Γ H Γ , (8)
where H Γ is a Heaviside distribution, equal to 1 in the liquid phase and equal to 0 in the gas phase, with the following definitions for density and viscosity jump conditions:

[ρ] Γ = ρ liq -ρ gas ,

[µ] Γ = µ liqµ gas .

(

) 10
Moreover, considering a divergence-free condition on the velocity field, the following unsteady advection equation for the density field is obtained from mass conservation:

Dρ Dt = 0. ( 11
)
From that equation, we can deduce the usual propagation equation for interface motion that can be used with Level Set methods:

Dφ Dt = 0, ( 12
)
where φ is a signed distance function from the interface. Note that another similar equation can be derived for Volume Of Fluid methods:

DF Dt = 0, ( 13
)
where F is the local volume fraction of the considered phase in a computational cell. The balance of normal stresses at the interface involves a jump of pressure and a jump of normal viscous stresses:

[p] Γ = σ κ + 2[µ] Γ ∂ V n ∂n . (14)
Let us remark now that the viscous-stress tensor can be split into two parts:

∇ • (2µD) = µ∇ • (2D) + 2D • ∇µ. ( 15
)
If viscosity is a piecewise constant, by using the definition of the Heaviside derivative, it can be shown that:

∇µ = [µ] Γ nδ Γ . (16)
Finally, we observe that the previous splitting ( 15) is a decomposition of the global viscous term in a continuous part and a singular part (zero everywhere except at the interface):

∇ • (2µD) = µ V + 2[µ] Γ D • nδ Γ (17)
This highlights that the contribution of the viscous terms on the pressure jump is already included in the left-hand side of eq. ( 17). Thus, we can suppose that computing the divergence of the overall viscous stress tensor allows including the contribution of viscous effects in the pressure jump. A more rigorous proof of that assertion is provided in Section 2.4 of this paper.

To summarize, the momentum equation can be defined in a whole-domain formulation leading to eq. ( 6) with a density and a viscosity field defined respectively by eqs. [START_REF] Hnat | Spherical cap bubbles and skirt formation[END_REF] and [START_REF] Hong | On boundary condition capturing for multiphase interfaces[END_REF], or in a jump-condition formulation leading to eq. ( 5) with additional jump conditions (9), [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF], [START_REF] Leveque | Immersed interface methods for Stokes flow with elastic boundaries or surface tension[END_REF].

About the link between jump conditions and singular source terms

Let us consider a volume Ω, bounded by a surface Σ , and defined as Ω 1 U Ω 2 , where the two subdomains Ω 1 and Ω 2 are separated by a closed surface Γ with Γ ∩ Σ = ∅, as illustrated in Fig. 1. We define n 1 and n 2 as the normal vectors at the surface Γ pointing respectively towards Ω 1 and Ω 2 .

Let us now consider a piecewise continuous function ψ which respects the Poisson equation,

ψ = 0 ( 18
)
inside both subdomains Ω 1 and Ω 2 , with the following jump condition at the interface:

[ψ] Γ = a( x Γ ), (19)
where a( x Γ ) is a function of the interface coordinate x Γ . The jump condition (19) involves that the following surface integral is equal to zero

Γ [ψ] Γ -a( x Γ ) n 2 dS = 0 ( 20
)
Considering a homogeneous Dirichlet boundary condition on the external boundary Σ of the computational field, it can be easily demonstrated that

Σ ψ n Σ dS = 0, (21)
where n Σ is the normal vector pointing outwards the external boundary Σ . Thus

Σ ψ n Σ dS + Γ [ψ] Γ -a( x Γ ) n 2 dS = 0. (22)
Moreover, by remarking that

Γ [ψ] Γ -a( x Γ ) n 2 dS = Γ ψ n 2 dS + Γ ψ n 1 dS -Γ a( x Γ ) n 2 dS ( 23
)
we obtain the following equality:

Σ ψ n Σ dS + Γ ψ n 2 dS + Γ ψ n 1 dS = Γ a( x Γ ) n 2 dS. ( 24
)
Eq. ( 24) can also be written as

Ω1 ∇ψdΩ + Ω2 ∇ψdΩ = Ω ∇ψdΩ = Ω a( x Γ ) n 2 δ Γ dΩ, ( 25
)
which leads to the following local formulation:

∇ψ = a( x Γ ) n 2 δ Γ . (26)
Finally, after differentiation, we obtain:

ψ = ∇ • a( x Γ ) n 2 δ Γ . (27)
We have therefore shown that adding the following source term ∇ • (a( x Γ ) n 2 δ Γ ) to the right-hand side of a Poisson equation imposes a jump condition [ψ] Γ = a( x Γ ) on its solution ψ . This is also valid for a Poisson equation like eq. ( 18) with a right-hand side different from zero. Moreover, following a similar approach, it can also be demonstrated that a scalar field that is solution of the following equation ψ = b( x Γ )δ Γ , [START_REF] Trontin | Direct numerical simulation of a freely decaying turbulent interfacial flow[END_REF] with b( x Γ ) another function depending on interface coordinate, will also fulfills the following jump condition:

[ n 2 • ∇ψ] = b( x Γ ). ( 29
)
More accurate demonstrations on the derivation of jump conditions can be found in [START_REF] Leveque | Immersed interface methods for Stokes flow with elastic boundaries or surface tension[END_REF].

Projection method for whole-domain formulation and jump-condition formulation

Let us consider now a first-order explicit projection method for the single-phase incompressible Navier-Stokes equation. It consists in prescribing a first intermediate velocity field which does not respect the divergence-free condition (predictor step):

V * 1 = V n -t V n • ∇ V n -υ V n -g , (30)
with υ the kinematic viscosity.

Next, a Poisson equation must be solved to determine the pressure field:

p n+1 = ρ n+1 ∇ • V * t . (31)
Then, in a last step, the final divergence-free velocity field can be computed using the pressure field determined at the previous step:

V n+1 = V * -t ∇ p n+1 ρ n+1 . ( 32
)
Applying the projection method to an incompressible two-phase flow is less straightforward. Depending on whether the whole-domain formulation or the jump-condition formulation is used, the intermediate velocity in the first step of the projection method has to be different.

Let us first consider a whole-domain formulation. An explicit temporal discretization of the capillary and viscous terms leads quite naturally to the following predictor step

V * 2 = V n -t V n • ∇ V n - ∇ • (2µD n ) ρ n+1 - σ κ nδ Γ ρ n+1 -g , (33)
with the corresponding Poisson equation

∇ • ∇ p n+1 ρ n+1 = ∇ • V * 2 t = ∇ • V n -t V n • ∇ V n - ∇ • (2µD n ) ρ n+1 - σ κ nδ Γ ρ n+1 -g . (34)
Next, the correction step gives:

V n+1 = V * 2 -t ∇ p n+1 ρ n+1 . ( 35
)
Let us now consider the application of the projection method to a jump-condition formulation of the Navier-Stokes equation.

We have to solve eq. ( 5), so the first step of the projection method is similar to eq. ( 30) for a single-phase flow. As we want to impose the appropriate jump condition at the interface on the pressure field (eq. ( 14)), we can show, using the results presented in Section 2.2, that it can be done by solving the following Poisson equation for the pressure with a right-hand side split into a continuous part and a singular part:

∇ • ∇ p n+1 ρ n+1 = ∇ • V * 1 t + ∇ • (σ κ + 2[µ] Γ ∂ V n ∂n ) nδ Γ ρ n+1 . ( 36
)
Using the splitting of the viscous term given in eq. ( 17) and noting that n

• D • n = ∂ V n
∂n , we find the following relation between V * 2 and V * 1 :

V * 2 = V * 1 + t (σ κ + 2[µ] Γ ∂ V n ∂n ) nδ Γ ρ n+1 . ( 37
)
The use of eq. ( 37) allows demonstrating that eq. ( 34) and eq. ( 36) are formally identical:

∇ • ∇ p n+1 ρ n+1 = ∇ • V * 2 t = ∇ • V * 1 t + ∇ • (σ κ + 2[µ] Γ ∂ V n ∂n ) nδ Γ ρ n+1 . ( 38
)
Then, by considering eq. (37), we can deduce the following correction step for a jump-condition formulation, leading to the same field V n+1 than in eq. ( 35):

V n+1 = V * 1 - t ρ n+1 ∇ p n+1 -σ κ + 2[µ] Γ ∂ V n ∂n nδ Γ . (39)
To summarize, we have shown that both a whole-domain formulation and a jump-condition formulation can be used when applying a projection method to an incompressible two-phase flow. These two methods are formally equivalent, but they lead to a different definition of the predicted velocity field in the first step of the projection method. In [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF], the intermediate velocity field V * 1 is used, whereas in [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF][START_REF] Chang | A level set formulation of eulerian interface capturing methods for incompressible fluid flows[END_REF][START_REF] Gueyffier | Volume-of-fluid interface tracking with smoothed surface stress methods for three-dimensional flows[END_REF][START_REF] Sussman | A Level Set approach for computing solutions to incompressible two-phase flow[END_REF] it is the intermediate velocity field V * 2 . These two methods will be respectively named in the rest of this paper the Ghost-Fluid Primitive viscous Method (GFPM) and the Delta-Function Method (DFM).

Projection method for intermediate formulations

Other definitions of the predicted velocity field are possible. For example, we propose the following two definitions, the respective interests of which will be discussed later:

V * 3 = V n -t V n • ∇ V n - ∇ • (2µD n ) ρ n+1 -g , (40) V * 4 = V n -t V n • ∇ V n - ∇ • (µ∇ V n ) ρ n+1 -g . ( 41
)
In the same way as in the previous section, they lead to the definition of these two Poisson equations, formally equivalent to eq. ( 34) and eq. ( 36):

∇ • ∇ p n+1 ρ n+1 = ∇ • V * 3 t + ∇ • σ κ nδ Γ ρ n+1 , (42)
∇ • ∇ p n+1 ρ n+1 = ∇ • V * 4 t + ∇ • (σ κ + [µ] Γ ∂ V n ∂n ) nδ Γ ρ n+1 . ( 43
)
The following relations have been used to determine eq. ( 43):

∇ • (2µD) = ∇ • µ ∇ V + ∇ T V = ∇ • (µ∇ V ) + µ∇ • ∇ T V + ∇µ • ∇ T V , (44)
∇ • (2µD) = ∇ • (µ∇ V ) + [µ] Γ ∇ T V . nδ Γ , (45)
since ∇.∇ T V = 0 for an incompressible flow, and

n.∇ T V . n = ∂ V n ∂n . (46)
Then, the two final corrector steps, respectively suitable with eq. (40) and eq. ( 41), are:

V n+1 = V * 3 - t ρ n+1 ∇ p n+1 -σ κ nδ Γ , (47)
V n+1 = V * 4 - t ρ n+1 ∇ p n+1 -σ κ + [µ] Γ ∂ V n ∂n nδ Γ . (48)
The intermediate velocity field V * 3 has already been used in Refs. [START_REF] François | A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework[END_REF][START_REF] Popinet | A front-tracking algorithm for accurate representation of surface tension[END_REF][START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF] in the framework of Sharp Interface methods. In the rest of this paper, the Ghost-Fluid Conservative viscous Method (GFCM) will refer to the numerical method proposed in [START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF]. Finally, the use of V * 4 leads to an original numerical method to deal with the viscous terms, that we name the Ghost-Fluid Semi-Conservative viscous Method (GFSCM).

The four methods DFM, GFPM, GFCM and GFSCM described above are therefore theoretically equivalent, but they lead to different numerical approximations. The next objective of this paper is to discuss about their respective interests and numerical accuracy.

Spatial discretization of surface tension

When a whole-domain formulation [START_REF] Chang | A level set formulation of eulerian interface capturing methods for incompressible fluid flows[END_REF][START_REF] Gueyffier | Volume-of-fluid interface tracking with smoothed surface stress methods for three-dimensional flows[END_REF][START_REF] Lafaurie | Modelling merging and fragmentation in multiphase flows with SURFER[END_REF][START_REF] Sussman | A Level Set approach for computing solutions to incompressible two-phase flow[END_REF][START_REF] Unverdi | A front-tracking method for viscous, incompressible multi-fluid flows[END_REF] is used, the surface-tension singular term is accounted for when V * 2 is computed. The intermediate velocity field is next differentiated to compute the right-hand side of the pressure Poisson equation. Therefore, it is well suited to a smooth discretization of the surface-tension term that involves a fictitious interface thickness. Its application to a sharp-interface description would be trickier. This approach is that of the Delta-Function Method, known as Continuum Surface Force method. In this framework, the numerical discretization of discontinuous terms makes use of the following regularized smoothed Heaviside function defined from the Level-Set function φ:

H ε (φ) =      0 if φ < -ε 1 2 1 + φ ε + sin( πφ ε ) π if |φ| ≤ ε 1 if φ > +ε, ( 49
)
where ε is the fictitious interface thickness, which is equal to 2 or 3 cell sizes. The surface-tension term can then be computed using the following approximation:

σ κ nδ Γ ∼ = σ κ∇ H ε . (50)
At the opposite, a projection method using a jump-condition formulation with V * 1 as intermediate velocity field is well suited for Sharp-Interface methods. Thanks to this formulation, artificial smoothing of interface singularities can be avoided, allowing a sharp-interface description. It has been used in the framework of Ghost-Fluid Method in [START_REF] Liu | A boundary condition capturing method for Poisson's equation on irregular domains[END_REF], where authors propose a first-order discretization to impose a jump condition to the solution as well as to the normal derivative of the solution when solving the Poisson equation

∇ • (β∇u) = f , (51)
i.e. with the following interface jump conditions

[u] = a( x Γ ) and β ∂u ∂n = b( x Γ ).
Considering that capillary effects can be computed as a jump condition on pressure, an original Sharp-Interface method for solving the Navier-Stokes equation in the framework of Level-Set methods was introduced in [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF]. This method has been widely used to perform direct numerical simulations of two-phase flows [START_REF] Menard | Coupling Level Set/VOF/Ghost Fluid methods: validation and application to 3D simulation of the primary break-up of a liquid jet[END_REF][START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF][START_REF] Tanguy | Application of a Level Set method for simulation of droplet collisions[END_REF][START_REF] Trontin | Direct numerical simulation of a freely decaying turbulent interfacial flow[END_REF][START_REF] Van Poppel | A Ghost Fluid Level Set methodology for simulating multiphase electrohydrodynamic flows with application to liquid fuel injection[END_REF]. Let us now briefly summarize this approach.

The discretization of a 2D Poisson equation as eq. ( 51), which involves a discontinuity of the solution and of its normal derivative at the interface, writes:

β i+1/2 j ( u i+1 j -u i j x ) -β i-1/2 j ( u i j -u i-1 j x ) x + β i j+1/2 ( u i j+1 -u i j y ) -β i j-1 2 ( u i j -u i j-1 y ) y = f i j + g i j + h i j , (52)
where interpolations of the diffusion coefficient β on the cell borders are computed as the harmonic average of its value β + in the region where the Level Set function is positive and its value β -in the region where the Level Set function is negative; for example, if the interface crosses a mesh segment between x i and x i+1 :

β i+1/2 j = β + β - β + θ + β -(1 -θ) , (53) with θ
= |φ i+1 j | |φ i+1 j | + |φ i j | . (54)
The source terms g i j and h i j in the right-hand side of eq. ( 52) match with the discontinuities of u and of β ∂u ∂n ; they can be computed as

g i j = g R i j + g L i j + g T i j + g B i j , ( 55
)
h i j = h R i j + h L i j + h T i j + h B i j , (56)
where each of these terms are non-zero when an interface crosses the corresponding grid segment, and superscripts R, L, T , B denote respectively right, left, top and bottom sides. Their expressions are:

g R i j = ± β i+1/2 j a Γ x 2 g L i j = ± β i-1/2 j a Γ x 2 g T i j = ± β i j+1/2 a Γ y 2 g B i j = ± β i j+1/2 j a Γ y 2 , ( 57
) h R i j = ± β i+ 1 2 j b Γ θ β ± x h L i j = ∓ β i-1 2 j b Γ θ β ± x h T i j = ± β i j+ 1 2 b Γ θ β ± y h B i j = ∓ β i j-1 2 b Γ θ β ± y . (58)
It is noteworthy that an analogy can be established by using the results demonstrated in Section 2.2 and the discretization described in this section. For sake of simplicity, let us take a diffusion coefficient β equal to unity. In Section 2.2, it has been shown that solving

u = f with [u] = a( x Γ ) and ∂u ∂n = b( x Γ ) is equivalent to solve u = f + ∇. a( x Γ ) nδ Γ + b( x Γ )δ Γ .
Comparing this expression with discrete equation (52), we deduce that source terms g i j and h i j , that appear in the framework of the Ghost-Fluid Method, are respectively numerical approximations of ∇.(a

( x Γ ) nδ Γ )| i j and b( x Γ )δ Γ | i j : ∇ • (a nδ)| i j = g i j + O ( x) ( 59
) bδ| i j = h i j + O ( x) (60)
Such an analogy establishes the link between the Continuum Surface Force approach based on numerical approximations of δ, and the Ghost-Fluid Method that accounts for source terms in the discretization of Poisson equations.

Spatial discretization of viscous terms

More details are provided in this section about the spatial discretization of viscous terms for the different methods described in the previous sections. In the following, for sake of simplicity, all developments are presented in 2D with Cartesian coordinates and a uniform grid. The generalization to axisymmetric coordinates or 3D geometry is straightforward.

Conservative viscous formulation

If a conservative viscous formulation ∇.(2µD) is considered, as it is the case in the GFCM or in the DFM, the projection of the divergence of the viscous-stress tensor leads to the following relations:

∇ • (2µD)| x = ∂ ∂x 2µ ∂u ∂x + ∂ ∂ y µ ∂u ∂ y + ∂ v ∂x (61) ∇ • (2µD)| y = ∂ ∂x µ ∂u ∂ y + ∂ v ∂x + ∂ ∂ y 2µ ∂ v ∂ y ( 62
)
Let us remark that these two relations depend on both velocity components. This is why a fully implicit temporal discretization of viscous terms is more complex when this formulation is used. Suitable approximations of eq. ( 61) lead to the following spatial discretizations, D x and D y denoting the differential operators:

∇ • (2µD) • e x | i+1/2 j ∼ = D x (2µD x u)| i+1/2 j + D y (µD y u)| i+1/2 j + D y (µD x v)| i+1/2 j (63) D x (2µD x u)| i+1/2 j = 2µ i+1 j ( u i+3/2 j -u i+1/2 j x ) -2µ i j ( u i+1/2 j -u i-1/2 j x ) x (64) D y (µD y u)| i+1/2 j = µ i+1/2 j+1/2 ( u i+1/2 j+1 -u i+1/2 j y ) -µ i+1/2 j-1/2 ( u i+1/2 j -u i+1/2 j-1 y ) y (65) D y (µD x v)| i+1/2 j = µ i+1/2 j+1/2 ( v i+1 j+1/2 -v i j+1/2 x ) -µ i+1/2 j-1/2 ( v i+1 j-1/2 -v i j-1/2 x ) y ( 66
)
The following spatial discretization for eq. ( 62) is

∇ • (2µD) • e y | i j+1/2 ∼ = D x (µD y u)| i j+1/2 + D x (µD x v)| i j+1/2 + D y (2µD y v)| i j+1/2
(67)

D x (µD y u)| i j+1/2 = µ i+1/2 j+1/2 ( u i+1/2 j+1 -u i+1/2 j y ) -µ i-1/2 j+1/2 ( u i-1/2 j+1 -u i-1/2 j y ) x (68) D x (µD x v)| i j+1/2 = µ i+1/2 j+1/2 ( v i+1 j+1/2 -v i j+1/2 x ) -µ i-1/2 j+1/2 ( v i j+1/2 -v i-1 j+1/2 x ) x (69) D y (2µD y v)| i j+1/2 = 2µ i j+1 ( v i j+3/2 -v i j+1/2 y ) -2µ i j ( v i j+1/2 -v i j-1/2 y ) y ( 70
)
This spatial discretization for the conservative viscous approach can be used both in the framework of Smoothed-Interface methods (e.g. in the DFM) [START_REF] Chang | A level set formulation of eulerian interface capturing methods for incompressible fluid flows[END_REF][START_REF] Sussman | A Level Set approach for computing solutions to incompressible two-phase flow[END_REF] or Sharp-Interface methods (e.g. in the GFCM) [START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF]. The difference lies in the way viscosity is approximated on the border of a cell that is crossed by an interface. In the DFM, eq. ( 8) is used with the approximation of the Heaviside function given by eq. ( 49).

In the GFCM proposed by Sussman et al. [START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF], the interpolation used to compute viscosity on the cell border is quite similar to that proposed by Liu et al. [START_REF] Liu | A boundary condition capturing method for Poisson's equation on irregular domains[END_REF] to interpolate the diffusion coefficient of a Poisson equation (given here by eq. ( 53)). Thereafter, this kind of interpolation will be used for all the computations that involve the Ghost-Fluid Method (i.e. GFCM, GFPM and GFSCM).

Ghost-Fluid Primitive viscous Method

The GFPM of Kang et al. [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF] is now briefly described. When this approach is used, the jump condition on the normal viscous stresses is added to the capillary forces when solving the Poisson equation for pressure as shown in eq. ( 36).

This approach for the viscous term can be summarized by considering the following projections for the Laplacian of the velocity:

µ V | x = µ ∂ 2 u ∂x 2 + ∂ 2 u ∂ y 2 , (71) µ V | y = µ ∂ 2 v ∂x 2 + ∂ 2 v ∂ y 2 . ( 72
)
The corresponding discretizations can be expressed by using operators D x and D y :

(µ V ) • e x | i+1/2 j ∼ = µ D 2 x u| i+1/2 j + D 2 y u| i+1/2 j , (73)
(µ V ) • e y | i j+ 1 2 ∼ = µ D 2 x v| i j+ 1 2 + D 2 y v| i j+ 1 2 . (74)
In [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF], remarking that the velocity gradients are strongly discontinuous at the interface because of viscosity jumps, it was proposed to compute each component of the Laplacian of the velocity µ V (eqs. ( 73) and ( 74)) by computing the different terms of the vector ∇.(µ∇ V ) and subtracting their jump at the interface; on the x-axis,

µD 2 x u| i+1/2 j = D x (µD x u)| i+1/2 j -[µD x u]δ Γ | i+1/2 j , ( 75
)
µD 2 y u| i+1/2 j = D y (µD y u)| i+1/2 j -[µD y u]δ Γ | i+1/2 j , (76)
and on the y-axis, µD 2

x v| i j+1/2 = D x (µD x v)| i j+1/2 -[µD x v]δ Γ | i j+1/2 , ( 77
)
µD 2 y v| i j+1/2 = D y (µD y v)| i j+1/2 -[µD y v]δ Γ | i j+1/2 . ( 78
)
The appropriate jump conditions for eqs. ( 75), ( 76), ( 77) and ( 78) are given in the following matrix by:

[µ ∂u ∂x ] [µ ∂u ∂ y ] [µ ∂ v ∂x ] [µ ∂ v ∂ y ] = [µ]∇ V 0 t T 0 t + [µ] n T n∇ V n T n -[µ] 0 t T 0 t ∇ V n T n, ( 79
)
where t is the unit vector tangent to the interface.

To compute accurately eqs. ( 75), ( 76), ( 77), (78) in the framework of the Ghost-Fluid Method [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF], suitable sharp approximations, eq. ( 58), of the Dirac distribution are required.

Note that the jump conditions in eq. ( 79) are only valid for a continuous velocity field between the two phases: this approach cannot be used when phase changes or mass transfers occur.

Ghost Fluid Semi-Conservative viscous Method

In the GFSCM, updating V * 4 requires computing the following viscous terms:

∇ • (µ∇ V )| x = ∂ ∂x µ ∂u ∂x + ∂ ∂ y µ ∂u ∂ y , ( 80
) ∇ • (µ∇ V )| y = ∂ ∂x µ ∂ v ∂x + ∂ ∂ y µ ∂ v ∂ y . ( 81
)
That can be done in a similar way as for the conservative viscous approach, except for the second-order cross-derivatives, which are non-zero only at the interface and are therefore included in the pressure jump when the Poisson equation is solved. As pointed out previously, this approach can be used with an implicit temporal discretization of the viscous terms and also when a phase change occurs.

Implicit temporal discretization of viscous terms

Further elements are provided now to deal with an implicit temporal discretization of viscous terms.

When the incompressible Navier-Stokes equations for a single-phase flow are considered, the projection method leads to the following semi-implicit algorithm for the predictor step

V * 1 -tυ V * 1 = V n -t V n • ∇ V n , (82)
which requires solving an independent linear system for each velocity component. When a two-phase flow with a discontinuous viscosity is considered, the use of a semi-implicit projection method to solve the incompressible Navier-Stokes equations is more complex. Indeed, in 2D, the following two coupled equations must be solved simultaneously, leading to a more complicated linear system as it is carried out in [START_REF] Sussman | A coupled level and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows[END_REF][START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF]:

ρ n+1 i+1/2 j u * i+1/2 j -t D x 2µD x u * + D y µD y u * + D y µD x v * i+1/2 j = ρ n+1 i+1/2 j u n i+1/2 j -t V n • ∇ u n | i+1/2 j , (83) ρ n+1 i j+1/2 v * i j+1/2 -t D x µD y u * + D x µD x v * + D y 2µD y v * i j+1/2 = ρ n+1 i j+1/2 v n i j+1/2 -t V n • ∇ v n | i j+1/2 .
(84) Some previous works report the following temporal discretization to split the resolution of this system:

ρ n+1 i+1/2 j u * i+1/2 j -t D x 2µD x u * + D y µD y u * i+1/2 j = ρ n+1 i+1/2 j u n i+1/2 j -t V n • ∇ u n i+1/2 j + t D y µD x v n i+1/2 j , (85) ρ n+1 i j+1/2 v * i j+1/2 -t D x µD x v * + D y 2µD y v * i j+1/2 = ρ n+1 i j+1/2 v n i j+1/2 -t V n • ∇ v n i j+1/2 + t D x µD y u n i j+1/2 . (86)
Nevertheless, it is doubtful that the numerical scheme composed of eqs. ( 85), (86) really alleviates the time step constraint imposed by viscosity since the last terms on the right is treated explicitly.

Another way to carry out an implicit treatment of the viscous terms for incompressible flows is proposed in [START_REF] Kim | Application of a fractional-step method to incompressible Navier-Stokes equations[END_REF]. This is based on an approximate factorization which requires solving tridiagonal matrices rather than inversion of a large sparse matrix. This alternative method has not been tested in this work.

At the opposite, as it has been previously shown, for incompressible flows the second-order cross-derivatives are perfectly compensated everywhere except at the interface where viscosity is discontinuous. Therefore, the GFSCM can be used with benefits with an implicit treatment of the viscous terms, which leads to the following discretization:

ρ n+1 i+1/2 j u * i+1/2 j -t D x µD x u * + D y µD y u * i+1/2 j = ρ n+1 i+1/2 j u n i+1/2 j -t V n • ∇ u n | i+1/2 j , (87)
ρ n+1 i j+1/2 v * i j+1/2 -t D x µD x v * + D y µD y v * i j+1/2 = ρ n+1 i j+1/2 v n i j+1/2 -t V n • ∇ v n | i j+1/2 . ( 88
)
The resulting linear systems can be solved as easily as in the single-phase case: the matrices of discretization are symmetric, definite and positive, which allows the use of an efficient and simple black-box solver. Note that an implicit temporal discretization of the viscous terms can also be used with the GFPM [START_REF] Hong | On boundary condition capturing for multiphase interfaces[END_REF].

We conclude by some comments on the way to define the density: ρ n+1 i+1/2 j and ρ n+1 i j+1/2 . The definition must be compatible with the computation of the correction step, which itself relies on the manner density is computed when the Poisson equation on pressure is solved. Accurate results are obtained by using the following general expression for ρ n+1 i+1/2 j and ρ n+1 i j+1/2 in the grid cells that are crossed by an interface:

1 ρ n+1 = 1 ρ + ρ - 1 ρ + θ + 1 ρ -(1 -θ) = 1 ρ -θ + ρ + (1 -θ)
.

(89)

Summary about the four different methods to deal with viscosity

Table 1 gives an overview of the characteristics of the four different methods to deal with viscous terms in the simulations of two-phase flows with Level-Set methods, and describes the main steps for implementing it.

Level-Set methods for interface capturing

Interface motion is captured with a Level-Set method [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 flow[END_REF], which consists in solving a convection equation for the Level-Set function φ:

∂φ ∂t + ( V • ∇)φ = 0, (90)
with V the velocity at the interface. A reinitialization step [START_REF] Sussman | A Level Set approach for computing solutions to incompressible two-phase flow[END_REF] can next be performed to insure that the φ function remains a signed distance in the computation domain. This can be done by solving iteratively the specific PDE:

∂d ∂τ = sign(φ) 1 -|∇d| , ( 91
)
where d is the reinitialized distance function, τ a fictitious reinitialization time and sign(φ) a smoothed signed function defined in [START_REF] Fedkiw | A non-oscillatory Eulerian approach to interfaces in multimaterial flows (The Ghost Fluid Method)[END_REF][START_REF] Sussman | A Level Set approach for computing solutions to incompressible two-phase flow[END_REF]. At the end of every physical time step, two temporal iterations of the redistancing equation (91) are solved using a second-order TVD-Runge-Kutta scheme.

The signed distance function allows preserving a good accuracy when computing interface geometrical properties as the normal vector or the curvature with the following simple relations:

n = ∇φ |∇φ| , (92)
κ(φ) = -∇ • n.
(93)

V * 2 = V n -t ( V n • ∇) V n -∇•(2µD n ) ρ n+1 -σ κ nδΓ ρ n+1 -g (33)
One complex coupled parabolic linear system for all the velocity components

∇ • ∇ p n+1 ρ n+1 = ∇• V * 2 t (34) V n+1 = V * 2 -t ∇ p n+1 ρ n+1 (35)
"Jump condition formulation"

Ghost-Fluid Primitive viscous Method GFPM Kang et al. [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF] Hong et al. [START_REF] Hong | On boundary condition capturing for multiphase interfaces[END_REF] V *

1 = V n -t ( V n • ∇) V n -µ V n ρ n+1 -g (30)
Several simple independent linear systems with explicit jump conditions on the components of the viscous tensor

∇ • ∇ p n+1 ρ n+1 = ∇• V * 1 t + ∇ • (σ κ+2[µ]Γ ∂ Vn ∂n ) nδΓ ρ n+1 (36) ⇔    ∇ • ∇ p n+1 ρ n+1 = ∇• V * 1 t [p n+1 ] Γ = σ κ + 2[µ] Γ ∂ Vn ∂n V n+1 = V * 1 -t ρ n+1 ∇ p n+1 -σ κ + 2[µ] Γ ∂ Vn ∂n nδ Γ (39) Ghost-Fluid Conservative viscous Method GFCM Sussman et al. [26] V * 3 = V n -t ( V n • ∇) V n -∇•(2µD n ) ρ n+1 -g (40)
One complex coupled parabolic linear system for all the velocity components

∇ • ∇ p n+1 ρ n+1 = ∇• V * 3 t + ∇ • σ κ nδΓ ρ n+1 (42) ⇔ ∇ • ∇ p n+1 ρ n+1 = ∇• V * 3 t [p n+1 ] Γ = σ κ V n+1 = V * 3 -t ρ n+1 (∇ p n+1 -σ κ nδ Γ ) (47) Ghost-Fluid Semi-Conservative viscous Method GFSCM Present paper V * 4 = V n -t ( V n • ∇) V n -∇•(µ∇ V n ) ρ n+1 -g (41)
Several simple independent linear systems without any jump conditions on the components of the viscous tensor

∇ • ∇ p n+1 ρ n+1 = ∇• V * 4 t + ∇ • (σ κ+[µ]Γ ∂ Vn ∂n ) nδΓ ρ n+1 (43) ⇔    ∇ • ∇ p n+1 ρ n+1 = ∇• V * 4 t [p n+1 ] Γ = σ κ + [µ] Γ ∂ Vn ∂n V n+1 = V * 4 -t ρ n+1 ∇ p n+1 -σ κ + [µ] Γ ∂ Vn ∂n nδ Γ (48)
The curvature is firstly computed on all the grid nodes with simple finite difference schemes. Next, if the Ghost Fluid Method is used, the curvature is interpolated on the interface. For example, if it crosses the segment [x i , x i+1 ] the following interpolation is used:

κ I = κ i+1 θ + κ i (1 -θ) (94)

Temporal and spatial discretizations

All convective terms and spatial derivatives of the redistance equation are discretized by means of a fifth-order WENO 5 scheme [START_REF] Jiang | Efficient implementation of weighted essentially non-oscillatory schemes[END_REF]. First-order explicit or semi-implicit temporal integration methods have been described in previous sections. These elementary steps are finally included in a second-order TVD-Runge-Kutta scheme for the overall temporal integration of the PDE system in order to improve the computation stability. A classical time-step constraint [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF][START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF][START_REF] Sussman | A Level Set approach for computing solutions to incompressible two-phase flow[END_REF] accounting for the stability conditions on convection, viscosity (in case of an explicit temporal discretization) and surface tension is imposed.

t conv = x max V (95) t visc = 1 4
x 2 max(υ liq , υ gaz )

(96)

t surf _tens = 1 2 ρ liq x 3 σ (97)
If an implicit temporal discretization is used with the GFSCM, an explicit part depending on the viscosity jump condition remains in Eqs. (43), (48). Therefore, the corresponding term involves another temporal constraint which is different from Eq. ( 96), since it depends on a viscosity jump condition instead of a maximum kinematic viscosity value. By performing several numerical test-cases with zero surface tension and very low gravity (in order to remove the stability constraints due to surface tension and convection), we have observed that the following time step restriction is sufficient to ensure the numerical stability:

t visc = ρ x 2 [µ] (98)
with ρ the average density of the two phases. Whereas this constraint depends also on x 2 , it can be much more lenient than the classical one, Eq. ( 96), because ρ

[µ] ∼ 1 υ liq is much higher than 1 υ gas for many typical gas-liquid systems such as air-water or the one used in Section 3.3. For example, it is shown in Section 3.3 that the implicit GFSCM allows performing the computation with a much larger time step in comparison to its explicit counterparts. Finally, the global time step restriction can be computed with the following relation:

1 t = 1 t conv + 1 t visc + 1 t surf _tens (99)

Numerical results

Relevant benchmarks have been designed to assess and compare the four different numerical methods discussed in this paper: DFM, GFPM, GFCM and GFSCM.

In the first benchmark, computations of a spherical bubble rising under gravity are performed for a range of Reynolds number varying from 20 to 100.

The second benchmark consists of numerical simulations of bubble shape-oscillations in the absence of gravity. Such a physical problem is characterized by an oscillation Reynolds number which compares the oscillations frequency to the damping rate. The accuracy of the simulations is evaluated from comparisons with the theoretical prediction of the frequency and the damping rate of the oscillations for Reynolds numbers varying from 20 to 300. In addition, a study of the spatial convergence of the simulations is carried out in order to determine how mesh-grid requirements evolve when the Reynolds number is varied.

Finally, the third benchmark focuses on the rise of a non-spherical bubble in a very viscous liquid as in Refs. [START_REF] Bonometti | Transition from spherical cap to toroidal bubbles[END_REF][START_REF] Gueyffier | Volume-of-fluid interface tracking with smoothed surface stress methods for three-dimensional flows[END_REF][START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF]. These three benchmarks consist of axisymmetric computations.

A spherical bubble rising under gravity

Currently, performing well-resolved simulations of bubbly flows or droplet sprays is a challenge for the direct numerical simulation of two-phase flows. Many issues must be overcome to achieve such ambitious targets. One can expect that the development of High Performance Computing hardware and software will allow performing simulations involving a growing number of bubbles or drops. However well-resolved 3D simulations remain costly and a particular attention must be paid to the development and the assessment of accurate numerical methods to solve the fluid dynamics equations. In particular, a computation which is efficient to describe accurately the dynamics with half the number of grid cells on a bubble diameter makes possible to compute eight times more bubbles with the same computational resources. Actually, the gain would even be significantly greater since, for explicit methods, a larger time step can be used on a coarser grid. The test-case proposed in this section allows us to assess the different discretizations of the viscous terms. Simulations are performed with an axisymmetric coordinate system. The radial and axial sizes of the computational domain are defined relatively to the bubble radius R bubble : l r = 8R bubble and l z = 4l r . Considering a spherical bubble, reference results on the drag coefficient are provided by boundary-fitted numerical simulations from [START_REF] Oliver | Flow about a fluid sphere at low and moderate Reynolds numbers[END_REF] when Re ∞ ≤ 50 or by Moore's theory [18] for higher Reynolds numbers.

Our computations are performed for low Bond numbers to preserve bubble sphericity. It is worth mentioning that the case of a rising spherical bubble is particularly advised to assess the discretization of the viscous terms because, in this situation, the drag force mainly depends on viscous dissipation and not on pressure effects, unlike the case of a strongly deformed bubble.

The physical properties of liquid and gas are kept identical to water and air for all simulations presented in this section:

ρ liq = 1000 kg m -3 , µ liq = 0.00113 kg m -1 s -1 , ρ gas = 1.226 kg m -3 , µ gas = 1.78 × 10 -5 kg m -1 s -1 . Simulation conditions are summarized in Table 2 for five configurations corresponding to Re ∞ = 20, 40, 60 for a Bond number Bo = 0.025, and to Re ∞ = 80, 100 for Bo = 0.0125. Bubbles are initially static and evolve under gravity until the terminal velocity is reached.

The boundaries of the computational domain are defined as walls (except on symmetry axis) and set sufficiently far from the bubble to have a negligible influence on the results.

Simulations are performed with three different mesh grids (32 × 128, 64 × 256, 128 × 512, 256 × 1024), which describe respectively one bubble diameter with 8, 16, 32 or 64 mesh points. Fig. 2 shows a snapshot of the interface location and the fully-developed vorticity in the frame moving with the bubble for the case Re ∞ = 60 with the thinnest grid. Results of all the simulations are summarized in Tables 3, 4, 5, 6 and 7, which correspond respectively to computations carried out using the Ghost Fluid Conservative viscous Method (GFCM), explicit and implicit Ghost-Fluid Semi-Conservative viscous Method (GFSCM), Ghost-Fluid Primitive viscous Method (GFPM) and Delta-Function Method (DFM). Many conclusions can be drawn from these tables. First, all methods fail to provide accurate results on the coarsest grid regardless the value of Re ∞ : this mesh is clearly insufficient to capture accurately viscous effects. When the intermediate grid is used, the GFCM and the explicit and implicit GFSCM provides very good results: whatever Re ∞ , deviations between reference results and numerical simulations are less or equal to 3% (except for Re ∞ = 80 when the explicit GFSCM is used). The GFPM also provides acceptable results on this

Table 3

Error on the terminal Reynolds number of rising bubbles for different computational grids with the explicit GFCM.

Grids

Re ∞ =

Fig. 1 .

 1 Fig. 1. Sketch of a two-phase flow domain Ω 1 U Ω 2 with an immersed interface Γ .

Fig. 2 .

 2 Fig. 2. Streamlines, interface location and vorticity field (s -1 ) in the moving frame of a rising spherical bubble at Re ∞ = 60 (grid 128 × 512) using the GFCM.

Table 1

 1 Summary of four different methods to deal with viscous terms with Level-Set methods.

	Framework	Method and reference papers	Steps of the explicit projection method and number of the	Implementation of an
			corresponding equations in the reference paper	implicit method
	"Whole domain	Delta-Function Method
	formulation"	DFM
		Sussman et al. [3,24,25]

Table 2

 2 Bubble radius, gravity acceleration, surface tension and final time of simulation for a spherical bubble rising at various Reynolds numbers. Re ∞ = 20 Re ∞ = 40 Re ∞ = 60 Re ∞ = 80 Re ∞ = 100

	R bubble (µm) g z (m s -2 ) σ (N m -1 )	19.55 -9144.4 0.56	44.65 -1755.2 0.56	70.45 -705.6 0.56	98.95 -357.5 1.12	127.9 -213.9 1.12
	t f (s)	0.0004	0.00135	0.002375	0.0048	0.007

  20 Re ∞ = 40 Re ∞ = 60 Re ∞ = 80 Re ∞ = 100

	32 × 128	-22.0%	-34.0%	-38.8%	-63.0%	-67.4%
	64 × 256	+2.5%	+1.8%	-1.3%	-1.3%	-3.0%
	128 × 51	+0.5%	+0.8%	-1.0%	+1.9%	+1.9%
	256 × 1024	-0.9%	-0.4%	-1.95%	-1.02%	-0.65%

Table 4

 4 Error on the terminal Reynolds number of rising bubbles for different computational grids with the explicit GFSCM.

	Grids	Re ∞ = 20	Re ∞ = 40	Re ∞ = 60	Re ∞ = 80	Re ∞ = 100
	32 × 128	-23.5%	-32.5%	-39.0%	-62.5%	-67.0%
	64 × 256	-0.5%	+2.5%	+1.3%	+4.3%	+3.0%
	128 × 512	-3.5%	+0.6%	0.7%	+5.4%	+6%
	256 × 1024	-2.15%	+3.25%	+4.2%	+6.75%	+8.2%

Table 5

 5 Error on the terminal Reynolds number of rising bubbles for different computational grids with the implicit GFSCM.

	Grids	Re ∞ = 20	Re ∞ = 40	Re ∞ = 60	Re ∞ = 80	Re ∞ = 100
	32 × 128	-18.6%	-25.0%	-30.2%	-52.9%	-55.1%
	64 × 256	-1.2%

Acknowledgements

The authors gratefully acknowledge funding by ANR (French Research National Agency) in the frame of the IDHEAS (Interaction of Droplets with HEAted Surface) program.

grid since the discrepancy between numerical simulations and reference solutions is less than 7%; the error increases with the Reynolds number, which is probably due to the difficulty to calculate accurately boundary layers that becomes thinner at high Re ∞ . By contrast, the DFM provides a lower accuracy than the other methods with errors on Re ∞ which vary from 15% to 30% and increase with Re ∞ . Finally, considering results obtained with the most refined grids, it clearly appears that both the GFCM and the GFPM have converged to the correct solution. The explicit and implicit GFSCM provide results with larger errors than the two previous methods with errors varying between 2% and 8%. The DFM seems to converge slowly towards the theoretical solution and still shows deviations varying from 4% to 6% (according to Re ∞ ) with 64 mesh points on a bubble diameter.

Shape oscillations of a bubble

Considering small-amplitude-shape-oscillations of a bubble in the absence of gravity, theoretical solutions of the linearized Navier-Stokes equations for an incompressible flow exist and give expressions for the frequency and the damping rate of the oscillations [START_REF] Miller | The oscillations of a fluid droplet immersed in another fluid[END_REF][START_REF] Prosperetti | Normal mode analysis for the oscillations of a viscous liquid drop in an immiscible liquid[END_REF]. We have took advantage of this to design an accurate and stringent test-case for direct numerical simulations. While the oscillation frequency depends essentially on capillary effects, we have focused on the damping rate in order to test the accuracy of different numerical approaches about viscosity. Like in the previous section, simulations are performed with an axisymmetric coordinate frame. The dimensions of the computational field are: l r = 4R bubble , l z = 2l r . Wall boundary conditions are used on all domain borders (except on symmetry axis) and confinement effects have been checked to be negligible. Densities and viscosities, corresponding to an air/water system, are the same as in Section 3.2, and the surface tension is σ = 0.07 N m -1 . We propose a parametric study by varying the Reynolds number of oscillation, by Re osc = σ R bubble ρ liq µ liq , from 20 to 300. (Note that the Reynolds number of oscillation is the inverse of the Ohnesorge number.) In our simulations, the initial condition corresponds to a deformed interface described by a Legendre polynomial of order 2 with the following Level Set function

where (r, θ) are the spherical coordinates, P 2 (cos θ) = 1 2 (3 cos 2 θ -1) and ǫ is the dimensionless initial amplitude of the deformation, chosen low enough (ε = 0.075) to maintain the oscillations in the linear regime. The computed shape of the interface is extracted at each time step and decomposed into spherical harmonics in order to obtain the time evolution of amplitude a 2 (t) of harmonic l = 2: r(θ, t) = R bubble + a 2 (t)P 2 (cos θ).

(101) Then, the frequency ω 2 and the damping rate β 2 of the oscillations are obtained from a 2 (t) since a 2 (t) = εR bubble cos(ω 2 t)e -β 2 t (102)

The values issued from the simulations are compared with the theoretical values ω th 2 and β th 2 that are computed by solving numerically the characteristic equation given by Prosperetti (eq. (33) of Ref. [START_REF] Prosperetti | Normal mode analysis for the oscillations of a viscous liquid drop in an immiscible liquid[END_REF]). The damping rate is a global parameter that requires an accurate calculation of the viscous terms all over the computation domain.

Snapshots of the interface location and vorticity are shown in Fig. 3 and Fig. 4 for Re osc = 20 and Re osc = 200, respectively. It is clear from these figures that the vorticity is localized inside the bubble, whereas the flow remains potential outside. For a shape-oscillating bubble, both the frequency and the damping rate are mainly controlled by the potential flow in the liquid phase. The asymptotic development of the damping rate given by Miller and Scriven [START_REF] Miller | The oscillations of a fluid droplet immersed in another fluid[END_REF] shows that, for a bubble, the viscous dissipation in the boundary layers is responsible for only 1 to 3% of the total dissipation in a present

Table 8

Error on the damping rate β 2 for an oscillating bubble with different computational grids with the explicit GFCM. , where the deviation of the numerical simulations relative to the theoretical predictions have been reported. The error on the frequency of oscillation is very small (<1%) whatever the numerical method and the mesh grids; it has therefore not been reported. A detailed analysis of the errors on the damping rate leads to the following conclusions. First, the GFCM and the GFPM accurately match with the exact solution when refined grids are used. Both methods present comparable abilities on this test-case, with errors less than 2% in the range 20 ≤ Re osc ≤ 200 (the GFPM is just slightly better than the GFCM). For the highest Reynolds numbers of oscillation considered, the error seems to decrease monotonously and the order of convergence can be roughly estimated to 1.5. That is not surprising considering that, even if the discretization of viscous terms remains of the second order far from the interface, the subcell viscosity interpolation in the cells that are crossed by an interface is only of the first order. The DFM also provides good results in this case and converges well, but exhibits a different behavior than the GFCM or the GFPM. For example, when coarse grids are used, the error is relatively low (between 5 and 10%) whatever the Re osc , contrary to the two previous methods where the error was clearly increasing with Re osc to become larger than 20% at the highest Re osc . The rate of convergence of the DFM is not as clear as for the GFCM or the GFPM. Finally, the explicit GFSCM provides poor results with errors up to 12% in some cases with the thinnest grid. Although this discrepancy is not negligible, a more detailed analysis of the temporal signal a 2 (t), presented in Fig. 5 for the most disadvantageous case for the GFSCM (Re osc = 100, grid 256 × 512), allows us to realize that this underestimation of 12% of the damping rate does not induce a so significant difference compared to the other methods.

In addition, it is worth mentioning that the GFCM has been successfully used in Ref. [START_REF] Lalanne | Effect of rising motion on the damped shape oscillations of drops and bubbles[END_REF] to simulate the shape oscillations of rising drops and bubbles in order to determine the effect of the rising motion on the frequency and the damping rate of the oscillations.

Deformed rising bubble at Re

To complete the assessment of the different numerical methods, we propose a third test-case which consists in computing the rise of a bubble in a very viscous liquid. This test-case is based on experimental data [START_REF] Hnat | Spherical cap bubbles and skirt formation[END_REF] and has already been used for the validation of numerical simulations, see Refs. [START_REF] Bonometti | Transition from spherical cap to toroidal bubbles[END_REF][START_REF] Gueyffier | Volume-of-fluid interface tracking with smoothed surface stress methods for three-dimensional flows[END_REF][START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF] for example. Its interest lies in the strong deformation of the bubble, which involves the development of a recirculation zone in the bubble wake. Here, unlike the case of a spherical rising bubble, the drag force results both from pressure and viscous effects. Simulations are performed in axisymmetric coordinates and the computational domain is defined by l r = 8R bubble and l z = 4l r . We used the following physical properties for the fluids: ρ liq = 875.5 kg m -3 , µ liq = 0.118 kg m -1 s -1 , ρ gas = 1 kg m -3 , µ gas = 0.001 kg m -1 s -1 , σ = 0.0332 N m -1 ,

All results obtained for this test-case are summarized in Table 12, the simulations being carried out with a grid (128 × 512), which counts 32 points in a bubble diameter. We compare the GFCM, GFPM, DFM and implicit GFSCM. A snapshot of interface location and streamlines around and inside the bubble are plotted in Fig. 6; it clearly shows the existence of a recirculation zone in the bubble wake. Temporal evolutions of the Reynolds number have been plotted in Fig. 7 to describe the results obtained during the transient stage. The reference terminal rising velocity issued from experiments corresponds to Re ∞ = 19.6. We observe that all methods provide satisfactory results with errors on the Reynolds number varying from 4.1% for the GFCM to 6.9% for the GFSCM. We can conclude that all tested methods have a comparable accuracy for this test-case. Fig. 8 highlights the interest of the implicit GFSCM on a simulation where the time step restriction due to viscosity (96) is dominant. This time step restriction can be removed if the implicit GFSCM is used. Therefore the numerical stability is maintained by imposing the time step restriction due to capillary forces (97) and the modified time step restriction due to the viscosity jump (98) which have both comparable values in this case.

The cost of a temporal iteration is approximately the same for an explicit discretization or the implicit GFSCM. Indeed, the linear systems resulting from the implicit discretization of the viscous terms are diagonally dominant due to the other time step restrictions. Therefore, their solutions are computed much faster than for the pressure Poisson equation. Thus, for the simulations which are considered in Fig. 8 the implicit GFSCM will run 10 to 20 faster than an explicit version of the GFSCM, GFCM, GFPM or DFM.

Conclusions

In this paper, we have proposed an extensive study to assess different numerical methods for the computation of the viscous terms in the framework of direct numerical simulation of two-phase flows with Level-Set methods. An important feature of this work is to provide the theoretical evidence that the Delta-Function Method described by Sussman et al. in [START_REF] Sussman | A Level Set approach for computing solutions to incompressible two-phase flow[END_REF], the original Ghost-Fluid Method for incompressible two-phase flows proposed by Kang et al. in [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF], and another version of the Ghost-Fluid Method proposed later by Sussman et al. in [START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF] are formally equivalent. For sake of clarification, these three methods are respectively renamed in this paper the Delta-Function Method (DFM), the Ghost-Fluid Primitive viscous Method (GFPM), and the Ghost-Fluid Conservative viscous Method (GFCM). It is shown that all these different methods result from different choices of the predictive velocity V * , but they all lead to solve the same set of formal equations. Thanks to this proof, we have designed a novel method to compute viscous effects in two-phase flows, which is named the Ghost-Fluid Semi-Conservative viscous Method (GFSCM). The main interest of this new method lies in its ease of implementation for the implicit temporal integration of the viscous terms.

Next, the respective accuracy of these different methods has been compared from many simulations based on three relevant test-cases: a spherical bubble rising at different Reynolds numbers, damped shape-oscillations of a bubble, and the rise of a non-spherical bubble at low Reynolds number. Important conclusions have been drawn from these benchmarks. The GFCM and the GFPM succeed all tests and clearly appear as the most accurate methods, with comparable results. The DFM presents a poor efficiency in the case of a spherical rising bubble because of its slow spatial convergence, but works well for the other tests. Finally, the GFSCM provides results in simulations with rising bubbles but less satisfactory results for bubble shape oscillations, with a precision on the damping rate of about 10%. Nevertheless, thanks to its simplicity to compute an implicit discretization of the viscous terms and its ability to treat the case of discontinuous velocities at the interface (due for example to mass transfers or phase changes), this approach is promising to deal with challenging complex problems for which slight inaccuracies can be tolerated.