Consistent pressure Poisson splitting methods for incompressible multi-phase flows: eliminating numerical boundary layers and inf-sup compatibility restrictions

For their simplicity and low computational cost, time-stepping schemes decoupling velocity and pressure are highly popular in incompressible flow simulations. When multiple fluids are present, the additional hyperbolic transport equation in the system makes it even more advantageous to compute different flow quantities separately. Most splitting methods, however, induce spurious pressure boundary layers or compatibility restrictions on how to discretise pressure and velocity. Pressure Poisson methods, on the other hand, overcome these issues by relying on a fully consistent problem to compute the pressure from the velocity field. Additionally, such pressure Poisson equations can be tailored so as to indirectly enforce incompressibility, without requiring solenoidal projections. Although these schemes have been extended to problems with variable viscosity, constant density is still a fundamental assumption in existing formulations. In this context, the main contribution of this work is to reformulate consistent splitting methods to allow for variable density, as arising in two-phase flows. We present a strong formulation and a consistent weak form allowing standard finite element spaces. For the temporal discretisation, backward differentiation formulas are used to decouple pressure, velocity and density, yielding iteration-free steps. The accuracy of our framework is showcased through a wide variety of numerical examples, considering manufactured and benchmark solutions, equal-order and mixed finite elements, first- and second-order stepping, as well as flows with one, two or three phases.


