

Version française abrégée 
L'idée clé pour la construction de cette nouvelle méthode ( ) repose sur les décompositions rapides de Helmholtz–Hodge en domaine borné proposées dans [[2]] et analysées dans [[3]]. L'étape de correction vectorielle en vitesse v est rendue purement cinématique. Elle correspond à une méthode de projection à divergence nulle approchée pour obtenir la composante irrotationelle du champ de vitesse prédit . L'étape dynamique de correction du gradient de pression est alors obtenue en introduisant le concept de densité inertielle calculée sur chaque arête du maillage principal et reliant les sommets des inconnues scalaires p , ϕ et ρ . On montre qu'il existe satisfaisant l'équation de consistance locale, , et permettant de recouvrer un gradient de pression entre les sommets. En pratique, comme la condition de compatibilité associée est naturellement satisfaite pour les maillages de type MAC, le calcul de est effectué facilement et de façon rapide à partir de ρ et d'une fonction de partage α , tandis que ϕ est reconstruit à partir de son gradient connu sur les arêtes.
We consider the numerical solution to unsteady incompressible flows with variable mass density and viscosity. This is typically the case of incompressible or lowMach dilatable and/or multimaterial multiphase flows [[13]]. The design of efficient numerical methods, i.e. accurate, fast and robust, is still a challenging issue [[11]], especially for multiphase flows with strong stresses [[4]]: large mass density or viscosity ratios, large surface tension, open boundary conditions, fluid–structure interactions, or multiphysics problems. Many refined techniques have been developed for the fronttracking or frontcapturing of moving and deformable phase interfaces, e.g., [[10], [12], [13], [14]], and most of them use the very popular timestepping projection methods of Chorin–Teman for the flow solver.
Here, we focus on the important feature of velocity–pressure coupling with fully vector splitting methods [[1], [4]], which eliminate the pressure Poisson equation of scalar projection methods and thus the pressure boundary condition. The main objective of this work is to show that the velocity correction step in the vector penalty–projection ( ) method [[1]] can be made fully kinematic and completely independent of the mass density. This is obtained by using fast Helmholtz–Hodge decompositions studied in [[2], [3]]. A similar result is obtained in [[9]] with a scalar projection method. Then, it is shown that the dynamic step of pressure gradient correction can be easily made in a locally consistent way if the spatial discretization satisfies a discrete compatibility condition that is verified for edgebased generalized MACtype unstructured meshes, already introduced in [[1], [7]].


