Access to the PDF text

Free Article !

Comptes Rendus Mathématique
Volume 354, n° 11
pages 1124-1131 (novembre 2016)
Doi : 10.1016/j.crma.2016.06.007
Received : 11 May 2016 ;  accepted : 27 June 2016
A kinematic vector penalty–projection method for incompressible flow with variable density
Une méthode cinématique de pénalité–projection vectorielle pour l'écoulement incompressible à densité variable

Philippe Angot a , Jean-Paul Caltagirone b , Pierre Fabrie c
a Aix-Marseille Université, Institut de mathématiques de Marseille – CNRS UMR 7373, Centrale Marseille, 39, rue Frédéric-Joliot-Curie, 13453 Marseille cedex 13, France 
b Université de Bordeaux & IPB, Institut de mécanique et d'ingénierie de Bordeaux – CNRS UMR 5295, 16, avenue Pey-Berland, 33607 Pessac, France 
c Université de Bordeaux & IPB, Institut de mathématiques de Bordeaux – CNRS UMR 5251, ENSEIRB–MATMECA, 33400 Talence, France 


In this Note, we present a new version of the vector penalty–projection splitting method described in [[1]] for the fast numerical computation of incompressible flows with variable density and viscosity. We show that the velocity correction can be made completely independent of the mass density ρ . Hence, this step is purely kinematic using the fast Helmholtz–Hodge decompositions proposed in [[2]]. Then, it is shown that the dynamic step of pressure gradient correction can be fast and locally consistent on edge-based generalized MAC-type unstructured meshes that naturally verify the compatibility condition in the proposed discrete setting. By the way, a new accurate front-tracking Lagrangian-advection technique is also introduced for multiphase flows.

This new method preserves the fully vector formulation of both the prediction and correction steps of the original scheme, the primary unknowns being   and ρ by advection, since the pressure Neumann–Poisson problem remains eliminated. The efficiency of the present method is demonstrated through numerical results on sharp test cases.

The full text of this article is available in PDF format.

On présente dans cette Note une nouvelle version de la méthode de splitting par pénalité–projection vectorielle décrite dans [[1]] pour le calcul des écoulements incompressibles à masse volumique et viscosité variables. Le principal résultat est de rendre la correction vectorielle de vitesse complètement indépendante de la masse volumique ρ . Cette étape devient donc purement cinématique et correspond à une décomposition rapide de Helmholtz–Hodge proposée dans [[2]]. On montre que l'étape dynamique de correction du gradient de pression peut être rapide et localement consistante sur des maillages généralisés de type MAC non structurés.

The full text of this article is available in PDF format.
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 low-Mach 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 multi-physics problems. Many refined techniques have been developed for the front-tracking or front-capturing of moving and deformable phase interfaces, e.g., [[10], [12], [13], [14]], and most of them use the very popular time-stepping 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 edge-based generalized MAC-type 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 time-space 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 semi-discrete scheme ( )

The fast ( ) method, proposed here using a first-order 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) Divergence-free 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 front-tracking of density (see Sect. 2.3):ρn+1−ρnδt+vn+1⋅∇ρn=0 in Ω. In the sub-step (d ), the scalar potential   is theoretically equal to   up to an additive constant, as a by-product of the (VPP) sub-steps  . However, to avoid round-off 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 edge-linking pressure vertices.

