Daniela Garajeu

Marc Medale

An XF EM model for incompressible two-fluid flows with arbitrary high contrasts in material properties

Keywords: Two-fluid flows, strong and weak field discontinuities, Extended Finite Element Method

published or not. The documents may come

An XFEM model for incompressible two-fluid flows with arbitrary high contrasts in material properties

Daniela Garajeu, Marc

Introduction

Two-fluid flows designate fluid flows of two immiscible fluids, which are therefore separated by a sharp interface whose position is usually part of the problem unknowns to be determined. They take place in a large collection of fluid mechanics problems ranging from environmental and geophysical to industrial processing situations.

Many numerical methods have been developed to deal with such problems, but their accuracy and computational efficiency mainly depend on their capabilities to deal with two key-points: i) the topological complexity of the interface(s) and its (their) potential changes in the course of time; ii) the material discontinuities across the interface(s) and their related velocity and pressure field discontinuities. Depending on the two-fluid flow problem to be solved two kinematic descriptions are usually used, each one having its own virtues and drawbacks. For sharp material contrasts between the two fluids but quite simple interface topology the mixed Lagrangian-Eulerian description or Arbitrary Lagrangian-Eulerian (ALE) description is well suited [START_REF] Hirt | An arbritary lagrangian-eulerian computing method for all flow speeds[END_REF][START_REF] Hughes | Lagrangian-eulerian finite element formulation for incompressible viscous flows[END_REF], provided the interface topology does not change in the course of time. This approach requires

Ω 2 Ω 1 Γ(t) ρ2 μ2 ρ1 μ1
(a) Interface fitting mesh. to build a mesh such that some element edges coincide with the interface (no crossed elements) dealing with the material discontinuity in a direct way, cf. figure 1(a). For fixed interface problems this approach provides accurate solutions at low computational cost. However, for moving interface problems a repositioning of mesh nodes is necessary at each time-step and domain re-meshing could be required to sustain accuracy and could become computationally expensive especially in three-dimensional problems. So, in these cases or in those where it is hard to build a mesh that coincides with the interface(s), i.e. multiple interfaces, the ALE approach is no longer suitable and an Eulerian description is more appropriate, cf. figure 1(b). Therefore, the two-fluid flow is now computed on a fixed mesh and the interface can move over it and does not coincide any more with any mesh lines or edges. The interface location and movement are described according to Lagrangian (front tracking) or Eulerian (front capturing) methods. In front or interface tracking methods a set of discrete particles are set along the interface and are advected by the local flow, as in the pioneer "Marker and Cells" method [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluids with a free surface[END_REF]. Accordingly, an interface grid or mesh can also be built instead of discrete particles to better compute high curvature interfaces [START_REF] Univerdi | A front-tracking method for viscous, incompressible, multi-fluid flows[END_REF][START_REF] Esmaeeli | Direct numerical simulations of bubbly flows. Part I: Low Reynolds number arrays[END_REF] when surface tension enters the game. These methods are easy to implement, but to preserve accuracy in cases of strong flow heterogeneities markers or mesh points deposited along the interface have to be redistributed leading here again to high computational costs in 3D. Moreover, front tracking methods are not well suited to handle topology changes of the interface(s) (coalescence of two interfaces, interface separations, etc.), so one has to consider in these latter cases interface capturing methods, in which the interface is implicitly represented by scalar field defined in the whole computational domain [START_REF] Hirt | Volume of fluids (vof) method for the dynamic of free boundaries[END_REF]. The interface is associated with a particular value of the scalar field (level set), which is advected by the flow in an Eulerian framework [START_REF] Osher | Front propagating with curvature dependent speed: Algorithms based on hamilton-jacobi formulations[END_REF][START_REF] Thompson | Use of pseudo-concentrations to follow creeping viscous flows during transient analysis[END_REF]. Interface capturing methods are basically Eulerian-Eulerian ones, so they are nicely well adapted for problems with several complex interfaces, in which topology changes are implicitly undertaken. But, mass conservation is a recurrent issue of this class of methods.

Let us now review the different ways to deal with material discontinuities along with their related velocity and pressure field ones in incompressible twofluid flows. Two kinds of models can be distinguished in Eulerian approaches: continuous one-fluid models and discontinuous two-fluid ones. The basic idea behind the former approach is to consider only one fluid over the whole domain, whose material properties depend on the fluid volume fraction in every cell or mesh element. So, in cells where the two fluids are simultaneously present average values (algebraic or harmonic) of each fluid properties are computed.

Therefore, material and fields discontinuities across the interface are modeled as continuously varying fields whose steepness depends on mesh resolution. This approach is very simple and computationally efficient as far as the mesh size in transition elements is small enough to properly account for a continuous spatial representation of averaged material properties. It is noteworthy that this condition is related to the contrast in material properties: the higher this contrast the finer the mesh size in transition elements. So, when this condition involves unaffordable mesh refinement in 3D computations, one has to move to two-fluid models that account for discontinuous material properties, along with the related discontinuities they induce in velocity and pressure fields. A class of methods have emerged in the Finite Element Method (F EM ) to supplement the classical continuous polynomial approximations with extra discontinuous functions that satisfy the local partition of unity [START_REF] Melenk | The partition of unity finite element method: Basic theory and applications[END_REF] and account for the convenient discontinuities across the interface. These methods stand as the Extended Finite Element Method (XF EM ), they have been designed to approximate discontinuous fields thanks to enrichment functions directly related to the known physical phenomenon inducing the discontinuity [START_REF] Chessa | An extended finite element method for two-phase fluids[END_REF][START_REF] Moes | Discontinous enrichment in finite elements with a partition of unity method[END_REF][START_REF] Fries | The extended/generalized finite element method: An overview of the method and its applications[END_REF]. This is the way we have followed in this paper to deal with arbitrary contrasts in material properties arising in two-fluid flow problems, however to reach this goal one still has to overcome the issue related to the computation of domain integrals within mesh elements crossed by the interface. Indeed, the presence of two material properties in the same element leads to discontinuous integrands, that could be strongly discontinuous for fluids with high contrasts in material properties. In this case, Gauss quadrature over the entire two-fluid elements no longer produces accurate results, therefore preliminary splitting of these twofluid elements enables to achieve numerical integration over two homogeneous sub-domains where integrands are piecewise continuous.

The paper is organized as follows: section 2 presents the eulerian model for incompressible two-fluid flows. It relies on the Navier-Stokes equations on each fluid supplemented with the Young-Laplace equation along the interface, modeled by the Level-Set method. Section 3 describes the spatial discretization of the weak form performed by standard finite element method for the continuous part of the fields and with XF EM enrichments for the discontinuous one. The main numerical integration techniques enabling to properly deal with elements crossed by the interface are presented in section 4. Several applications are presented in section 5 to validate velocity and pressure field enrichments for three-dimensional problems. Finally, the paper closes with a discussion that emphasizes the main difficulties encountered with the XF EM approach, the solutions we have adopted to overcome them and few issues that still seem to remain open questions beyond this work.

An Eulerian model for two-fluid flows

Let us consider a domain Ω in which takes place the incompressible fluid flow of two immiscible fluids separated by an interface Γ, that defines two complementary material subdomains Ω 1 and Ω 2 of Ω (cf. fig. 1(b)). The interface between the two fluids locates surface tension, denoted τ , along with the related discontinuity of pressure field (if the interface has a nonzero curvature), as well as a discontinuity of the velocity gradient (if µ 1 = µ 2 or ∂τ ∂s = 0, where s stands for the curvilinear abscissa along the interface).

Governing equations of fluid flow

In each subdomain Ω α (α = 1, 2) the fluid flow is governed by the incompressible Navier-Stokes equations. Mass conservation equation and momentum read for two-fluid flows in moving interface problems described in an eulerian framework as follows:

∂ρ α ∂t + ∇ • (ρ α v α ) = 0. (1)
ρ α ∂ v α ∂t + v α • ∇ v α = ∇ • σα + ρ α g, (2)
where t denotes time, ρ α and µ α are density and dynamic viscosity of fluid α, assumed constants in each subdomain Ω α , g is the gravitational acceleration and σα is the Cauchy stress tensor defined by a Newtonian behavior, p α and v α are pressure and velocity vector, respectively:

σα = -p α Ī + 2 µ α εα , with εα = 1 2 ∇ v α + ∇ T v α .
The closures at the interface are two-fold: on the one hand velocity continuity is assumed for immiscible fluids (at least one liquid) and on the other hand the Young-Laplace equation provides the jump of stress vector across the interface induced by surface tension:

[ v] Γ(t) = 0 [ σ n] Γ(t) = τ χ n + ∂τ ∂s t g , (3)
where χ is the interface curvature, n is the normal vector to the interface (arbitrarily oriented from Ω 1 to Ω 2 ) and t g is the tangent vector along the interface.

Finally, initial and boundary conditions, specific to each problem must be supplemented to close the problem.

Level Set description of the interface

The interface Γ between the two fluids is an an implicit curve (in 2D) or surface (in 3D), defined as the zero-level of a level set function:

Φ( x, t) = 0,
where Φ( x, t) is the signed distance to the interface (figure 2). The domain Ω is divided into two subdomains Ω 1 and Ω 2 by the interface where the level set function changes sign, arbitrarily defined positive in the former and negative in the latter. The scalar function Φ is continuous in space and can be time dependent. If so, its time evolution is modeled by the pure advection equation of the Level Set method [START_REF] Osher | Level set methods: an overview and some recent results[END_REF][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], which requires as any hyperbolic equations a stabilization technique to be solved by classical Bubnov-Galerkin Finite Element approximations:

n 2 1 t 0 0 0 Φ( x, t) = 0 on Γ Φ( x, t) < 0 on Ω 1 Φ( x, t) > 0 on Ω 2
∂Φ( x, t) ∂t + v α • ∇Φ( x, t) = 0. (4)

Weak integral forms and fields approximations

The specificity of incompressible fluid flows stands in its peculiar velocitypressure coupling, which is undertaken in the present work by the unconditionally stable projection algorithm [START_REF] Guermond | Calculation of incompressible viscous flows by an unconditionally stable projection fem[END_REF], later on extended to free-surface flows and open boundary conditions [START_REF] Rabier | Computation of free surface flows with a projection fem in a moving mesh framework[END_REF]. The resulting algorithm consists in solving iteratively the momentum and mass conservation equations in two steps: i) an advection-diffusion step, in which the velocity field is advanced in time by solving for the momentum equation with a provisional pressure field; ii) a projection step, in which the pressure field is updated so that the velocity field will satisfy the incompressibility constraint at the end of the step.