The kinematic vector penalty–projection method for the Navier–Stokes equations 
Let us consider the incompressible Navier–Stokes model problem in the timespace domain with homogeneous Dirichlet boundary condition for velocity v on :
(1){ρ(∂tv+(v⋅∇)v)−div(2μ(ρ)d(v))+∇p=f in (0,T)×Ω,divv=0 in (0,T)×Ω,∂tρ+v⋅∇ρ=0 in (0,T)×Ω,vΓ=0 on (0,T)×Γ,ρ(t=0)=ρ0,v(t=0)=v0 in Ω. Here, Ω is a connected and bounded domain in ( ) with a Lipschitz continuous boundary . The strain rate tensor is denoted by , the mass density by ρ , the dynamic viscosity by and the pressure by p . The force field f usually includes the gravity force and the capillary surface tension supported by the interface Σ between the phases, where is a unit normal vector on Σ, κ its local curvature and σ the surface tension coefficient.
 
 The time semidiscrete scheme ( ) 
The fast ( ) method, proposed here using a firstorder linearly implicit scheme reads with classical notations, for all such that , being the time step, the penalty parameter, and starting from the initial conditions , and in Ω:
(2){(a) Standard prediction step:ρn(v˜n+1−vnδt+(vn⋅∇)v˜n+1)−div(2μ(ρn)d(v˜n+1))+∇pn=fn in Ω,v˜Γn+1=0 on Γ,(b) Divergencefree velocity penalty–projection (VPP): purely kinematic stepεvˆn+1−∇(divvˆn+1)=∇(divv˜n+1) in Ω,vˆΓn+1=0 on Γ,(c) Velocity correction: vn+1=v˜n+1+vˆn+1 in Ω,(d) Find the inertial density ρ‾n such that (see Sect. 2.2):∇(ρ‾nϕn+1)=ρnvˆn+1, with ϕn+1 reconstructed from its gradient vˆn+1:=∇ϕn+1 in Ω,(e) Explicit locally consistent pressure gradient correction: dynamic step∇(pn+1−pn)=−ρnδtvˆn+1=−1δt∇(ρ‾nϕn+1) or (pn+1−pn)=−ρ‾nϕn+1δt in Ω,(f) Advection and fronttracking of density (see Sect. 2.3):ρn+1−ρnδt+vn+1⋅∇ρn=0 in Ω. In the substep (d ), the scalar potential is theoretically equal to up to an additive constant, as a byproduct of the (VPP) substeps . However, to avoid roundoff errors in the numerical computations when ε is very small up to machine precision, typically , since we have , it is shown in [[1], [2]] that it is far better to directly reconstruct locally from its known gradient . This is easily carried out by using the circulation theorem with the Stokes formula along each edge, which gives a very fast numerical procedure. A similar technique is also detailed in Section 2.2 to locally calculate on each edgelinking pressure vertices.
Since with (c ), it is important to notice that the substeps (d ) and (e ) in ((2)) allow us to write the effective velocity and pressure gradient corrections as in [[1], [4]], being the solution by the (VPP) method of substep (b ), which is now purely kinematic and completely independent of ρ :
(3){ρnvn+1−v˜n+1δt+∇(pn+1−pn)=0 in Ω, with ∇(pn+1−pn)=−ρnδtvˆn+1 in Ω. Then, if the calculation of the pressure itself is not necessary for physical reasons, the present method reads more simply:(4)(KVPPε){ρn(v˜n+1−vnδt+(vn⋅∇)v˜n+1)−div(2μ(ρn)d(v˜n+1))+∇pn=fn in Ω,v˜Γn+1=0 on Γ,εvˆn+1−∇(divvˆn+1)=∇(divv˜n+1) in Ω,vˆΓn+1=0 on Γ,Velocity correction: vn+1=v˜n+1+vˆn+1 in Ω,Explicit pressure gradient correction: ∇(pn+1−pn)=−ρnδtvˆn+1 in Ω,ρn+1−ρnδt+vn+1⋅∇ρn=0 in Ω. Another improvement of the present method is that the substep (b ) in ((2)) can be also enforced by the Dirichlet condition and not only with the natural normal component as in [[1], [2], [3]]. Indeed, we can prove by a vanishing viscosity method that it yields a wellposed problem with a unique solution when the boundary Γ is of class or if Ω is a convex domain.
Now by summing the prediction and correction steps using and ((3)), we get that the unique solutions and at each time step satisfy the equations below, even if the practical scheme is never solved in this way but within two steps:
(5){ρn(vn+1−vnδt+(vn⋅∇)v˜n+1)−div(2μ(ρn)d(v˜n+1))+∇pn+1=fn in Ω,vΓn+1=0 on Γ,∇(divvn+1)=εvˆn+1, and ∇(pn+1−pn)=−ρnδtvˆn+1 in Ω,ρn+1−ρnδt+vn+1⋅∇ρn=0 in Ω. This formally shows that the ( ) method does produce the same numerical solutions as the original ( ) one, and thus gives an time accuracy of both velocity and pressure for the firstorder scheme with , as numerically verified in [[1]].
Moreover, the velocity–pressure coupling in ((2)) can be interpreted as a new twostep artificial compressibility scheme for variable density flows satisfying the pressure correction:
(6)pn+1=pn−ρ‾nϕn+1δt in Ω.
 
 Construction of the inertial density field on a staggered generalized MACtype mesh 
The main important issue of the present ( ) method ((2)) is to show the existence at each time step of the socalled inertial density . Indeed, since the irrotational component of is calculated by the (VPP) substep (b ), this also gives easily as a byproduct the scalar potential ϕ such that as detailed further. Then, the inertial density should verify the local consistency condition below, which is required to recover a pressure gradient:
(7)∇(ρ‾ϕ)=ρ∇ϕ=ρvˆ in Ω, with ρ‾=ρ on Γ. A necessary, and also sufficient condition if the domain Ω is simplyconnected, to get such that reads , which gives the local compatibility condition because :(8)rot(ρvˆ)=ρrotvˆ+∇ρ∧vˆ=∇ρ∧vˆ=0, in Ω. This condition is naturally satisfied in the proposed discrete setting below for edgebased generalized MACtype discretisations, since is calculated between two vertices along any edge where the velocity components are located, see Fig. 1. Then from [[2], [3]], since for and defined up to an additive constant, we can choose a continuous potential on the compact domain , which yields .


 Fig. 1. Edgebased generalized MACtype unstructured mesh. Topology of the 3D primal mesh with vertices, edges, faces and an interface Σ: p ,ρ ,ϕ unknowns located at all vertices a or b and velocity components v ⋅t on each edge [a ,b ]. Zoom 
In practice, the compatibility equation in ((7)) can be solved locally on the spatial mesh with a very fast numerical technique. Here, we consider the structured MAC staggered grid [[10]] or the edgebased generalized unstructured staggered mesh of MACtype [[5], [7]] in the context of dispersed multiphase flows. The edgebased generalized MAC mesh shown in Fig. 1 is constructed in the spirit of Whitney's finite elements, see, e.g., [[6]]. A typical configuration at time , in 2D for the sake of clarity, for a triangular primal mesh is shown in Fig. 2. The scalar unknowns p , ϕ , ρ , or are located at the vertices of the primal mesh. The velocity components , as well as the gradients of scalar potentials or are located at the midpoint of any edge linking two neighbour vertices and oriented by a unit vector t . We have verified in [[1]] that the underlying space discretisation yields secondorder errors in for the norm of both velocity and pressure, h being the mesh step. It is also easy to show that the following important properties are exactly satisfied at the discrete scale: and for any scalar ϕ or vector potential ψ . This was numerically verified up to machine precision in [[2]].


 Fig. 2. Detection of an interface represented by a Lagrangian marker chain on the primal mesh. Left: topology of the primal mesh with vertices, edges, cells and the interface – Right: intersections of the connected marker chain cutting across some primal edges. Zoom 
The discrete interface Σ between two phases is represented by a chain of Lagrangian markers having its own connectivity table and defining a polygon with N marker sides (or surfacemarkers based on triangles in 3D). Each marker point belongs to a single cell in the primal mesh and this association is made by a fast algorithm of raytracing [[12]], which also easily determines the intersections of the connected marker chain with the edges separating two neighbour cells. Then, a phase function ξ is defined at each vertex taking the value 1 or 0, depending on whether the point lies inside or outside the polygonal chain.
The calculation of ϕ and is made by reconstruction of a scalar potential by integrating its known gradient along all the edges in the primal mesh. Starting from one point where arbitrarily, we have along any edge , which gives the value when is already known, and so on:
(9)∫abvˆ⋅tdx:=∫ab∇ϕ⋅tdx=ϕb−ϕa, on any edge [a,b]. Now, the calculation of must take into account the density variations and thus the position of the interface Σ. From one side, using the generalized average formula, there exists that is constant along the segment such that:∫abρvˆ⋅tdx=ρ‾∫abvˆ⋅tdx=ρ‾(ϕb−ϕa)=∫ab∇(ρ‾ϕ)⋅tdx. Hence, satisfies the compatibility condition ((7)) along the edge . From another side, denoting the intersection point with the marker chain by and the distance , we have also:∫abρvˆ⋅tdx=∫acρvˆ⋅tdx+∫cbρvˆ⋅tdx=(ρac−a+ρbb−c)vˆ⋅t=(ρac−a+ρbb−c)b−a∫abvˆ⋅tdx=(αρa+(1−α)ρb)(ϕb−ϕa), with α:=c−ab−a. Comparing the two expressions, we get associated with the edge as a weighted average:(10)ρ‾[a,b]=αρa+(1−α)ρb, on any intersected edge [a,b];0≤α≤1.
Moreover, it is explained in the context of Discrete Mechanics [[7], [8]] that the inertial density and the mass density ρ have different physical meanings.
A kinematic scalar projection (KSP) method, similar to the version proposed in [[9]], is now simply obtained by replacing the (VPP) substep (b ) in ((2)) by the standard technique for the divergencefree projection of consisting in solving for ϕ the Neumann–Poisson problem: with on Γ.
 
 A new fronttracking Lagrangianadvection method 
A new simple but accurate Lagrangian fronttracking technique is now introduced for the advection of interface Σ in the substep ((2)(f )) and the updating of density from the velocity field . The algorithm that ensures a good mass conservation of the different phases is briefly described below.
a )  Calculate the barycentric velocity associated with each marker point x from the velocity components on the edges bordering the primal cell where the marker lies. 
b )  Move the markers such that by calculating the new position with the Heun RungeKutta explicit scheme (RK2 or RK4 with the KVPP method of secondorder in time): xn+1=xn+δt2(vbn(xn)+vbn+1(xn+δtvbn(xn))). 
c )  Detect the cells in the primal mesh that are crossed by the updated marker chain with a raytracing technique issued from computer graphics procedures [ [12]] and, according to that, update the phase function ξ at the vertices. 
d )  Calculate the intersection points between the marker chain segments and the edges of the crossed cells in the primal mesh. 
e )  From , calculate the dividing function α on each edge oriented by t and cut across by Σ. 
f )  Update the density ρ , the corresponding viscosity and the inertial mass density , on any intersected edge . 
g )  Compute the local curvature at each marker point x using the osculator circle crossing three consecutive points. 
h )  Compute the force source term modelling the capillary effects on Σ to be included in the force balance on any intersected edge. 
i )  Solve for the flow at time with the method of velocity–pressure coupling. 