Since   with (c ), it is important to notice that the sub-steps (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 sub-step (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)(K-VPPε){ρ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 sub-step (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 well-posed 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 first-order scheme with  , as numerically verified in [[1]].

Moreover, the velocity–pressure coupling in ((2)) can be interpreted as a new two-step 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 MAC-type mesh

The main important issue of the present ( ) method ((2)) is to show the existence at each time step of the so-called inertial density  . Indeed, since the irrotational component   of   is calculated by the (VPP) sub-step (b ), this also gives easily as a by-product 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 simply-connected, 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 edge-based generalized MAC-type 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

Fig. 1. 

Edge-based generalized MAC-type unstructured mesh. Topology of the 3-D 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 ].


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 edge-based generalized unstructured staggered mesh of MAC-type [[5], [7]] in the context of dispersed multiphase flows. The edge-based 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 2-D 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 mid-point 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 second-order 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

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.


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 surface-markers based on triangles in 3-D). Each marker point belongs to a single cell in the primal mesh and this association is made by a fast algorithm of ray-tracing [[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=(ρa|c−a|+ρb|b−c|)vˆ⋅t=(ρa|c−a|+ρb|b−c|)|b−a|∫abvˆ⋅tdx=(αρa+(1−α)ρb)(ϕb−ϕa), with α:=|c−a||b−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) sub-step (b ) in ((2)) by the standard technique for the divergence-free projection of   consisting in solving for ϕ the Neumann–Poisson problem:   with   on Γ.

A new front-tracking Lagrangian-advection method

A new simple but accurate Lagrangian front-tracking technique is now introduced for the advection of interface Σ in the sub-step ((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 Runge-Kutta explicit scheme (RK2 or RK4 with the K-VPP method of second-order in time):
c )
Detect the cells in the primal mesh that are crossed by the updated marker chain with a ray-tracing 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 time-step δ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 3-D. 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 second-order 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

Fig. 3. 

Second-order 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.


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

Fig. 4. 

Laplace uniform capillary pressure p c =400Pa in a disk droplet of radius R =2.510−3m for a constant surface tension σ =1N/m: the velocity field is zero in both cases with no parasite current. Left: unstructured mesh non-fitted to the interface-markers circle – Right: unstructured mesh fitted to the interface.


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

Fig. 5. 

Dynamics of the air bubble of diameter d =0.01m, mass density ρ g =1.1768kg/m3, dynamic viscosity μ g =1.8510−5 Pas in a vertical cavity of width l =410−2m and height h =10−1m. The bubble rises in the liquid steel with ρ l =104kg/m3, μ l =10−3Pas 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.



From the physical point of view, the main conclusion of this work is that the velocity correction in a time-stepping 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.


Angot P., Caltagirone J.-P., Fabrie P. A fast vector penalty–projection method for incompressible non-homogeneous or multiphase Navier–Stokes problems Appl. Math. Lett. 2012 ;  25 (11) : 1681-1688 [cross-ref]
Angot P., Caltagirone J.-P., Fabrie P. Fast discrete Helmholtz–Hodge decompositions in bounded domains Appl. Math. Lett. 2013 ;  26 (4) : 445-451 [cross-ref]
Angot P., Caltagirone J.-P., Fabrie P. Analysis for the fast vector penalty–projection solver of incompressible multiphase Navier–Stokes/Brinkman problems Numer. Math. 2016 ; submitted for publication. hal-01194345
P. Angot, J.-P. Caltagirone, P. Fabrie, A spectacular solver of low-Mach multiphase Navier–Stokes problems under strong stresses, in: Turbulence and Interactions, Proceedings of the TI 2015 Int. Conference, in: Notes on Numerical Fluid Mechanics and Multidisciplinary Design, Springer, in press, hal-01248744.
Blair Perot J. Discrete conservation properties of unstructured mesh schemes Annu. Rev. Fluid Mech. 2011 ;  43 : 299-318
Bossavit A. Computational Electromagnetism  San Diego, CA, USA: Academic Press (1998). 
Caltagirone J.-P. Discrete Mechanics  London: ISTE Ltd and J. Wiley & Sons (2015). 
Caltagirone J.-P. Modélisation des effets capillaires en mécanique des milieux discretshal-012516702016
Caltagirone J.-P., Vincent S. A kinematics scalar projection method (KSP) for incompressible flows with variable density Open J. Fluid Dyn. 2015 ;  5 : 171-182 [cross-ref]
Harlow F.H., Welch J.E. Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface Phys. Fluids 1965 ;  8 (12) : 2182-2189 [cross-ref]
Hysing S., Turek S., Kuzmin D., Parolini N., Burman E., Ganesan S., Tobiska L. Quantitative benchmark computations of two-dimensional bubble dynamics Int. J. Numer. Methods Fluids 2009 ;  60 : 1259-1288 [cross-ref]
Sarthou A., Vincent S., Caltagirone J.-P., Angot P. Eulerian–Lagrangian grid coupling and penalty methods for the simulation of multiphase flows interacting with complex objects Int. J. Numer. Methods Fluids 2008 ;  56 (8) : 1093-1099 [cross-ref]
Tryggvason G., Scardovelli R., Zaleski S. Direct Numerical Simulations of Gas-Liquid Multiphase Flows  Cambridge, UK: Cambridge University Press (2011). 
Wang Y., Simakhina S., Sussman M. A hybrid level set-volume constraint method for incompressible two-phase flow J. Comput. Phys. 2012 ;  231 : 6438-6471 [cross-ref]

© 2016  Académie des sciences@@#104156@@
EM-CONSULTE.COM is registrered at the CNIL, déclaration n° 1286925.
As per the Law relating to information storage and personal integrity, you have the right to oppose (art 26 of that law), access (art 34 of that law) and rectify (art 36 of that law) your personal data. You may thus request that your data, should it be inaccurate, incomplete, unclear, outdated, not be used or stored, be corrected, clarified, updated or deleted.
Personal information regarding our website's visitors, including their identity, is confidential.
The owners of this website hereby guarantee to respect the legal confidentiality conditions, applicable in France, and not to disclose this data to third parties.
Article Outline