Weak integral forms

The Finite Element Method is used to transform the governing equations into algebraic systems tractable on modern parallel computers. Its first step consists in building weak integral form(s) of the problem, that is then discretized in space by finite element approximations (classical and XF EM ) and in time by first order finite difference scheme (backward Euler). The unconditionally stable projection algorithm, adapted to present two-fluid flow problems leads to the following three weak integral forms on material domains Ω α (α = 1, 2):

Ωα ρ α v * -v t ∆t • δ v dΩ + Ωα ρ α v * ∇ v * • δ v dΩ + 2 Ωα µ α ε * : δ ε dΩ -Ωα ∇(2p t -p t-∆t ) δ v dΩ + ∂Ωα σ n • δ v dΣ + Ωα ρ α g e z • δ v dΩ = 0 Ωα ∇(p t+∆t -p t ) • ∇δp dΩ + 1 ∆t Ωα ρ α ( ∇ • v * ) δp dΩ = 0 (5) Ωα Φ t+∆t -Φ t ∆t δΦ dΩ + Ωα v • ∇Φ δΦ dΩ = 0 ∀ ( v * , v t ) ∈ U v , ∀ δ v ∈ W v , ∀ p ∈ U p , ∀ δp ∈ W p , ∀ Φ ∈ U Φ , ∀ δΦ ∈ W Φ
where U and W are appropriate functional spaces for trial and weighting functions, ∆t is the time step and starred quantities are based one the intermediate velocity vector ( v * ), which is not required to satisfy the incompressibility constraint in the advection-diffusion step.

Spatial discretization of weak integral forms

The weak integral forms of eq. ( 5) are continuous in space, so one has next to approximate all variable fields in order to build algebraic systems to be solved on a computer. The physical domain Ω and its boundary ∂Ω are approximated by a mesh, which is not aligned with the interface Γ, so that some mesh elements are crossed by it. Therefore, one uses two different approximation spaces depending on the location of the mesh element with respect to the interface: the standard finite element basis for the continuous part of the fields and XF EM enrichments for the discontinuous one (mesh elements crossed by the interface).

F EM approximations in single-fluid elements

For mesh elements that are not crossed by the interface (single-fluid elements), velocity, pressure and Level Set fields are approximated with standard Finite Element basis functions (Bubnov-Galerkin aproximations), using iso-parametric Lagrange elements: Q 9 and Q 4 quadrilateral elements in two dimensions; H 27 and H 8 hexahedral elements in three dimensions. Thus, the field approximations read:

• piecewise quadratic approximations for the velocity field and its weighting function:

v( x, t) = 9 i=1 N Q9 i ( x) V i (t), in 2D and v( x, t) = 27 i=1 N H27 i ( x) V i (t), in 3D
• piecewise linear approximations for the pressure field and its weighting function:

p( x, t) = 4 i=1 N Q4 i ( x) P i (t), in 2D and p( x, t) = 8 i=1 N H8 i ( x) P i (t), in 3D
• piecewise quadratic approximations for the Level Set field and its weighting function:

Φ( x, t) = 9 i=1 N Q9 i ( x) Φ i (t), in 2D and Φ( x, t) = 27 i=1 N H27 i ( x) Φ i (t), in 3D (6)

XF EM approximations in two-fluid elements

On the other hand, mesh elements crossed by the interface (two-fluid elements) are made-up of two subdomains with different material properties, which induce discontinuities across the interface in velocity and pressure fields or in their gradients, as follows:

• the velocity field is continuous across Γ, but its normal gradient is discontinuous across the interface if it exists any contrast in material properties between fluids, cf. figure 3(a). This is referred to as weak discontinuities (kink in the field);

• the pressure field is discontinuous across Γ if surface tension exists and the interface has a nonzero curvature, as stated by the Young-Laplace To take into account these discontinuities, the continuous piece-wise polynomial approximations of the standard F EM are enriched with additional discontinuous functions, which bear the required discontinuities. To start with, we have followed previous XF EM works for incompressible fluid flow problems [START_REF] Legay | Strong and weak arbitrary discontinuities in spectral finite elements[END_REF][START_REF] Gross | An extended pressure finite element space for twophase incompressible flows with surface tension[END_REF],

in which velocity and pressure fields along with their weighting functions are approximated as follows:

v( x, t) = 9/27 i=1 N Q9/H27 i ( x) V i (t) + 4/8 j=1 N Q4/H8 j ( x) ψ v ( x) W v j (t), (7) p
( x, t) = 4/8 i=1 N Q4/H8 i ( x) P i (t) standard approximation + 4/8 j=1 N Q4/H8 j ( x) ψ p ( x) W p j (t) enrichment (8)
where ψ v ( x) and ψ p ( x) are the enrichment functions for velocity and pressure fields. These enrichments satisfy the local partition of unity [START_REF] Melenk | The partition of unity finite element method: Basic theory and applications[END_REF][START_REF] Babuska | The partition of unity method[END_REF] of Q 4 elements in two dimensions and H 8 elements in three dimensions. They are associated with additional degrees of freedom W v j , W p j for the velocity and pressure fields, respectively. As discontinuities depend on the interface location inside the spatial domain, implicitly defined by Φ( x, t) = 0, the enrichment functions are obviously defined using the Φ function:

-for the velocity field enrichment, one selects the |Φ( x, t)| function, which is continuous but with discontinuous gradient across the interface; -for the pressure field enrichment, one selects the Sign(Φ( x, t)) function, which is discontinuous across the interface.

Moreover the enrichment functions must vanish at the element boundary nodes to ensure that no additional contribution is introduced to adjacent elements which share these nodes. For the velocity field this condition is fulfilled by the shifted enrichment [START_REF] Moes | A computational approach to handle complex microstructure geometries[END_REF]:

ψ v ( x, t) = |Φ( x, t)| - 4/8 i=1 N Q4/H8 i ( x) |Φ i (t)|, (9)
where Φ i (t) are nodal values of the Φ function at time t. This enrichment function is depicted in fig. 4(a) for a rectilinear interface. For the pressure field

ψ p ( x, t) = SignΦ( x, t) - 4/8 i=1 N Q4/H8 i ( x) SignΦ i (t). (10)
The additional degrees of freedom introduced by the enrichment can be easily taken into account with an in-house finite element code. For elements crossed by the interface the global approach is the same as for elements not crossed, Finally, the standard finite element approximation ( 6) is also used for the Level Set function in these two-fluid elements as it is continuous over the whole computational domain Ω.

Resulting algebraic systems and solution algorithm

Appropriately inserting standard F EM approximations ( 6) and XF EM ones eq. ( 7)-( 10) and into the three weak integral forms (5) results in the three following algebraic systems:

([M v ] + ∆t[K v (V * )]) [V * ] = [M v ][V t ] + ∆t[F v ], (11)
[K p ] [P ] = [F p ], (12)
([M Φ ] + ∆t[K Φ ]) [Φ t+∆t ] = [M Φ ][Φ t ]. (13)
where

[V * ] and [V t ]

Specific computations in two-fluid elements

The velocity and pressure weak integral forms involve the calculation of various contributions that are discretized in space with XF EM approximations [START_REF] Osher | Front propagating with curvature dependent speed: Algorithms based on hamilton-jacobi formulations[END_REF] in mesh elements crossed by the interface (two-fluid elements). Therefore, specific computations have to be implemented to accurately perform numerical integrations in these two-fluid elements as follows: i) Identify every two-fluid elements, denoted Ω Γ e ; ii) Select appropriate enriched field approximations; iii)

Split every of two-fluid elements in two geometrical sub-domains Ω Γ e1 and Ω Γ e2 conforming to the interface Γ, inheriting this way piece-wise homogeneous material properties; iiii) Perform numerical integration of domain and boundary integrals over the two sub-domains (Ω Γ e1 and Ω Γ e2 ).

Partition of two-fluid elements

The partition of two-fluid elements enables to perform accurate numerical integrations on homogeneous piece-wise continuous sub-domains. This partition is strongly related to the interface geometry, so to proceed efficiently separate implementations have been specifically designed for the related finite elements (Q 9 in 2D and H 27 in 3D), using ad-hoc strategies. To split two-fluid elements, one first look for intersecting points between the interface and the element sides.

Due to piece-wise quadratic approximation of the Level Set function [START_REF] Hirt | Volume of fluids (vof) method for the dynamic of free boundaries[END_REF] within each element, the searched intersecting points solve for the following second order equation in the parametric element:

Φ( ξ) = 9/27 i=1 N Q9/H27 i ( ξ) Φ i = 0, with ξ j = ±1
whose coefficients depend on the nodal values Φ i (t). The sub-domain partition inside a two-fluid element is built using these intersecting points.

Two-dimensional implementation

In two dimensional cases the two-fluid elements are partitioned with biquadratic triangles T 6 in order to match the two material sub-domains Ω Γ e1 and 265 Ω Γ e2 . A recursive algorithm has been designed to locate 2 n + 1 (n ∈ N) "equidistant" points along the interface portion into this two-fluid element. So, the partition is performed with the number of points along the interface satisfying a user-provided accuracy, as shown in fig. 5.

Three-dimensional implementation

In three-dimensional cases the geometry of the interface can be complex, as it is quadratic in each space directions. Many intersecting cases have to be managed, thus, we have implemented several partitioning techniques devoted to the most generic cases. The first one is based on an a priori partitioning of the H 27 hexahedra reference element into a pre-selected patch of tetrahedra not conforming to the interface. We have implemented either a 5 tetrahedra patch, cf. fig. 6(a), which leads to an asymmetric partition, or the symmetrical 24 tetrahedra patch, cf. fig. 6(b). Obviously, a subset of these first level tetrahedra are crossed by the interface requiring to split them into second level tetrahedra to end-up with a tetrahedra partitioning conforming to the interface, as depicted in fig. 7. The outcome of such first level partitioning of an hexahedron into preselected tetrahedra patches is to achieve a computationally efficient procedure. However, the asymmetric 5 tetrahedra patch can generate some relative inaccuracy in numerical integration compared to the symmetric one. The reason is that the finite element basis of tetrahedron is incomplete, so numerical integration errors spatially vary in the frame of reference of the master element. Therefore, it turned out to be more efficient and accurate to split two-fluid elements with a second approach, which uses as many as possible H 27 hexahedra subdomains (cf fig. 8) and otherwise tetrahedra T E 10 . This mixed splitting technique presents the best compromise between geometrical flexibility 290 and computational accuracy and efficiency.

Computations of domain integrals

For single-fluid elements (not crossed by the interface, Ω e ), the element contribution I e to the global weak integral form is computed by numerical integration on a reference parametric element Ω 0 , using Gauss quadratures (with npg integration points, of ξ pg coordinates in the reference element frame and ω pg weighting coefficients):

I Ωe = npg pg=1 ω pg f ( ξ pg ) |J( ξ pg )|.
where f ( x) is a continuous function that depends on standard finite element basis functions and/or their derivatives, J is the jacobian of the geometric transformation from the reference element to the physical one.

295

For two-fluid elements (crossed by the interface, Ω Γ e ), the numerical integration is undertaken differently from the previous case as integrands are discon- tinuous. So, according to procedures described in section 4.1 these elements are split into geometrical sub-domains T e (triangles T 6 in two dimensions, tetrahedra T E 10 or hexahedra H 27 in three dimensions), where all integrands are piece-wise continuous, cf. fig. 9. So, the element integral results in the sum on all sub-domains of quadratures performed on npg T integration points of the reference element T 0 associated with every split subdomain:

I Ω Γ e = Te npg T pg=1 ω T pg f ξ T pg |J( ξ T pg )| |J T ( a T pg )|,
where J T is the jacobian of the transformation from T 0 to T e in Ω 0 , ξ T pg are the coordinates of the integration points a T pg of the reference element T 0 in the space of Ω 0 and ω T pg the corresponding weights. The partitioning method into homogenous sub-domains is only used for integration purpose, so it does not introduce any additional degrees of freedom. However, the price to pay is that they requires a higher number of integration points than single-fluid elements. Indeed, the integral over two-fluid elements use XF EM enrichment functions along with a second coordinate transformation jacobian J T that both increase the polynomial degree of integrands.

T 0 a b J T ¨¨¨¨¨¨¨B Ω 02 Ω 01 Ω 0 Γ ξ η J $ $ $ $ $ $ $ X

Computations of boundary integral along the interface

The boundary integrals that appear into the weak integral form (5) come generically from two boundary subsets: i) the external domain boundary ∂Ω (where Neumann boundary conditions are prescribed); ii) the interface Γ between the two material sub-domains Ω 1 and Ω 2 , where surface tension induces a momentum source term according to [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluids with a free surface[END_REF]. The related integral reads:

I Γ = Γ τ χ( x Γ ) n( x Γ ) + ∂τ ∂s t g ( x Γ ) • δv( x Γ ) dΣ. ( 14
)
For an interface given as an implicit curve Φ( x) = 0, the unit normal vector to the interface (arbitrarily oriented from Ω 1 to Ω 2 ) is given by the normalized

gradient of Φ: n( x ∈ Γ) = ∇Φ || ∇Φ||
. A direct integration of the first term appearing in the boundary integral [START_REF] Chang | A level set formulation of eulerian interface capturing methods for incompressible fluid flows[END_REF] requires the computation of the interface curvature χ, which involves second order derivatives of Φ with respect to global coordinates. So, it is more efficient and accurate to use alternate approaches that only relies on first order derivatives. In two dimensions the Frenet relationship expresses the curvature as the derivative of the tangent vector t g (s) along the interface:

d t g ds (s) = χ(s) n(s), (15)
Integrating by parts the related terms ends up with only first order derivatives:

I 2D Γ = - Γ τ t g • ∂ δv ∂s ds, ( 16
)
where s is the curvilinear abscissa along the interface. In three dimensions, despite it no longer exists any relation in the form of eq. ( 15), integration by parts can however be performed thanks to the Laplace-Beltrami differential [START_REF] Chessa | An enriched finite element method and level sets for axisymmetric two-phase flow with surface tension[END_REF][START_REF] Demlov | An adaptive finite element method for the Laplace-Beltrami operator on implicitly defined surfaces[END_REF] that similarly makes only first order derivatives to appear. Using the previous definition of the unit normal along the interface, it reads [START_REF] Gross | Finite element discretization error analysis of a surface tension force in two-phase incompressible flows[END_REF][START_REF] Gross | An extended pressure finite element space for twophase incompressible flows with surface tension[END_REF][START_REF] Fries | The intrinsic XFEM for two-fluid flows[END_REF]:

I 3D Γ = - 3 i=1 Γ τ Ī -n Γ ⊗ n T Γ e i • ∇ Γ δv i dΣ, ( 17
)
where e i is the i th unit vector of the cartesian frame, δv i is the i th component of the test function associated with the velocity vector.

Numerical results

The developed numerical model has been first validated on several physical configurations that have analytic solutions, enabling us to evaluate the accuracy of our numerical results. The first of such considered configuration is the gravitational two-fluid flow in a vertical channel, in which the interface is vertical so only one velocity component does not reduce to zero. In this case, any contrast of fluid densities and/or dynamic viscosities induces a discontinuity of velocity gradient across the interface. The second considered configuration is a capillary two-fluid flow, in which surface tension induces a pressure discontinuity across the interface.

Discontinuity of velocity gradient across an interface

A simple configuration that brings into play a discontinuity of velocity gradient across the interface is a steady gravitational two-fluid flow in a vertical channel, in which the interface is vertical. Indeed, any contrast in material properties of the two fluids located on each side of the interface induces a discontinuity in the normal derivative to the interface of the velocity field (discontinuous velocity gradient). Plane and cylindrical shapes of the vertical interface have been considered in this configuration to assess the capabilities of the developed numerical model.

Plane interface

For parallelepiped vertical channels in which the interface is parallel to one pair of lateral walls, the steady-state analytic solution can be obtained by integrating the related momentum equation ( 2), as it simply reduces to a second order partial differential equation of the vertical velocity component in each fluid:

µ α ∂ 2 w α ∂x 2 + ∂ 2 w α ∂y 2 = ρ α g,
with no-slip boundary condition along the vertical channel walls and velocity continuity at the interface (w 1 (x I ) = w 2 (x I )), the interface being located at x I in the x direction. As the interface is plane (zero curvature) the stress vector is continuous across it, but if the two fluids have different material properties a discontinuity of the normal derivative of the velocity field takes place across the interface. Further simplifying the problem assuming an infinite channel width in the y direction, the 1D analytic solution reads:

w 1 (x) = ρ 1 µ 1 g x 2 2 + C 11 x for x ∈ [0, x I ] ; w 2 (x) = ρ 2 µ 2 g (x 2 -L 2 ) 2 + C 21 (x -L) for x ∈ [x I , L] ; (18)
C 21 =   µ 2 (2ρ 2 -ρ 1 ) + µ 1 ρ 2 (2 L 2 x 2 I -1) 2µ 2 [µ 1 (x I -L) -µ 2 x I ]   gx 2 I ; C 11 = C 21 (1 - L x I ) + ρ 2 µ 2 (1 - L 2 x 2 I ) - ρ 1 µ 1 g x I 2 .
This analytic solution eq. ( 18) is piece-wise quadratic, so it could be exactly represented by piece-wise Q 2 finite element approximations provided the mesh coincides with the interface, i.e., it exists only single-fluid elements. When this latter condition is no longer satisfied, i.e., it exists two-fluid elements, only approximate solutions are obtained with piece-wise Q 2 finite element approxi-mations. Therefore, this basic gravitational two-fluid flow is one of the simplest configuration to evaluate the accuracy of our numerical model.

To start with, the two salient ingredients of our numerical model (partition of two-fluid elements and appropriate enrichment of approximate fields) have been evaluated on a cubical domain [0, 1] 3 , in which a vertical plane interface is located at x I = 0.5. Fluid properties are set to

ρ 1 = ρ 2 = 1/g, µ 1 = 1, µ 2 = 50
for densities and dynamic viscosities, respectively. At first, the computational domain is discretized with the minimal uniform mesh not conforming to the interface, made-up of only three the F EM approximations completely fail to reproduce the analytic solution in most of the computational domain. Indeed, piece-wise continuous polynomial approximation has no capability to resolve any discontinuity in velocity gradient within an element. Therefore, the stress balance across the interface located inside any two-fluid element cannot be satisfied. As a consequence, the velocity profile computed on such a coarse mesh significantly departs from the analytic solution, cf. fig. 10(a). On the other hand, XF EM approximations accurately reproduce the discontinuity of velocity gradient across the interface, see figure 10(b). So, it is able to recover the parabolic analytic solution in the two side single-fluid elements and it produces a very good approximation in the two-fluid element. In the latter, the analytic solution is piece-wise quadratic, whereas XF EM approximations are piece-wise cubic on each side of the interface (it results from the product of a linear part eq. ( 7) and a quadratic one eq. ( 9)). Refining uniformly the mesh in the x direction while keeping one twofluid element symmetrically crossed by the interface improves the accuracy of both computations. The related root mean square error (RM SE) computed on all mesh nodes for various mesh resolutions are reported in table 1 for both F EM approximation (with split two-fluid element) and XF EM one. For F EM computations, mesh refinement slowly improves accuracy, but the solution quality is desperately poor even for the finest mesh considered, for which the RM SE is only O(10 -4 ). On the other hand, for XF EM approximation the RM SE is impressively very good O(10 -9 ) from the coarsest minimal mesh (3 × 1 × 1)

H
and slightly drops down to numerical accuracy O(10 -10 ) for the finest mesh considered. It turns out that mesh refinement provides an improvement mainly within the two-fluid element, since a very accurate solution is already achieved elsewhere in single-fluid elements even from the coarsest mesh.

Once the accuracy of the developed XF EM model has been assessed, let us now investigate its capabilities to deal with arbitrarily high contrasts in material properties even on very coarse meshes. Various ratios of dynamic viscosities have been undertaken and it turns out that the accuracy of the computed XF EM solution only decreases very slightly when increasing this ratio. For a one thousand ratio in dynamic viscosities the RM SE is only four times larger than that for a fifty ratio. The x-profiles of the vertical velocity component are plotted in figure 11(a) for µ 1 = 1, µ 2 = 1000, ρ 1 = ρ 2 = 1/g and in figure 11(b) for Table 1: Gravitational two-fluid flow in a vertical parallelepiped channel with vertical plane interface located at x I = 0.5, ρ 1 = ρ 2 = 1/g, µ 1 = 1, µ 2 = 50. Root mean square error

H 27 mesh (N x × 1 × 1) h x = 1/N x RM SE F EM RM SE XF EM
(RM SE = 1 n n i=1 (w i -w ana i ) 2
, n stands for total mesh nodes) for F EM (with split two-fluid element) and XF EM approximations for various not conforming uniform mesh resolutions. interface locations, ratios in material properties and mesh resolutions. Among them, representative x-profiles of the vertical velocity component are depicted in figure 12 for an interface located at x I = 0.31, a ratio of dynamic viscosities 415 of 100, two ratios of densities, but roughly ten times larger densities than in previous cases. Here again, the RM SE of our XF EM solutions are O(10 -9 )

µ 1 = 1, µ 2 = 500, ρ 1 = 1/g, ρ 2 =
for a 50 × 1 × 1 H 27 not conforming uniform mesh.

Cylindrical interface

To further illustrate our model capabilities let us now consider gravitational this problem has no analytic solution two independent approximate solutions have also been considered to assess the quality of our XF EM computations: i) the approximate solution computed by standard F EM on a refined mesh (1827 H 27 elements) conforming to the interface (only single-fluid elements); ii) the analytic solution of an axisymmetric problem in which the interface coincides with the actual one (at r I = 3.8) and whose solid wall is the inscribed cylinder of the parallelepiped (r w = 6). Among these two latter approximate solutions, the former is definitely the most reliable since it accounts for the actual geometry, whereas the latter considers the outer flow into a pipe instead of a parallelepiped channel. Moreover, the mesh of the standard F EM fits the interface, so the stress balance across it is naturally satisfied by adjacent elements in this case.