Introduction
In numerical methods for incompressible flows, the assumption that both density and viscosity are constant is often mentioned just in passing, if even mentioned at all. It is, however, a rather "fragile" assumption that can be violated, for example, in non-Newtonian, two-phase or non-isothermal flows [1]. In fact, this violation can have severe implications in both theory and numerical practice. A variable viscosity, for instance, demands a complete reformulation of projection methods decoupling velocity and pressure in time [2][3][4]. A non-constant density, on the other hand, requires redefining even the projection operator used in such meth-B Douglas R. Q. Pacheco douglas.r.q.pacheco@ntnu.no 1 Department of Mathematical Sciences, Norwegian University of Science and Technology, Trondheim, Norway 2 Institute of Structural Analysis, Graz University of Technology, Graz, Austria ods [5]. In particular, for flows with more than one fluid phase, density and viscosity fields are not only variable, but discontinuous-which certainly brings additional challenges to the discretisation. Nowadays, there are various numerical techniques with different levels of sophistication to handle such issues. Sharp two-phase flow simulations can be achieved, for instance, with enriched [6], extended [7] or unfitted [8] finite element methods. For the temporal discretisation, there are also different methods to cope with the variable density, especially in split-step schemes [5,[9][10][11][12][13].
Split-step, time-splitting or fractional-step methods are interchangeable terms to designate time-stepping schemes decoupling the computation of velocity and pressure. Their idea is to transform an incompressible flow system into a series of Poisson and convection-diffusion subproblems that are simpler to solve. Especially for multi-phase flows, which feature conservation laws of different nature, split-step methods are more popular [5,[9][10][11][12][13][14][15][16][17][18][19][20] than monolithic schemes [21][22][23]. Among the most common time-splitting schemes are so-called projection or pressure correction methods, in which incompressibility is enforced by projecting the velocity onto a weakly divergence-free space. Such methods, however, have some well-known shortcomings. The earliest non-incremental variant, for example, induces completely non-physical pressure boundary conditions [24]. This can be partially remedied by so-called rotational variants [2,4,11], which in turn induce an inf-sup compatibility restriction [25,26] on the finite element spaces (equal-order discretisations of pressure and velocity become unstable). In fact, these issues are not exclusive to multi-phase flows, arising even in the Stokes problem.
For single-phase flows, Henshaw and Petersson [27], Johnston and Liu [28] and Liu [29] proposed an alternative family of split-step schemes overcoming the limitations of classical projection methods. In their version, explicit incompressibility is replaced by a consistent pressure Poisson equation (PPE) obtained directly from the balance of momentum. Their schemes do not only reduce the computational cost, but also break the finite element compatibility restrictions without inducing any spurious pressure boundary layers. The price to pay is a slightly more complicated right-hand side, which can however be handled with hardly any additional implementation effort. Nonetheless, constant viscosity and constant density are two strong assumptions in the classical derivations [28,29]. While we have recently extended consistent splitting methods to the case of variable viscosity [30], handling more than one phase requires also allowing for density variations, which has a stronger impact on the formulation. This is thus the main topic of the present work.
Our main goal is to devise a fully consistent pressure Poisson framework for incompressible flows with variable density, with a focus on two-phase flows. We can derive an elliptic equation for the pressure by first dividing the Navier-Stokes momentum equation by the flow density, adding a consistency term and then taking the divergence of the result. Similarly, pressure boundary conditions can be obtained from the momentum equation and its boundary data. This gives us an equation to compute the pressure from the velocity and density fields, and at the same time guarantees, on the continuous level, a divergence-free velocity-without requiring projections or corrections. This, combined with semi-implicit temporal discretisations, converts the rather complex mass-momentum system into a series of much simpler subproblems: a Poisson equation with variable coefficient, a scalar advection problem and a vector-valued convection-diffusion equation that may be decomposed into scalar components. Furthermore, when considering a finite element spatial discretisation, our schemes allow the same polynomial degree to be used for all quantities, simplifying implementation.
The rest of this article is organised in the following way. After the basic problem setup and notation are briefly described in Sect. 2, Sect. 3 presents a PPE-based reformulation that implies incompressibility without explicitly enforcing the divergence-free constraint on the velocity. When the goal is to use standard finite element elements for the spatial discretisation, we need a C 0 -compatible variational formulation, which is the topic of Sect. 4. Discretisation and stabilisation aspects are addressed in Sects. 5 and 6. Finally, several numerical examples are given in Sect. 7 to test our approach in problems with increasing complexity. Let us first remark that, differently from most works on twophase flow simulation, ours does not intend to address the challenges related to interface conditions, jumps, forces, etc. We focus instead on the time-splitting aspect of the problem, that is, on showing how one can construct simple, minimally invasive, finite-element-suitable schemes to decouple all system unknowns in time, without sacrificing boundary accuracy or inducing inf-sup compatibility restrictions. For this reason, we focus on gravity-driven flow regimes in which surface forces can be neglected. How one handles interface-related matters, however, is immaterial to the present discussion, as different techniques can be incorporated straightforwardly.

Problem setup and notation
Let us consider a spatial domain ⊂ R d , d = 2 or 3, with a Lipschitz boundary := ∂ decomposed into three nonoverlapping regions D , N and S . On the first two we prescribe velocities and normal tractions, respectively, while the latter describes a free slip boundary. Imposing slip or freeslip conditions is important in multi-phase flows to allow "wetting", since a standard no-slip condition would preclude fluid-fluid interfaces from moving along walls. Under these conditions, the standard setting for the incompressible Navier-Stokes system with non-constant density reads equipped with initial and boundary conditions where n is the outward unit normal vector on ∂ , (ρ 0 , u 0 , , v, t) are appropriate data, g denotes a gravitational field, (u, p, ρ) are the flow velocity, pressure and density, respectively, and S is the viscous stress tensor, which for Newtonian fluids takes the form with μ := ρν and ν denoting the dynamic and kinematic viscosities, respectively. The free-slip condition can be replaced by a slip condition if we appropriately modify the right-hand side in (10), see for instance Refs. [31,32]. Also note that (10) can be written in terms of tangent vectors [33]. Since the mass equation (1) is hyperbolic, boundary conditions for ρ are only prescribed on in , which is the part of ∂ where u · n < 0.

A consistent pressure Poisson reformulation
Regardless of the spatial discretisation method used, efficiently solving (1)-(10) in a fully coupled way is a challenging task. However, because of the saddle-point structure of the velocity-pressure subsystem, decoupling the unknowns is not straightforward here. Such structure also induces an inf-sup or Ladyzhenskaya-Babuška-Brezzi (LBB) compatibility restriction when discretising p and u with finite elements [25,26]. To overcome both limitations, our idea is to replace the incompressibility constraint (3) by a consistent pressure Poisson problem: in which D,S := D ∪ S , ∂(·) ∂n := n · ∇(·), β is a nonnegative parameter to be defined later (see Sect. 5.2), and with v · n = 0 on S and I denoting the d × d identity tensor.
Parameter χ in Eqs. (12) and (16) is a positive constant which we will later show that should be set as χ = 2. It is still not evident why (12)-(16) would be a suitable replacement for the usual condition ∇ · u = 0, but for now let us take a closer look at the equations. Equation (12) can be obtained by adding χμ∇(∇ · u) to the left-hand side of Eq. (2), dividing both sides by ρ, taking minus the divergence of the result and finally adding β∇ · u to the left-hand side. Both added terms are consistent with the original system, whose solution u is divergence-free. The absence of ∂ t u in Eq. (12) is also consistent, as Moreover, the Neumann boundary condition (14,16) for the pressure is, up to an additional χμ∇(∇ · u), simply the momentum equation projected onto the normal direction n and restricted to D,S . The Dirichlet boundary condition (13,15) for the pressure can be obtained by dotting the traction boundary condition (8) with n and subtracting 2μ∇ · u from the left-hand side, which is again consistent. The reason for the added terms will become clear soon, as we show the equivalence between the standard Navier-Stokes system and the PPE-modified one.
The first side of the equivalence, i.e., that (12)-(16) can be derived from the momentum-incompressibility system, has just been outlined. Slightly more involved will be proving the other side of the equivalence: does the pressure Poisson problem, combined with the momentum equation, imply incompressibility? To show that, we divide Eq. (2) by ρ, apply the divergence and add the result to the PPE (12). This leaves us with or, by introducing ϕ := ∇ · u, with the convection-diffusionreaction equation What is missing now are boundary and initial conditions for this equation. Neumann data can be obtained by dotting (2) with n, restricting the result to D,S and then subtracting (14), which gives χμn · ∇(∇ · u) = 0, that is, since χμ = 0. Similarly, to get Dirichlet data for ϕ we multiply the traction boundary condition (8) by n and add the result to (13), yielding 2μ∇ · u = 0, or simply Thus, with zero Dirichlet and Neumann boundary conditions, and provided that the partial differential equation (17) has ϕ ≡ 0 as its unique solution, that is, ∇ · u = 0 in × [0, T ], as we wanted to show. In summary, combining the momentum equation (2) with the pressure Poisson problem (12)-(16) implies incompressibility without the need to enforce ∇ · u = 0 explicitly or to project u onto a divergence-free space. The main advantage of this PPE-based reformulation is that we now have an invertible operator for the pressure, as opposed to the original saddle-point structure. From the spatial discretisation standpoint, this eliminates the LBB compatibility restriction and allows equal-order finite element spaces [30]. For the temporal discretisation, the PPE gives us an invertible equation to compute p from (ρ, u) without resorting to the artificial boundary conditions typical of projection methods.
The pressure Poisson equation, as presented in (12), has third-order derivatives on the velocity. Fortunately, we can eliminate them by using certain identities. Notice that Now, by choosing χ = 2, we can get rid of the third-order derivatives: since ∇ · (∇ × v) = 0 for any vector v. We are then left with only first-and second-order derivatives in (22). Now, the purpose of the χν∇(∇ · u) term artificially (but consistently) added in Eq. (12) becomes clear, and it is in fact twofold: it adds diffusion to the divergence equation (17) and allows us to eliminate the third-order derivatives in the PPE by simply setting χ = 2. Although both consistency and incompressibility hold for any χ > 0, this is not intended here as a tunable parameter. We could in principle choose χ = 2, but this would only leave us with unnecessary higher-order derivatives. Since χ features only in the PPE, it adds no diffusion to the momentum equation, but only to Eq. (17), whose unique solution is ∇ · u = ϕ ≡ 0 for any positive χ .

Weak pressure Poisson formulation
The second-order velocity derivatives in (22) cannot be handled by standard Lagrangian finite elements, so we need a weak formulation containing only first-order derivatives. This is not straightforward, and in fact Li et al. [13] used finite differences to evaluate the higher-order derivatives at the nodal points, then interpolated those values in a continuous source term-which is hardly practical in unstructured meshes. To allow a more general approach suitable for C 0 finite elements in arbitrary meshes, we will next show how to reduce the order of differentiation, starting from (21) and using some vector calculus relations. We first multiply the PPE (12) by a continuous test function q with zero value on N . Then, we integrate by parts all terms but β∇ · u, and enforce the Neumann boundary condition (14) to get where we can further use integration by parts to write but (24) and Moreover, Therefore, collecting the results from (23)- (26), we arrive at a weak formulation: find p satisfying the Dirichlet condition p = ζ on N (13), such that for all test functions q being zero on N . Now that we are left with only first-order derivatives, C 0 finite elements can be used to discretise all unknowns and test functions.

Discretisation
This section will address discretisation aspects in both space and time. Although we focus on finite elements for the spatial discretisation, the fractional-step schemes presented here are also valid for other frameworks such as finite differences and finite volumes.

Temporal discretisation
As usual in split-step methods, we employ backward differentiation formulas (BDFs) to approximate the temporal derivatives, combined with matching extrapolation rules for linearisation. The idea is to have iteration-free schemes completely decoupling velocity, pressure and density. For concision, we shall limit the presentation to first-and second-order temporal discretisations, but the schemes can be straightforwardly extended to higher-order versions with possibly adaptive stepping (see Ref. [30] for details). Let t > 0 be a finite time-step size and υ n denote a certain quantity υ at time t = t n . For the first-order case, the backward finite difference BDF1 and corresponding extrapolation read: whereas to go one order higher we use BDF2: For the mass equation, this results in for BDF1, or (31) for BDF2. Either way, we are left with a linear partial differential equation that still needs to be discretised in space.

Spatial discretisation
It is well-known that the Galerkin finite element method, in its standard form, is not ideal for problems with low or zero diffusion, such as (30) or (31). Therefore, when solving the mass transport, it is common practice to employ residual-based [32,34] or entropy [11] stabilisation techniques. However, we can leverage the absence of second-order derivatives to use a least-squares (instead of Galerkin) finite element formulation, as done by Pyo and Shen [10]. For an equation of the form it consists of seeking u ∈ X fulfilling inflow conditions, such that with X being an appropriate finite element space. The leastsquares finite element method is suitable for hyperbolic problems without requiring stabilisation terms or tunable factors. The other two equations in our system-momentum and pressure Poisson-on the other hand, can be tackled by a standard Galerkin finite element discretisation. Yet, some stabilisation may be required for high-Reynolds flows, as addressed in Sect. 6.4. For concision, we shall state the weak momentum equation considering S = ∅; details on the implementation of (free-)slip conditions can be found for instance in Refs. [31,35]. The following sets are considered: where subscript h denotes the finite element interpolation of a certain quantity or geometry, and S k h ( h ) denotes the standard space of continuous, piecewise polynomial shape functions of degree k in the triangulation of h . Notice that, in general, we will formally have R m (t), since the inflow boundary depends instantly and locally on n · u(t).
Another discretisation-related matter to be discussed is the choice of the weighting parameter β in the pressure Poisson equation (27). On the continuous level, all we need is β ≥ 0 to guarantee incompressibility, otherwise a negative β could spoil the unique solvability of (17) with non-zero eigenmodes of ϕ = ∇ · u. On the discrete level, however, β should be chosen appropriately if optimal convergence is desired. First, notice that the term β∇ · u penalises large violations of discrete incompressibility, so that β works as a penalty-like parameter. While Li [36] proposed an expression based solely on the spatial resolution, Pacheco et al. [37] showed that in fact β can be taken as the inverse of the well-known PSPG parameter: in which h denotes the element size and (α 1 , α 2 , α 3 ) are tunable parameters. For stationary problems we take α 3 = 0, whereas α 1 = α 2 = 0 when using LBB-compatible spaces (e.g. Taylor-Hood elements). To reduce the number of userdefined parameters, we set where α is a single tunable parameter, k and l are the orders of the velocity and pressure discretisations, respectively, and δ kl refers to the Kronecker delta.

The Dirichlet boundary condition for the pressure
We now briefly discuss a technical issue involving the pressure Dirichlet boundary condition (13). For problems with N = ∅, we have where t is some given normal traction. Since on the discrete level the stress tensor S = 2μ∇ s u will be discontinuous, so will ζ across elements. This is an issue when using continuous finite elements for the pressure, but there are various ways to circumvent that. The most common approach is to project ζ (in the L 2 ( N ) sense) onto a continuous space [4,29,30].
Another possibility is to notice that, although n can be discontinuous on non-planar surfaces, outflow boundaries are normally generated by truncating a larger physical domain, so that each outlet portion of the computational domain is planar. That being the case, the tangential projector (n ⊗ n − I) will be constant-or piecewise constant, if there are multiple outlets-so that we could simply perform an L 2 ( ) projection of S onto a continuous finite element space. An L 2 ( ) projection is simpler than an L 2 ( N ) projection, but also more expensive. Yet, we can pre-compute and store the LDL decomposition of the mass matrix used in the projection, or even lump it [30]. As we will soon see, projecting the stress tensor may also be useful in other parts of the solver. It is also possible to compute nodal values by simply averaging ζ on the boundary patch surrounding each node, although we have not tested this option. Regardless of the approach used, we shall henceforth consider that a continuous valueζ is used.

Fractional-step schemes
The main strength of split-step methods is decoupling pressure from velocity. With PPE-based schemes, the idea is to march in time by first computing the current velocity from previous pressures, then using the updated velocity (and density) to compute the current pressure, and so on for the next steps. To do that, the pressure term in the momentum equation is treated explicitly. On the other hand, when updating the pressure we can treat all terms implicitly in the PPE, since both the current velocity and the current density will already be known at that point.

Initialisation
In general, the incompressible Navier-Stokes system is not equipped with initial conditions for the pressure. This can be an issue depending on the temporal discretisation considered, as discussed for instance by Rang [38]. In the pressure Poisson framework, on the other hand, we naturally compute this initial value by seeking p 0 ∈ Q lζ 0 such that for all q ∈ Q l 0 . This is why it is crucial, in deriving the PPE, to divide the momentum equation by ρ before taking the divergence. If we were to apply ∇· directly to the momentum equation, the transient term would result in ∇ · (ρ∂ t u) = ρ∂ t (∇ · u) + ∂ t u · ∇ρ = ∂ t u · ∇ρ ≡ 0 , instead of ∇ · (∂ t u) = ∂ t (∇ · u) ≡ 0. Consequently, the integral on D in Eq. (36) would be replaced by a volume integral depending on the acceleration ∂ t u. In that case, initialising the time-marching would not be straightforward, since the initial acceleration is not given for flow equations. That would require using a forward discretisation of ∂ t u, resulting in an implicit step. Our formulation, on the other hand, requires only the time derivative of the normal data v · n on D , which may sound similar to the situation just described, yet is fundamentally different. Since v is given, its time derivative is either known exactly (e.g. when artificially ramping up inflow data), or can be approximated by a forward finite difference-which does not lead to an implicit coupling, since v·n is known a priori at all times. This discussion is applicable also to the initial-boundary value problem in strong form, hence it holds regardless of the spatial discretisation. Thus, using only given quantities, we are able to compute p 0 and start marching in time.

First-order stepping
As mentioned in Sect. 5.1, our first-order stepping combines BDF1 with linear extrapolation (u n+1 ≈û n+1 = u n and p n+1 ≈p n+1 = p n ). Given (ρ n , u n , p n ), we march in time by executing the following steps:

Viscosity update
Compute ν n+1 based on the updated phase distribution described by ρ n+1 .

Velocity update
Find for all w ∈ W k 0 .

Pressure Poisson step
Find p n+1 ∈ Q lζ n+1 such that for all q ∈ Q l 0 .

Viscosity update
Compute ν n+1 based on the updated phase distribution described by ρ n+1 .

Velocity update
Find for all w ∈ W k 0 .
In both schemes, the pressure Poisson step is fully implicit, that is, the current pressure is computed directly from the current velocity and density fields. Similarly, there is no need to extrapolate the density when updating the velocity, since we already know ρ n+1 when computing u n+1 . If the order of the steps is inverted, with u updated before ρ, we need to replace ρ n+1 byρ n+1 in the velocity step, but on the other handû n+1 could be replaced by u n+1 in the convective step.

Convection-dominated flows
At high Reynolds numbers, the Galerkin finite element approximation of the momentum balance suffers from convective instabilities, unless a very fine mesh is used. A simple remedy for that is the streamline upwind/Petrov-Galerkin (SUPG) stabilisation. On the other hand, dominant convection also means that we can treat the viscous term explicitly (or semi-implicitly) without incurring a parabolic time-step restriction [28]. In the weak form of the momentum equation, the viscous term reads The last term has a coupling effect over the velocity components, so that treating it explicitly decouples the velocity update step into d scalar convection-diffusion equations (see Ref. [30] for details). This reduces computational costs and allows us to apply SUPG individually to each of the scalar equations [39]. Let us, for concision, write the ith component of u at t = t n+1 as u i , andυ n+1 asυ for any extrapolated quantity υ. In that case, our semi-discretised balance of momentum in the x i direction reads with α j denoting the coefficients of the finite difference discretisation (α 0 = −α 1 = 1/ t for BDF1 or α 2 = −α 1 /4 = α 0 /3 = 0.5/ t for BDF2). For the SUPG method, we add to the left-hand side of the weak momentum equation the nonlinear form where τ SUPG is a mesh-dependent parameter and w is a scalar test function. Notice that this stabilisation is in general consistent, as the added term is proportional to the residual of Eq. (42). However, this is not the case when using linear elements, which are unable to approximate the second-order derivatives in the viscous term. Although this is in principle not an issue for stability, the missing terms may lead to inaccuracy on flow regions where diffusion dominates. Therefore, for first-order elements we propose replacing (43) by whereŜ i j denotes the components of the L 2 ( ) projection of the viscous stress tensor S onto a continuous finite element space. In other words, we findŜ i j ∈ X h such that for all σ ∈ X h , in which X h can be taken, for instance, as the pressure finite element space (without boundary conditions). Projecting the viscous stress tensor onto a continuous finite element space is also usual in certain stabilisation [40,41] and projection methods [3]. Since S is symmetric, the "cost" of this projection are d(d + 1)/2 scalar mass-matrix systems, but there are two simple ways to basically eliminate the overhead. One is to store the LDL decomposition of the mass matrix during initialisation and re-use it for all time steps; alternatively, the mass matrix can be lumped, which is our approach herein. As for the stabilisation parameter, we use which is simply one among several possibilities [42].

Numerical examples
In this section, we assess the accuracy of our fractional-step framework through test cases with increasing complexity. One of the main advantages of our approach with respect to projection methods is allowing equal-order interpolation of all unknowns. Regarding the velocity-pressure pairing, we consider linear-linear ( P 1 P 1 ) and Taylor-Hood ( P 2 P 1 ) triangular elements; for numerical examples with other element classes, including three-dimensional ones, refer to our recent work on split-step schemes for non-Newtonian flows [30].
In all examples, a first-order discretisation (P 1 ) is used for the density.

Manufactured solution
We start with a simple problem having an analytical solution and a pure Dirichlet boundary. In the square domain = (0, 1) 2 , we consider the flow described by The accuracy of the first-order ( P 1 P 1 P 1 , BDF1) and second-order ( P 2 P 1 P 1 , BDF2) schemes is assessed through the spatial error at t = T = 1. Divergence penalty is not used in this example, i.e., we set β = 0 in the PPE. The velocity error is measured in the relative H 1 ( ) semi-norm, whereas for pressure and density we use the relative L 2 ( ) norm. Starting with a coarse orthogonal mesh of four triangles, seven levels of uniform spatial refinement are applied, with the time-step size also halved at each level via t = h/2. The results of the convergence study are shown in Figure 1. As expected, the first-order elements yield linear convergence, whereas the Taylor-Hood elements converge quadratically when combined with BDF2. Notice, however, that the exact solution for this problem is smooth, so that optimal convergence is attainable with standard discretisations. In flows with more than one phase, on the other hand, the spatial convergence will in general not be optimal unless velocity kinks and possible pressure discontinuities are appropriately resolved [7].

Rayleigh-Taylor instability
One of the most popular benchmarks for variable-density flow solvers is the Rayleigh-Taylor instability with regularised interface [5,43]. The standard setup considers two fluids initially at rest in = (0, a/2) × (−2a, 2a) and sub-(a) (b) Fig. 1 Uniform refinement study considering triangular elements and a problem with smooth exact solution ject to g = (0, −g) . Fluid 2 sits initially on top of fluid 1 and is denser (ρ 2 > ρ 1 ) but equally viscous (μ 1 = μ 2 ). The initial transition is regularised via where ι(x) = −0.1a cos(2π x/a) describes the interface line at t = 0. We nondimensionalise the problem through the following reference quantities: ν 1 for viscosity, ρ 1 for density, a for length and √ ag for velocity, which gives us τ = √ a/g as temporal scale. Then, the problem is fully parametrised through the Reynolds and Atwood numbers:  which are set as 1000 and 0.5, respectively. For the simulation, we use the first-order scheme with α = 0.01, t = τ/1000 and a uniform mesh with 6.4 × 10 5 elements. No SUPG stabilisation is used, which is possible due to the relatively fine mesh. The top and bottom walls are no-slip boundaries, while symmetry is enforced on the sides. showing good agreement with other numerical [22] and experimental [44] results Figure 2 shows the density field for different time instants, agreeing well with the results reported by Guermond and Salgado [11] (mind that they use a different time scale τ = τ √ At). For a quantitative comparison, Figure 3 shows how the height H of the rising bubble evolves in time. Our results agree well with the simulations reported by Tryggvason [43] and Guermond and Quartapelle [5].

Dam break
We now move to a more challenging setup involving air and water. The dam break experiment by Martin and Moyce [44] is a popular benchmark for two-phase [22,45] and free-surface [46] flow simulations. This is a high-Reynolds, gravity-driven flow and, as such, suits our application scope very well. The setup considered by Landet et al. [22] consists of a tank = (0, 5a) × (0, 3a) containing a rectangular water column (0, a) × (0, 2a) surrounded by air, initially at rest, with free-slip walls and a = 57.15 mm. We set ρ air = 3ρ water /2500 = 1.2 kg/m 3 and ν air = 15ν water = 1.5 × 10 −5 m 2 /s as material properties, and g = (0, −9.81) m/s 2 for gravity.
Due to the large density difference between air and water, for this example we replace the explicit mass convection (1) by the level-set equation with φ being a signed distance function: φ > 0 corresponds to air, φ < 0 to water, and φ = 0 to the interface. We can therefore write ρ = ρ water + ρ air 2 + ρ air − ρ water 2 sign φ , (b) Grid used for both P 2 P 1 P 1 and P 1 P 1 P 1 discretisations.

Fig. 6
Three-phase channel flow: setup and grid and analogously for the kinematic viscosity. Since both ∇ρ and ∇ν appear in the PPE (27), it is useful to regularise the discontinuity. A common approach is to replace the sign function by a hyperbolic tangent [47], which is also done for this example: An excellent discussion on different variants of the level-set method is provided by Olsson and Kreiss [48]. Most importantly, our formulation remains consistent regardless of how the density is treated, as the PPE is derived solely from the momentum equation.
To handle the dominant convection, especially in the aerodynamic flow, we also use SUPG stabilisation. The discretisation is again first-order, with 3 × 10 5 triangular elements and t = 10 −4 s. The main benchmark quantity here is the temporal evolution of the water column's height, H . We compare our solution to the measurements by Martin and Moyce [44] and the numerical results by Touré et al. [45], who used a modified level-set method with corrected mass conservation. As shown in Figure 4, all three results are in very good agreement.
To illustrate the complex flow emerging from this simple example, we show in Figure 5 the Euclidian norm of the velocity field for different instants in time. The downward motion of the water column "propels" the air upwards and, as the air goes past the column's corner, a vortex is created as in a backward-facing step flow. As the water column keeps going down, two adjacent air vortices can be seen. We also observe smaller vortices around the water surge front as it shoots upwards after meeting the right wall. In the same figure, the density plot reveals quite a sharp interface, in spite of the regularisation.

Three-phase channel flow
As our last test case, we investigate a three-phase flow in the straight channel = (0, L) × (0, H ). The three fluids are initially at rest, stacked up as illustrated in Figure 6 (a). They all have the same kinematic viscosity ν, but different dynamic viscosities since ρ 3 = 1.5ρ 2 = 3ρ 1 = 1.5. For a constant pressure gradient ∂ p ∂ x = k, we can calculate the analytical solution for the stationary Poiseuille-like profile: For concision, we omit the derivation, but it can be done in the standard way: solving the stationary momentum equation and enforcing continuity of velocities and viscous stresses at the interfaces. The main idea with this example is to illustrate how the penalty-like term β∇ · u (cf. Eq. (12)) can improve mass conservation on the discrete level. For this, we consider a test case with prescribed inflow, but an open outflow. When reaching a stationary solution, we would ideally like to see the same profile on the open outlet (x = L) as the one prescribed on the inlet (x = 0). So, we enforce the inflow profile u = (u(y) f (t), 0) , with ramping up smoothly from zero (for t ≤ 0) to unit value (for t ≥ t ). For the density, the inlet boundary condition is constant in time, obeying the initial phase distribution. Since prescribing zero (or constant) normal traction t does not allow reproducing a developed flow, we swap the traction boundary condition (8) by a pseudo-traction condition on the outflow boundary: (μ∇u)n − pn =t , witht being some normal pseudo-traction data, usually taken as zero in flows with one single outlet. To naturally impose the pseudo-traction, we switch from a stress-divergence formulation of the momentum equation to a generalised Laplacian approach [49]. In that case, the Dirichlet value (15) for the pressure has to be slightly adapted: Details on outflow boundary conditions and the generalised Laplacian formulation can be found in Refs. [30,49]. Then, by settingt = 0, g = 0, T = 2.5t = 1, L = H = 1, ν = 1 and k = −1, we have the complete problem setup. Fig. 7 Three-phase channel flow: solutions obtained with (α > 0) and without (α = 0) divergence penalty. The density field, on the right, shows no remarkable differences between all the schemes. On the left, the (prescribed) inflow and (computed) outflow profiles, normalised bȳ u := |0.5k H 2 /μ 1 |, showcase the improved mass conservation for α = 0.002 in the first-order scheme Fig. 8 Three-phase channel flow: how the divergence error and the density under-and overshoot are affected by the divergence penalty parameter α, using P 1 P 1 P 1 elements For the temporal discretisation, this time we use BDF2 with t = 0.02. If the interfaces were to lie completely on grid lines, the velocity kinks would be resolved by the spatial discretisation; to avoid that and be more realistic, we use the non-orthogonal mesh shown in Figure 6 (b). The first two test cases employ P 2 P 1 P 1 and P 1 P 1 P 1 elements, respectively, without any divergence penalty (α = 0, cf. Eq. (35)), whereas the third one has α = 0.002 and also first-order elements. Figure 7 depicts the density field and the velocity profiles at t = 1, when the developed steady-state solution is already established. The density distribution is similar in all cases, exhibiting interface oscillations due to the relatively coarse mesh for a multi-phase flow. In spite of that, the velocity profile shows no apparent oscillations. The plots reveal how even a small α can virtually eliminate the mass loss otherwise seen for the P 1 P 1 P 1 elements, whereas the second-order discretisation shows no relevant loss even without any penalty. We must remark, however, that finding a good value for α requires some tuning: setting α too large can lead to instability, by giving excessive weight to the velocity divergence in comparison to the pressure-stabilising term −∇ · (ρ −1 ∇ p). Figure 8 shows how α affects the divergence error and the density over-and undershoot at t = T . The density results show hardly any difference, as previously indicated in Figure  7. The divergence error, on the other hand, reduces by 34% when we change α from 0 to 1.4 × 10 −3 . While sensitivity analysis is out of the present scope, a detailed study on how this parameter influences the accuracy of different quantities is found in our previous work for single-phase flows [37].

Concluding remarks
In this work, we have devised, implemented and thoroughly tested a consistent split-step framework for flows with nonconstant density and viscosity. Most splitting schemes for multi-phase flows are based on pressure correction methods enforcing incompressibility via solenoidal projections. Those methods are very efficient but have well-known shortcomings, such as numerical boundary layers and/or spatial stability restrictions when using finite elements. Our method, on the other hand, is based on a consistent pressure-Poisson reformulation of the Navier-Stokes system, thereby avoiding those issues without incurring computational overhead. Moreover, the method is minimally invasive: it decomposes the variable-density Navier-Stokes system into simpler convection, diffusion and convection-diffusion problems, and the pressure Poisson step is consistent regardless of the technique-if any-applied to handle discontinuities. No spurious pressure boundary layers or finite element compatibility restrictions are induced.
As a matter of fact, consistent PPE-based schemes have recently been used in designing accurate, competitive solvers for challenging applications such as non-Newtonian flows, fluid-structure interaction and free-surface flows [30,36,[50][51][52]. Of course, a consistent pressure Poisson equation has a more complex right-hand side than a Leray projection does, but our analytical efforts in Sects. 3 and 4 ensure that this complexity does not translate into an implementation overhead. Indeed, in the end we get a weak formulation with the exact same left-hand side as in classical projection methods [5], and a right-hand side allowing standard finite element spaces-mixed or equal-order. Although we have focused on finite element methods, our PPE-based system can be tackled within basically any other framework. It is also straightforward to consider generalised Newtonian fluids, since our equations already account for variable viscosity.
What our formulation does not improve on, with respect to projection methods, is the pressure operator. In both cases we have −∇ · (ρ −1 ∇ p), which is not optimal in the presence of pressure discontinuities, as seen e.g. in surface tension mod-els. In particular, piecewise constant finite element spaces are not allowed for the pressure, which is a recurring limitation in split-step methods. To overcome that, we are currently working on an ultra-weak reformulation of the problem, with promising results so far [53].
Funding Open access funding provided by NTNU Norwegian University of Science and Technology (incl St. Olavs Hospital -Trondheim University Hospital).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.