Some validations with numerical experiments 
The numerical results are performed with a timestep δt satisfying a CFL value of .
 
 Perfect numerical verification of Laplace's law: the static equilibrium of a droplet 
The local curvature is calculated at each marker point x by using the osculator circle defined by x and its two neighbours. Thus, it should be exact when the interface Σ is a circle of radius R where the constant curvature equals or for a sphere in 3D. It is numerically verified up to machine precision. For an ellipse of radius a and b in the polar coordinates, the analytical value equals:
κ(θ)=ab(a2sin2θ+b2cos2θ)3/2, with θ∈[0,2π].
The secondorder accuracy in the norm is achieved with respect to the mean distance between two connected interface markers for calculating the local curvature , as shown in Fig. 3.


 Fig. 3. Secondorder space convergence rate for the local curvature of an ellipse with a =1 and b =0.75: error in L ^{2}norm versus number N of Lagrangian interface markers. Zoom 
The emblematic case of capillary effects is the static equilibrium of a cylindrical or spherical droplet [[8]]. Here, the issue is to attest to the ability of the whole numerical methodology to preserve the perfect equilibrium during an arbitrary time. Let us consider a cylindrical drop of radius R with a constant surface tension for two fluids of different densities. The computation is run with zero initial conditions, the boundary condition and the capillary force field on Σ with no gravity force. Because the curvature is , exactly computed in our case with the markers chain, the drop is likely to crash on itself with a speed equal to the capillary velocity. But since the flow is incompressible, the locally isovolume flow induces a reaction force exerted on the interface that exactly counterbalances the capillary force . Then, the static equilibrium is obtained instantaneously with a uniform capillary pressure inside the droplet which equals Laplace's law: for a disk, for a sphere and in general, whatever the densities.
The result shown in Fig. 4 is computed within one time iteration only and exhibits no parasite current. The exact uniform capillary pressure is obtained with a velocity field being exactly zero both inside and outside the drop up to machine precision.


 Fig. 4. Laplace uniform capillary pressure p _{c} =400Pa in a disk droplet of radius R =2.5⋅10^{−3}m for a constant surface tension σ =1N/m: the velocity field is zero in both cases with no parasite current. Left: unstructured mesh nonfitted to the interfacemarkers circle – Right: unstructured mesh fitted to the interface. Zoom 
 
 Numerical results for a sharp benchmark: air bubble dynamics in melted steel 