The vertical velocity profiles of all these approximate solutions are plotted versus a horizontal bisecting line in figure 14. A close inspection of this plot reveals

(b) Close-up in the vicinity of r = 0.

3.2 that our XF EM solution computed on the 6 × 6 × 3 H 27 uniform mesh (red dotted line) underestimates the reference solution (blue crosses), as the Q 2 approximation of the Level Set function is not accurate enough on this mesh.

Indeed, simply refining this mesh by a two factor in each horizontal directions (12 × 12 × 3 H 27 , red bold dots) significantly improves the solution quality, as the shape of the cylindrical interface is now much better resolved than on the coarsest mesh. On the other hand, the analytical solution of a related axisymmetric problem (green solid line) slightly overestimates the reference solution, especially for low radius of higher velocity. Close-up views in the vicinity of the profile extrema (r = 0, see figure 14(b)) and two-fluid interface (r = r I , see figure 14(c)) reveal that our XF EM solution computed on a 12 × 12 × 3 H 27 provides a very satisfactory result compared to the reference one computed with standard F EM on a 1827 H 27 mesh conforming to the interface.

To conclude this section dealing with discontinuities in velocity gradient across the interface, it turns out that all computations performed on these gravitational two-fluid flows in a vertical parallelepiped channel demonstrate that our XF EM implementation enables to get the right discontinuity of the normal derivative of velocity inside two-fluid elements, even for arbitrary high contrasts in fluid properties. Our implementation that properly deals with partition of two-fluid elements to perform accurate numerical integration along with appropriate enrichment of the velocity field produces very accurate results on astonishing coarse meshes. In proceeding that way the mesh size of two-fluid elements is no longer related to any contrast in fluid properties unlike in most averaging one-fluid models. So, accurate results can be achieved on relative coarse meshes, provided they faithfully approximate the interface shape in two-fluid elements by means of Q 2 finite element approximation of the Level Set function.

Discontinuity of pressure field across an interface

The physical problems we are interested in led us to account at the macroscopic scale for surface tension that acts at two-fluid interfaces. This section is therefore devoted to assess the capabilities of the present eulerian XF EM model to deal with such pressure discontinuity configurations. Let us consider a simple test case drawn from [START_REF] Gross | An extended pressure finite element space for twophase incompressible flows with surface tension[END_REF] in which a surface force density is assumed to act at a planar interface, disjoining this way the pressure jump across the interface from any tricky curvature computations. This configuration is realized in a closed container filled with two liquids at rest, separated by a planar vertical interface where a uniform and constant force density acts in the direc-tion normal to the interface, i.e. [ n σ n] Γ = T n . The latter produces a pressure jump across the interface whose magnitude is the applied force density itself [START_REF] Legay | Strong and weak arbitrary discontinuities in spectral finite elements[END_REF] (eq. ( 8) and [START_REF] Chessa | An extended finite element method for two-phase fluids[END_REF], denoted XF EM 1 ) in figure 15(b), the proposed pressure enrichment (eq. ( 8) and [START_REF] Gross | An extended pressure finite element space for twophase incompressible flows with surface tension[END_REF], denoted XF EM 2 ) in figure 15(c) and finally the theoretical step function. It turns out that F EM and XF EM 1 approximations produce very similar oscillatory solutions, instead of the theoretical step function. This is the expected behavior for F EM , in which the pressure approximation is piecewise linear inside elements and continuous from one to another, hence it has no way to account for any pressure jump. On the other hand, XF EM 1 has been designed to manage discontinuous pressure fields [START_REF] Legay | Strong and weak arbitrary discontinuities in spectral finite elements[END_REF][START_REF] Gross | An extended pressure finite element space for twophase incompressible flows with surface tension[END_REF]. Surprisingly, in this basic test case it fails to reproduce the expected pressure jump across the interface and it does not improve noticeably the F EM solution quality. The XF EM 1 pressure profile depicted in figure 15(b) presents a too tiny jump across the interface aside a large linear part that has the same slope as the F EM solution. Consequently, both pressure approximations produce a pressure gradient inside the two-fluid element, which induces in turn a purely spurious fluid flow at steady state.

([p] Γ = T n ),
Therefore, this basic test case points out an unexpected failure of the XF EM 1 enrichment to model pressure jumps, though it is widespread in the XF EM community of incompressible two-fluid flows. Thus, this disappointing result led us to evaluate alternate pressure enrichment functions. For this problem where the interface location only depends on one space dimension, the proposed enrichment function of eq. ( 19) outperforms the other tested ones. It reads:

ψ p ( ξ, t) = 1 -ξ 2 SignΦ( ξ, t) (19)
where ξ stands for the one dimensional coordinate along the normal direction to the interface in the parametric two-fluid element. The proposed XF EM 2 pressure approximation defined by eq. ( 8) and [START_REF] Gross | An extended pressure finite element space for twophase incompressible flows with surface tension[END_REF], depicted in figure 15(c), faithfully reproduces the pressure jump across the interface and furthermore computed pressures have an excellent accuracy at mesh nodes. However, owing to its piece-wise quadratic nature it nevertheless introduces numerical undershoot and overshoot within the two-fluid element on either side of the pressure jump. Interestingly, the enrichment locality prevents them to propagate outside the two-fluid element, unlike the oscillatory solution arising from both F EM and XF EM 1 approximations.

These three approximations are convergent in spatial resolution, so refining the mesh size improves their accuracy. For F EM and XF EM 1 the pressure gra-curacy, as from the coarsest mesh (O(10 -13 ) and O(10 -20 ), for pressure and velocity errors, respectively). This highly accurate pressure field solution computed at mesh nodes is the one expected from a nicely working XF EM approximation on such an academic problem. But the most outstanding feature of this pressure enrichment is that the resulting pressure field does not induce any spurious velocity field unlike the two other pressure approximations. Indeed, the computed pressure overshoot and undershoot arising on either side of the interface within the two-fluid element are perfectly antisymmetric with respect to the interface thanks to the chosen enrichment function of eq. ( 19). Therefore, they do not induce any artificial pressure gradient inside the two-fluid element, and results in a nicely welcome spurious-free velocity field. 8) and ( 10)) and proposed (eq. ( 8) and ( 19))

H 27 mesh hx = 1/Nx F EM XF EM 1 XF EM 2 Nx × 1 × 1 hy = hz =
XF EM approximations for various not conforming uniform mesh resolutions.

Several other test cases have been undertaken and it turns out that the proposed pressure approximation XF EM 2 provides a significant improvement with respect to XF EM 1 as it successfully catches the pressure discontinuity of any step function, regardless of the jump magnitude and mesh size of the two-fluid element. However, it is noteworthy that the introduced pressure enrichment function of eq. ( 19) was tailored for this one-dimensional interface configuration, but its generalization to any complex three-dimensional cases is not straightforward.

Discussion

Various other three-dimensional two-fluid flow configurations involving both weak discontinuities (velocity gradient discontinuity) and strong ones (pressure field discontinuity) have been computed with the present XF EM implementation. From this experience it turns out that the two salient ingredients to get accurate computations are the following: i) select appropriate enrichment functions [START_REF] Moes | Discontinous enrichment in finite elements with a partition of unity method[END_REF][START_REF] Chessa | An extended finite element method for two-phase fluids[END_REF][START_REF] Legay | Strong and weak arbitrary discontinuities in spectral finite elements[END_REF][START_REF] Fries | The extended/generalized finite element method: An overview of the method and its applications[END_REF] to model the desired discontinuity type from physical problem ; ii) perform appropriate numerical integrations for the specified problem data (physical, numerical and computational ones). Provided these two necessary conditions are simultaneously satisfied, then XF EM approximations enable to compute outstanding quality solutions regardless of discontinuity magnitude and mesh size of the two-fluid element, as far as the interface is accurately represented by Q 2 approximation of the Level Set function.

Appropriate enrichment functions

For velocity gradient discontinuity within two-fluid elements (weak discontinuity) the XF EM approximations drawn from the literature [START_REF] Moes | A computational approach to handle complex microstructure geometries[END_REF][START_REF] Legay | Strong and weak arbitrary discontinuities in spectral finite elements[END_REF], defined by eq. ( 7) and ( 9) and implemented in our 3D unconditionally stable projection algorithm framework [START_REF] Rabier | Computation of free surface flows with a projection fem in a moving mesh framework[END_REF] performed impressively well in representative test cases considered, such as in section 5.1. On the other hand, as mentioned in section 5.2

we didn't succeed to get the expected pressure field discontinuity (strong discontinuity) with approximations drawn from [START_REF] Legay | Strong and weak arbitrary discontinuities in spectral finite elements[END_REF][START_REF] Gross | An extended pressure finite element space for twophase incompressible flows with surface tension[END_REF], defined by eq. ( 8) and [START_REF] Chessa | An extended finite element method for two-phase fluids[END_REF]. It turns out that in the framework of our projection algorithm XF EM 1 approximations produce singular stiffness matrix of the Pressure Poisson Equation, eq. ( 12). Indeed, after close inspection thanks to extensive symbolic calculations, the determinant of the assembled stiffness matrix is analytically zero for the example considered in section 5.2, whatever the not conforming uniform mesh resolution and interface location within the two-fluid element. The related numerically discretized stiffness matrix has in turn a very close to singular determinant, whenever it is not exactly zero owing to round-off errors, resulting therefore in an awkward pressure solution. This explains why this XF EM 1 pressure approximations do not provide any improvement with respect of the standard F EM one, even degrading the velocity accuracy as reported in table 2. This was not explicitly mentioned in the paper in which this test case was formerly computed with the same pressure enrichment [START_REF] Gross | An extended pressure finite element space for twophase incompressible flows with surface tension[END_REF]. However, similar pressure oscillations localized in the vicinity of the interface were also present despite the authors selected highly refined meshes for this problem where the used approximations are supposed to represent the analytical solution. Fortunately, the proposed quadratic pressure enrichment function of eq. ( 19) resulting in XF EM 2 approximations solved the problem and produced the desired pressure jump across the interface. So, the proposed quadratic enrichment likely suggests that it should be selected in a separate approximation space than the standard one, unlike XF EM 1 . But, it is also noteworthy to recall that in the framework of incompressible fluid flow problems one has to care about the stability of velocity-pressure coupling. Indeed, velocity and pressure approximation spaces are closely linked by the inf -sup (LBB) stability condition. So, an XF EM pressure approximation should also be compatible with the related velocity-pressure coupling algorithm to lead to a reliable implementation.

Appropriate numerical integrations