The test case of air bubble dynamics rising in liquid water, where the density ratio is , viscosity ratio with the air–water surface tension , is reputed to be difficult to compute with the suitable mesh convergence property [[11], [13]], and hence considered as already discriminating.
Here, we show that the ( ) method can compute the case proposed in [[4]] of air bubbles arising in a liquid melted steel with a density ratio , viscosity ratio and surface tension . The motion is induced by the vertical gravity force with and the capillary force . The size of the mesh is Cartesian MAC cells, Lagrangian markers are used in the interface chain, and the mesh convergence is reached. The cylindrical air bubble is initialized asymmetrically at the bottom of the liquid column in order to create a lateral motion, and we have a symmetry boundary condition for the velocity components on the walls. A typical situation is shown in Fig. 5 where the interface is plotted over the pressure or the vertical velocity fields.


 Fig. 5. Dynamics of the air bubble of diameter d =0.01m, mass density ρ _{g} =1.1768kg/m^{3}, dynamic viscosity μ _{g} =1.8510^{−5} Pa⋅s in a vertical cavity of width l =4⋅10^{−2}m and height h =10^{−1}m. The bubble rises in the liquid steel with ρ _{l} =10^{4}kg/m^{3}, μ _{l} =10^{−3}Pa⋅s and σ =1.5N/m. Left: pressure field p ∈[−9235,0]Pa (p =0 at the bottom left) at time t =0.05s – Center: vertical velocity field v _{z } ∈[−0.48,1.55]m/s and streamlines at t =0.05s – Right: Some bubble positions and shapes during time and vertical velocity field v _{z } at final time t =0.2s. Zoom 
From the physical point of view, the main conclusion of this work is that the velocity correction in a timestepping method for the incompressible flow is fully kinematic and only concerned with the Helmholtz–Hodge decomposition of the predicted velocity. Hence, this step includes no physics inside. Conversely, the pressure gradient correction does link the mass density to induce the motion and perfectly reproduce the dynamics of the phenomena. This assertion, which was only a conjecture till now, is here validated by two physical examples.