Once appropriate enrichment functions have been selected, it is of first concern to perform appropriate numerical integrations to compute accurate XF EM solutions. As detailed in section 4 numerical integrations to be performed in two-fluid elements importantly contribute to accuracy. First of all, the splitting algorithm should preserve as far as possible problem symmetries. Secondly, numerical quadratures should be performed according to the polynomial order of integrands, which are now piece-wise continuous over split sub-domains. Indeed, accurate Gauss quadratures require npg points per spatial direction, with npg ≥ 2r -1, r being the polynomial order of the integrand. Finally, particular attention should be paid to ill-conditioning of algebraic systems resulting from XF EM approximations. Indeed, the condition number they produces can exceed by up to three orders of magnitude standard F EM one, as from first order hexahedral finite elements [START_REF] Sauerland | The stable XFEM for two-phase flows[END_REF], meanwhile this becomes even worse for present quadratic approximations. Not only ill-conditionning induces very low convergence rates with iterative solvers, but it also leads in extreme cases to inaccurate computations with direct solvers. Indeed, when it exists a huge scaling difference between algebraic system coefficients related to standard degrees of freedom compared to enrichment ones, diagonal scaling slightly improves the conditioning, but does not prevent from inaccurate XF EM solutions. Therefore, we have introduced an additional scaling coefficient associated with the enrichment functions for which a well chosen value ensures both good convergence and accurate discontinuity solution across the interface. The influence of this scaling coefficient on the solution accuracy is illustrated in figure 17 on the gravitational two-fluid flow example of section 5.1.1. In order to clearly bring to the fore its effect, computations are performed on the same minimal not conforming to the interface uniform mesh with the M U M P S direct solver [START_REF] Amestoy | MUMPS Web page[END_REF] interfaced with the P ET Sc library [22]. For the considered material property contrast (ρ 1 = ρ 2 = 1/g, µ 1 = 1, µ 2 = 50) and interface location (x I = 5 12 ) a roughly optimal value to get the best efficiency (requested accuracy and lowest computational cost) is S coef = 100, see figure 17(a). Departing away from this optimal scaling coefficient not only degrades the computation accuracy, as depicted in figure 17(b) for various scaling coefficients ranging on four orders of magnitude, but can ultimately lead to numerical divergence of the projection algorithm. Finally, the interface location within the two-fluid element also influences the optimal value of the scaling coefficient. Indeed, computing with the same scaling coefficients as previously but for an interface located at x I = 7 12 produces here again unacceptable solutions, see figure 17(c). Therefore, based on these results the scaling coefficient used in our computations has been adaptively related to both physical and geometrical considerations. On the one hand it accounts for the physical problem discontinuity (material property contrast for the discontinuity of velocity gradient, surface tension for pressure discontinuity, etc.), and on the other hand it also accounts for the interface configuration within two-fluid elements throughout the resulting volume ratio. It reads as follows:

S coef ∝ M in( P p 1 P p 2 , P p 2 P p 1 ) × M ax(

Ω Γ e1 Ω Γ e2 , Ω Γ e2 Ω Γ e1 ) (20)
where P p α is the material property of fluid α involved in the physical problem discontinuity. This scaling coefficient enabled us to get close to optimal values in the considered test cases. However, further analysis is obviously required to extend it to more complex geometrical configurations.

Conclusion

The developed XF EM implementation aims at computing incompressible

Figure 1 :

 1 Figure 1: Basic kinematic approaches for two-fluid flow problems.

Figure 2 :

 2 Figure 2: Interface defined as the zero-level of a level set function.

equation ( 3 )

 3 , cf. figure3(b). This kind of discontinuity is referred to as strong discontinuity. Continuous velocity field with discontinuous normal gradient across Γ.

  Discontinuous pressure field across Γ.

Figure 3 :

 3 Figure 3: Considered discontinuities of velocity and pressure fields across the interface Γ.

  (a) Enrichment with normal gradient discontinuity across the interface Γ. (b) Enrichment with jump across the interface Γ.

Figure 4 :

 4 Figure 4: Shifted enrichment functions that vanish at element boundary nodes.

  provided one introduces a family of basis functions made-up of the classical basis functions supplemented with the enrichment ones. In the present case, the standard basis function family {N Q9/H27 i } i=1,n is supplemented with the enrichment functions {N Q4/H8 i ψ} i=1,m , in the following way: { Ñi } i=1,n+m = {N where n = 9 for Q 9 and n = 27 for H 27 , meanwhile m = 4 for Q 4 and m = 8 for H 8 approximations, respectively.

  are global vectors of velocity degrees of freedom at intermediate and previous time-steps, [M v ] and [K v ] are global mass and stiffness matrices of the velocity algebraic system, respectively, [F v ] is the load vector of the discretized momentum equation. For the projection step, [P ] is the global vector of pressure degrees of freedom, [K p ] is the global stiffness matrix of the pressure algebraic system, [F p ] is the load vector of the discretized pressure poisson equation. Finally, for the algebraic system related to the Level Set function, [Φ t ] and [Φ t+∆t ] are global vectors of Level Set degrees of freedom at current and next time-steps, [M Φ ] and [K Φ ] are global mass and stiffness matrices of the Level Set algebraic system. This Finite Element model has been implemented in the P ET Sc environment [22], which give us access to High Performance Computations on parallel computers. The algebraic system associated with the momentum advectiondiffusion step is nonlinear, so it is iteratively solved at each time step of the solution algorithm by a Newton-Raphson algorithm using the SN ES environment. On the other hand, the algebraic systems associated with the pressure Poisson projection step and Level Set function are linear, so they are solved at each time step of the solution algorithm with the KSP routines. All the linearized systems are solved with the Bi-Conjugate Gradient Stabilized iterative solver, preconditioned with the Additive Schwartz Method, efficiently implemented in the P ET Sc library.

  Splitting with 5 points along Γ.

  Splitting with 9 points along Γ.

Figure 5 :

 5 Figure 5: Two-fluid element partitioning with different number of bi-quadratic triangles in 2D.

  (a) Splitting hexahedron into 5 tetrahedra. (b) Splitting hexahedron into 24 tetrahedra.

Figure 6 :

 6 Figure 6: First level partitioning of reference hexahedron into preselected tetrahedra patches.

  (a) Final splitting of the 5 tetrahedra patch. (b) Final splitting of the 24 tetrahedra patch.

Figure 7 :

 7 Figure 7: Second level partitioning of a hexahedron into conforming to the interface tetrahedra.

  (a) Hexahedron split into 2 hexahedra. (b) Hexahedron split into mixed tetrahedra and hexahedra.

Figure 8 :

 8 Figure 8: Partitioning of a hexahedron into conforming to the interface hexahedra and tetrahedra.

Figure 9 :

 9 Figure 9: Domain integration procedure in 2D two-fluid element split into T 6 triangles (quadrature points are only displayed in the upper-left triangle).

  27 finite elements in the x direction (two singlefluid elements apart from one two-fluid element) and only one element in the invariant y and z directions. The x-profiles of vertical velocity component are plotted for both standard F EM approximations (with split two-fluid element) in figure 10(a) and XF EM ones in figure 10(b), along with the analytic solution of eq. (18). Solely supplemented with a splitting of the two-fluid element F EM solution (split two-fluid element).

  XF EM solution.

Figure 10 :

 10 Figure 10: Gravitational two-fluid flow in a vertical parallelepiped channel with vertical plane interface located at x I = 0.5, ρ 1 = ρ 2 = 1/g, µ 1 = 1, µ 2 = 50. Vertical velocity plot: analytical solution (green solid line), computed solutions on 3 × 1 × 1 H 27 not conforming uniform mesh (red dotted line).

100 ρ 1 ,µ 1 = 1 ,µ 1 = 1 ,

 11111 along with the corresponding analytic solutions. Astonishingly, the RM SE remains O(10 -9 ) whatever being the 410 ratio in material properties and its related discontinuity in velocity gradient, despite our XF EM computations are performed on the minimal 3 × 1 × 1 H 27 not conforming uniform mesh. Many other tests were carried out for various µ 2 = 1000, ρ 1 = ρ 2 = 1/g, RM SE = 7.82877 10 -9 . µ 2 = 500, ρ 1 = 1/g, ρ 2 = 100 ρ 1 , RM SE = 7.92656 10 -9 .

Figure 11 :

 11 Figure 11: Gravitational two-fluid flow in a vertical parallelepiped channel with vertical plane interface located at x I = 0.5 for high contrasts in fluid properties. Vertical velocity plot: analytical solution (green solid line), XF EM solutions computed on a 3 × 1 × 1 H 27 not conforming uniform mesh (red dotted line).

  µ 1 = 1, µ 2 = 100, ρ 1 = 1, ρ 2 =1000, RM SE = 2.52914 10 -9 . Close-up of velocity profile (a) in the interface vicinity.

µ 1 = 1 ,

 11 µ 2 = 100, ρ 1 = 1, ρ 2 = 100, RM SE = 4.76565 10 -9 . Close-up of velocity profile (c) in the interface vicinity.

Figure 12 :

 12 Figure 12: Gravitational two-fluid flow in a vertical parallelepiped channel with vertical plane interface located at x I = 0.31 for high contrasts in fluid properties. Vertical velocity plot: analytical solution (green solid line), XF EM solutions computed on a 50 × 1 × 1 H 27 not conforming uniform mesh (red dotted line).

  two-fluid flows in a vertical parallelepiped channel in which the interface is cylindrical and aligned with lateral walls. For symmetry reasons only one quarter of the computational domain is modeled and a no-slip condition is prescribed at solid walls. The quarter domain is [0, 6] 3 , in which the cylindrical interface is set at a radius r I = 3.8, fluid properties are set to ρ 1 = ρ 2 = 1/g, µ 1 = 1, µ 2 = 500, for densities and dynamic viscosities, respectively. The computational domain is at first discretized with a 6 × 6 × 3 H 27 uniform mesh. The iso-values of computed vertical velocity component with our XF EM model are displayed in figure13(a) along with its profile over a horizontal plane in figure13(b).It turns out that the circular interface is poorly represented on such a coarse cartesian mesh by Q 2 approximation of the Level Set function. Therefore, to illustrate the solution sensitivity to mesh resolution we have also considered a twice refined mesh in the horizontal plane (12×12×3 H 27 uniform mesh). Since

  (a) Iso-values of vertical velocity component. (b) Vertical velocity profile over a horizontal plane.

Figure 13 :

 13 Figure 13: Gravitational two-fluid flow in a vertical parallelepiped channel with vertical circular interface located at r I = 3.8, µ 1 = 1, µ 2 = 500, ρ 1 = ρ 2 = 1/g. XF EM solution of the vertical velocity component computed on a 6 × 6 × 3 H 27 cartesian uniform mesh.

  ) Close-up in the vicinity of r I .

Figure 14 :

 14 Figure 14: Gravitational two-fluid flow in a vertical parallelepiped channel with vertical circular interface located at r I = 3.8, µ 1 = 1, µ 2 = 500, ρ 1 = ρ 2 = 1/g. Plot along bisecting line of vertical velocity component: XF EM solutions computed on a 6 × 6 × 3 H 27 uniform mesh (red dotted line) and 12 × 12 × 3 one (bold red dots), standard F EM solution computed on a 1827 H 27 mesh conforming to the interface (blue crosses); analytical solution of axisymmetric problem (solid green line).

1 .

 1 figure15(a), XF EM ones with pressure enrichment drawn from[START_REF] Legay | Strong and weak arbitrary discontinuities in spectral finite elements[END_REF] (eq.[START_REF] Thompson | Use of pseudo-concentrations to follow creeping viscous flows during transient analysis[END_REF]

  XF EM 1 .

  XF EM 2 .

Figure 15 :

 15 Figure 15: Uniform force density acting at a plane interface in a cubic container: interface located at x I = 0.5, Tn = 10 -3 , ρ 1 = ρ 2 = 1, µ 1 = µ 2 = 1. Pressure profiles in the normal direction to the interface: analytical solution (green solid line), F EM and XF EM solutions computed on a 3 × 1 × 1 H 27 not conforming uniform mesh (red dotted lines).

located at x I = 0. 5 ,

 5 Tn = 10 -3 , ρ 1 = ρ 2 = 1, µ 1 = µ 2 = 1. Relative root mean square error (RM SE/F ) for the pressure field, standard RM SE for the velocity field, computed for F EM (with split two-fluid element), standard (eq. (

Figure 17 :

 17 Figure 17: Gravitational two-fluid flow in a vertical parallelepiped channel with vertical plane interface, ρ 1 = ρ 2 = 1/g, µ 1 = 1, µ 2 = 50. Vertical velocity plot: analytical solution (green solid line), computed solutions on 3 × 1 × 1 H 27 non conforming uniform mesh for various scaling coefficients, dotted lines: S coef = 100 (red); 1 (purple); 0.5 (orange); 0.01 (blue).

  two-fluid flows with arbitrary high contrasts in material properties in 3D configurations. This model features a classical shift enrichment function drawn from the literature to account for velocity gradient discontinuities across the interface (weak discontinuity). On the other hand, for the pressure discontinuity across the interface (strong discontinuity) the related enrichment function drawn from the literature produces an analytically zero determinant in the stiffness matrix associated with the pressure Poisson equation involved in our projection algorithm. Therefore, we have been led to introduce a new quadratic enrichment function which satisfactorily overcomes the problem. Moreover, to achieve accurate numerical integrations for arbitrary high contrasts in material properties, a splitting of two-fluid elements is performed prior to any numerical integrations.Proceeding this way enables us to perform accurate numerical integrations with Gauss quadratures on piece-wise continuous integrand over each subdomain.Finally, an ad -hoc scaling coefficient that accounts for both physical and geometrical considerations has been introduced to alleviate ill-conditioning inherent to classical XF EM approximations.Various validations have been carried out on flow configurations for which analytic solutions are available and for all of them very good solution accuracy has been achieved. It is noteworthy that provided the considered approximations are able to represent the problem solution, this implementation enables to compute outstanding quality solutions regardless of discontinuity magnitude (arbitrary high contrast in material properties) and mesh size of the two-fluid element, as far as the interface is accurately represented by Q 2 approximation of the Level Set function. Moreover, we have shown for the first time in these cases that very accurate solutions are obtained even on coarse meshes, as from the minimal mesh not conforming to the interface. From our point of view, this XF EM feature is definitely a decisive advantage to efficiently deal with 3D incompressible two-fluid flows configurations with arbitrary high contrasts in material properties.

  where V i , P i and Φ i are velocity, Level Set and pressure nodal values at node i,

	respectively, N Q9 i , N Q4 i , N H27 i	and N H8 i	are the basis functions associated with
	node i of the related iso-parametric Lagrange elements Q 9 , Q 4 , H 27 and H 8 ,
	respectively.

Table 2 :

 2 Uniform force density acting at a plane interface in a cubic container: interface

Acknowledgment

This work was partially supported by the French government throughout its funding of the FUI-AAP14 M U SICAS project.

dient becomes steeper and steeper in the two-fluid element, but the theoretical step-shape pressure discontinuity is only poorly approximated until very tiny mesh sizes are considered ( dp dx ∝ T n /h x ) leading to expansive computational costs in 3D. Moreover, pressure oscillations are bounded within up to three elements on either side of the interface. On the other hand, XF EM 2 much accurately approximates the step-shape analytical solution as the overshoot and undershoot become more and more narrow. Indeed, they only spread over one half-width of the two-fluid element, meanwhile they have constant height extrema. These features are illustrated in figure 16, where the pressure profiles computed on a 9 × 1 × 1 H 27 non conforming uniform mesh are plotted for the three related approximations. Quantitative results associated with mesh size refinement are reported in table 2 for the three considered approximations and various uniform mesh resolutions non conforming to the interface. This table gathers the relative root mean square error (RM SE/T n ) of computed pressure fields, along with standard RM SE of velocity fields. As previously observed qualitatively, F EM and XF EM 1 have comparable overall solution qualities, same order of pressure errors (O(10 -2 ) on the finest mesh), despite XF EM 1 produces a three times less accurate velocity field (O(10 -8 ) on the finest mesh) than F EM . On the other hand, XF EM 2 produces an excellent solution ac-