{'fraction_non_alphanumeric': 0.04878730691270008, 'fraction_numerical': 0.01735747070198737, 'mean_word_length': 4.122568325575845, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Eliminating Spurious Velocities with a Stable Approximation of Viscous Incompressible Two-Phase Stokes Flow

We present a parametric finite element approximation of two-phase flow. This free boundary problem is given by the Stokes equations in the two phases, which are coupled via jump conditions across the interface. Using a novel variational formulation for the interface evolution gives rise to a natural discretization of the mean curvature of the interface. In addition, the mesh quality of the parametric approximation of the interface does not deteriorate, in general, over time; and an equidistribution property can be shown for a semidiscrete continuous-in-time variant of our scheme in two space dimensions. Moreover, on using a simple XFEM pressure space enrichment, we obtain exact volume conservation for the two phase regions. Furthermore, our fully discrete finite element approximation can be shown to be unconditionally stable. We demonstrate the applicability of our method with some numerical results which, in particular, demonstrate that spurious velocities can be avoided in the classical test cases.


Introduction
It is well-known that non-physical velocities can appear in the numerical approximation of two-phase incompressible flows. These so-called spurious currents appear in different representations of the interface, with and without surface tension. If surface tension effects are taken into account, a jump discontinuity in the pressure results, and this poses serious challenges for the numerical method. As the pressure jump is balanced by terms involving the curvature of the unknown interface, it is necessary to accurately approximate the interface, its curvature and the pressure. In this paper we introduce a new stable parametric finite element method with good mesh properties, which leads to an approximation of the interface, and the conditions on it, with the property that undesired spurious velocities are either small or vanish completely. Although Bänsch, [3], proved a stability result for a (highly nonlinear) parametric discretization of the Navier-Stokes equations with a free capillary surface, to our knowledge, our stability result, for a linear scheme, is the first in the literature for a parametric discretization of two-phase flow.
When approximating two-phase flows, one has to decide first of all on how to represent the interface. The most direct choice is an explicit tracking of the interface. In these tracking methods the interface is either triangulated or represented by a connected set of particles, which carry forces.
The interface is then transported using the flow velocity. Variants of these approaches are called front tracking methods, immersed interface methods or immersed boundary methods, see e.g. [39,27,3,30,19] for details.
In a second completely different approach the interface is captured implicitly by defining a function on the whole domain. In the volume of fluid (VOF) method the characteristic function of one of the fluid phases is used in order to evolve the interface, see e.g. [24,33,31]. In the level set method, instead of a characteristic function, the interface is represented as the zero level set of a smooth function, see e.g. [37,36,29,22] for this approach. Finally, in the phase field approach, instead of a sharp interface description, the interface is considered to be diffuse with a small interfacial layer in which a phase field variable rapidly changes from the different constant values in the two phases, see e.g. [1,25,16].
Spurious velocities have already been observed in the numerical approximation of one-phase incompressible fluid flow with external forces, see [20], where a projection method to deal with this phenomenon is also proposed. In two-phase flow with surface tension effects it is well-known that the imbalance between the discrete computation of the curvature of the interface and the pressure jump at the interface can create spurious velocity fields near the interface, even in situations where the exact solution has zero velocity. Several methods to compute the discrete curvature and different choices of enriching the pressure space have been introduced to reduce spurious velocities, see e.g. [32,33,26,17,38,18,23,41,2].
In this paper we propose a numerical method for two-phase incompressible Stokes flow based on a parametric representation of the interface. We use finite elements to approximate 1 Figure 1: Fitted and unfitted bulk finite element meshes for a circular interface. the velocity and the pressure, and the interface is approximated using a lower dimensional mesh. Here we allow both for a fitted approach, where the bulk mesh is adapted to the interface, and an unfitted approach, where the bulk and interface meshes are totally independent. Typical meshes for both approaches are shown in Figure 1. In this paper we use the unfitted approach for our numerical results, which means that we can avoid remeshings of the bulk mesh at every time step. As discussed for example in [18], unfitted bulk meshes for twophase flow with pressure jumps, due to surface tension effects, lead to a poor approximation of the pressure. We avoid this by using locally refined bulk meshes at the interface in practice. But we stress that all our theoretical results presented in this paper also hold for fitted meshes. Another strategy in the context of the unfitted approach, which can be combined with our proposed method, is to enrich the pressure finite element space by functions providing additional degrees of freedom close to the interface; and, in particular, allow for pressure jumps in the elements cut by the interface. Such an enrichment technique is an example of the extended finite element method (XFEM), and has been used for two-phase problems in the context of level set methods. A major drawback of this approach is that the resulting linear algebraic system is typically very ill-conditioned, because the linear independence of the finite element basis deteriorates as enrichment functions with very small support arise. Although strategies have been developed to reduce the problem of ill-conditioning, see e.g. [23,2,34], the computational effort still remains large due to the reconstruction of the XFEM basis as the interface moves. An XFEM approach is also possible within the context of our method, and will be discussed later.
A common problem in approaches which directly parameterize an evolving interface is that typically the mesh deteriorates, and often computations cannot be continued without remeshing the interface. Situations often occur in which distances between some interface mesh points or some angles in the interface triangulation become very small. In earlier work, the present authors introduced a new methodology to approximate curvature driven curve and surface evolution, see [5,4,6]. The method has the important feature that interface mesh properties remain good during the evolution. In fact, for curves semidiscrete versions of the approach lead to polygonal approximations where the vertices are equally spaced throughout the evolution. The approach has been successfully used for various geometric evolution equations, including surface diffusion, [5], and grain boundary motion, [8], and was recently applied to crystal growth phenomena, see [9,10]. In this paper we generalize the approach to two-phase incompressible Stokes flow. The generalization to the Navier-Stokes case will be considered in the forthcoming article [12].
Studies of other groups reveal that the main source of spurious velocities in two-phase flow with surface tension is the fact that discontinuous functions allowing for jumps at the interface are not in the pressure space for an unfitted bulk mesh. In addition, it has been observed that the size of the spurious velocities depend both on the calculation of the interface curvature in the surface tension term and on the approximation properties of the pressure space, see e.g. [18]. We address the issue of obtaining a good approximation of the curvature term by using the tangential degrees of freedom in the interface representation appropriately. This approach is based on our earlier work on geometric evolution equations, where we observed good approximation properties for the curvature also for surfaces in three spatial dimensions, see [6,7]. The issue concerning the approximation properties of the pressure space is taken into account by an adaptive refinement of the bulk mesh close to the interface. In addition, we extend the pressure space by one degree of freedom by the addition of the characteristic function of one of the phases. It turns out that this enriched space in a semidiscrete version leads to exact volume conservation (area conservation in 2d) for the two phase regions. This new pressure space also ensures that stationary spherical states can be computed exactly. In particular, it turns out that the approach eliminates spurious velocities in the standard test case of a spherical drop in equilibrium, see e.g. [18] for the difficulties other approaches have with this simple test case. Moreover, in more general situations spurious velocities either do not appear at all or are small. We also observe that the conditioning of the resulting linear systems is not so badly effected, in contrast to other XFEM approaches 2 which involve far more degrees of freedom.
Let us state a few properties of our scheme.
• The semidiscrete continuous-in-time version of our scheme is stable in the sense that for a closed system the total interface energy decreases in time at a rate given by the energy dissipated. Similarly, for the fully discrete scheme the rate of decrease is at least that given by the energy dissipated.
• If no outer forces act, then we can easily show that any discrete stationary solution must have zero velocity, i.e. we can prove that there are no stationary solutions with spurious velocities. In addition, we can prove the existence of such stationary solutions for our scheme.
• For the semidiscrete version of our scheme we obtain in two spatial dimensions that the interface mesh points are equally spaced. In three space dimensions the semidiscrete solutions are conformal polyhedral surfaces, see [6], which are known to have good mesh properties. In practice we also observe for the fully discrete scheme that the computed discrete interfaces have good mesh properties.
In particular, no remeshings of the discrete interface are necessary.
• In two spatial dimensions polygonal curves with equidistributed vertices on a circle have constant discrete mean curvature and lead to discrete solutions with zero velocity and a constant pressure jump across the interface. In three spatial dimensions we numerically compute stationary polyhedral approximations of a sphere with constant discrete mean curvature. These polyhedral surfaces lead to zero velocity solutions with a constant pressure jump. These solutions are the natural discrete analogues of the stationary circle/sphere, which are the only stationary solutions of this incompressible two-phase flow in the case when no outer forces act.
• Our scheme can be applied both in the fitted interface mesh approach, as well as in the unfitted approach. The latter avoids the repeated remeshing of the bulk mesh to ensure that it remains fitted to the interface mesh, while the former naturally captures pressure jumps at the interface with standard pressure finite element spaces.
The remainder of the paper is organized as follows. In Section 2 we describe the mathematical model of two-phase flow and introduce a suitable weak formulation. In Section 3 we propose our discretization and establish existence, stability and other theoretical results. Finally, in Section 4 we present some numerical results in two and three space dimensions, including some convergence experiments for stationary and expanding bubble problems.

Mathematical setting
We consider the governing equations for the motion of unsteady, viscous, incompressible, immiscible two fluid systems. For low Reynolds numbers one can neglect the inertia terms, and so the equations for the velocity u and the pressure p are given by where Ω + (t) and Ω − (t) are the time dependent regions occupied by the two fluid phases. For a given domain Ω ⊂ R d , where d = 2 or d = 3, we now seek a time dependent interface (Γ(t)) t∈[0,T ] such that Γ(t) is completely contained in Ω and separates it into the two domains Ω + (t) and Ω − (t). Here the phases could represent two different liquids, or a liquid and a gas. Common examples are oil/water or water/air interfaces, see Figure 2 for an illustration. For later use, we assume that (Γ(t)) t∈[0,T ] is a sufficiently smooth evolving hypersurface parameterized by x(·, t) : Υ → R d , where Υ ⊂ R d is a given reference manifold, i.e. Γ(t) = x(Υ, t). Then V := x t . ν is the normal velocity of the evolving hypersurface Γ, where ν is the unit normal on Γ(t) pointing into Ω + (t). We consider cases in which the viscosity of the two fluids can be different and introduce µ(t) = µ + X Ω + (t) + µ − X Ω − (t) , with µ ± ∈ R >0 denoting the fluid viscosities; where, here and throughout, X A denotes the characteristic function for a set A.
In order to define the conditions that have to hold on the interface Γ(t), we introduce the stress tensor where Id ∈ R d×d denotes the identity matrix and D( u) := 1 2 (∇ u + (∇ u) T ) is the rate-of-deformation tensor. Now on the free surface Γ(t), the following conditions have to hold: where γ > 0 is a positive constant and κ denotes the mean curvature of Γ(t), i.e. the sum of the principal curvatures of Γ(t). Here we have adopted the sign convention that κ is negative where Ω − (t) is locally convex. Moreover, as usual, [ u] + − := u + − u − and [σ ν] + − := σ + ν − σ − ν denote the jumps in the velocity and the normal component of stress across the interface Γ(t). Here and throughout we employ the shorthand notation g ± := g | Ω ± (t) for a function g : Ω × [0, T ] → R d ; and similarly for scalar and matrix-valued functions. To close the system we prescribe the initial data Γ(0) = Γ 0 and the boundary condition u = 0 on ∂Ω. Using the fact that the velocity is divergence free, we can rewrite the total system as which is appropriate for the weak formulation considered in this paper.
For the system (2a-g) an a priori energy bound holds, and our goal is to introduce a discretization that satisfies a discrete analogue. First, on noting (1) and (2e), we have that where, here and throughout, L d and H d−1 denote the Lebesgue measure in R d and the (d − 1)-dimensional Hausdorff measure, respectively. Then using the identity Due to the incompressibility condition, the volume of the two fluids is preserved. We have from (2b,f) that It is a further aim that our discretization also preserves the fluid volumes.
In order to compute a discrete version of the mean curvature for polyhedral surfaces, we note the identity where ∆ s = ∇ s . ∇ s is the Laplace-Beltrami operator on Γ(t), with ∇ s . and ∇ s denoting surface divergence and surface gradient on Γ(t), respectively. We note that the sign convention in (7) is such that κ < 0 where Ω − (t) is locally convex. It is possible to use a weak formulation of this identity, which was first suggested by Dziuk, [15], and was used in [3] for the Navier-Stokes equations with a free capillary surface. A variant of Dziuk's approach, which leads to good mesh properties, has been introduced in [5] for d = 2 and in [6] for d = 3, and this will be the basis of our novel weak formulation. The main novelty, which inherently leads to good mesh properties, is that in [5,6] the mean curvature is treated as a scalar and is discretized separately from the normal ν, in contrast to discretizing κ := κ ν directly as in [15,3,18], see also §3.6 below. Before introducing our finite element approximation of (2ag), we will state an appropriate weak formulation. With this in mind, we introduce the function spaces where we recall that Υ is a given reference manifold. Let (·, ·) and ·, · Γ(t) denote the L 2 -inner products on Ω and Γ(t), respectively. On recalling (3), a possible weak formulation of (2a-g) is then given as follows. Find time dependent functions u, p, x and κ such that holds for almost all times t ∈ (0, T ], as well as the initial condition (2g). Here we have observed that if p ∈ P is part of a solution to (2a-g), then so is p + c for an arbitrary c ∈ R. Moreover, we note that for convenience we have adopted a slight abuse of notation in (8a-d). In particular, in this paper we will identify functions defined on the reference manifold Υ with functions defined on , and we denote both functions simply as v. For example, x ≡ id is also the identity function on Γ(t).

Discretization
We consider the partitioning Let us now introduce the space discretization.

Finite element spaces
For simplicity we consider Ω to be a polyhedral domain. Then for all m ≥ 0, let T m be a regular partitioning of We also introduce S m 0 , the space of piecewise constant functions on T m . Then our approximation to the velocity and pressure on T m will be finite element spaces U m ⊂ U and P m ⊂ P. We require also the space P m := P m ∩ P. The velocity/pressure finite element spaces where · 1 := · 0 + ∇ · 0 denotes the H 1 -norm on Ω; see e.g. [21, p. 114]. Here we take the reduced pressure space P m in (9) because where n denotes the outer unit normal to Ω. For example, we may choose the lowest order Taylor-Hood element P2-P1, the P2-P0 element or the P2-(P1+P0) element on setting respectively. It is well-known that these choices satisfy the LBB condition (9). We only remark that results for P2-P1 and P2-(P1+P0) need the weak constraint that all simplices have a vertex in Ω, see [13] for the P2-(P1+P0) element.
The parametric finite element spaces in order to approximate x and κ in (8a-d) are defined as follows. Similarly to [6], we introduce the following discrete spaces, based on the work of Dziuk, [15]. Let Γ m ⊂ R d be a (d − 1)-dimensional polyhedral surface, i.e. a union of non-degenerate (d − 1)-simplices with no hanging vertices (see [14, p. 164 denoting the standard basis of W(Γ m ). For later purposes, we also introduce . Throughout this paper, we will parameterize the new closed surface Γ m+1 over Γ m , with the help of a parameterization X m+1 ∈ V(Γ m ), i.e. Γ m+1 = X m+1 (Γ m ). Moreover, for m ≥ 0, we will often identify X m with id ∈ V(Γ m ), the identity function on Γ m .
For scalar and vector functions v, w on Γ m we introduce the L 2 -inner product ·, · Γ m over the current polyhedral surface Γ m as follows If v, w are piecewise continuous, with possible jumps across the Given Γ m , we let Ω m + denote the exterior of Γ m and let Ω m We then partition the elements of the bulk mesh T m into interior, exterior and interfacial elements as follows. Let Clearly is a disjoint partition, which in practice can easily be found e.g. with the Algorithm 4.1 in [11].
Here we assume that Γ m has no self intersections, and for the numerical experiments in this paper this was always the case. In addition, we define the piecewise constant unit normal ν m to Γ m such that ν m points into Ω m + . Of course, in the case of a fitted bulk mesh it holds that T m Γ m = ∅.

Finite element approximation
Following a similar approach used by the authors in [11] for crystal growth phenomena, we consider an unfitted finite element approximation of (8a-d). On recalling (11), we introduce the discrete viscosity µ m ∈ S m 0 , for m ≥ 0, as Clearly, for a fitted bulk mesh T m , (12) reduces to µ m = µ + X Ω m + + µ − X Ω m − ∈ S m 0 . Our finite element approximation is then given as follows. Let Γ 0 , an approximation to Γ(0), be given.
Note that here ∇ s denotes the surface gradient on Γ m , and so it depends on m. We observe that (13a-d) is a linear scheme in that it leads to a linear system of equations for the unknowns ( U m+1 , P m+1 , X m+1 , κ m+1 ) at each time level.
Furthermore, we note that in line with our earlier work in e.g. [5,6] we employ mass lumping on the first terms in (13c,d). The latter is necessary in order to be able to prove good mesh properties, while the former is then enforced to yield a stable scheme. On the other hand, it will become clear in §3.4 below that for conservation of the phase volumes we must use true integration for the remaining term in (13c), which in turn enforces true integration for the third term in (13a) for stability reasons.

Existence and stability results
Before we can prove existence and uniqueness results for a solution to our approximation to (8a-d), we have to introduce the notion of a vertex normal on Γ m , which is given as a weighted sum of the neighbouring normals. We will combine this definition with a mild assumption that is needed in the existence proof.
Then we further assume that dim span{ ω m k } We stress that (A) is a very mild assumption that is only violated in very rare occasions. For example, it always holds for surfaces without self-intersections; see [5,6] for more details. Given the above definitions, we introduce the piecewise linear vertex normal function and, on recalling (10) The next theorem establishes the existence of a unique solution to (13a-d) under some natural conditions. For later developments we also consider the reduced system (13a,c,d) with U m replaced by Clearly, any solution ( In the absence of an LBB condition for Proof. As the system (13a-d) is linear, existence follows from uniqueness. In order to establish the latter, we consider the system: It immediately follows from (17) and Korn's inequality that U = 0. In addition, it holds that X = X c ∈ R d . Together with (16c) for U = 0, (15a) and the assumption (A) this immediately yields that X = 0, while (16d) with η = π m [κ ω m ], recall (15a), implies that κ = 0. Finally, it now follows from (16a) with U = 0 and κ = 0, on recalling (9), that P = 0. Hence there exists a unique solution ( U m+1 , P m+1 , The proof for the reduced system (13a,c,d), with U m replaced by U m 0 is analogous. We now demonstrate that the discretization (13a-d) satisfies an energy estimate, which corresponds to the computation (5) in the continuous case. In particular, we obtain unconditional stability for our scheme.
Proof. Choosing ξ = U m+1 ∈ U m 0 in (13a), χ = γ κ m+1 in (13c) and η = γ ( X m+1 − X m ) in (13d) yields that Hence (18) follows immediately, where we have used the result  , U h (t) and P h (t) similarly to S m k , U m and P m , with the corresponding interpolation operators I h k and discrete approximations µ h (t) ∈ S h 0 (t), which will depend on Γ h (t) via the analogue of (12). Then, where f h := I h 2 f (t). In (20a-d) we always integrate over the current surface Γ h (t), with normal ν h (t), described by the identity function X h (t) ∈ V(Γ h (t)). Moreover, ·, · h Γ h (t) is the same as ·, · h Γ m with Γ m and X m replaced by Γ h (t) and X h (t), respectively; and similarly for ·, · Γ h (t) .
Using the results from [6] it is straightforward to show that , which is the discrete analogue of (4) on noting (20d). It is then not difficult to derive the following energy bound for the solution ( U h , P h , X h , κ h ) of the semidiscrete scheme (20a-d): Clearly, (21) is the natural discrete analogue of (5). In addition, it is possible to prove that the vertices of Γ h (t) are well distributed. As this follows already from the equations (20d), we refer to our earlier work in [5,6] for further details. In particular, we observe that in the case d = 2, i.e. for the planar twophase problem, an equidistribution property for the vertices of Γ h (t) can be shown, while in the case d = 3 it can be shown that Γ h (t) is a conformal polyhedral surface; see also (27) below.

XFEM Γ for conservation of the phase volumes
In general, the fully discrete approximation (13a-d) will not conserve mass, which means in particular that the volume L d (Ω m − ), enclosed by Γ m will in general not be preserved. Clearly, given that volume conservation holds on the continuous level, recall (6), it would be desirable to preserve the volume of the two phases also on the discrete level.
For the semidiscrete approximation (20a-d) from Remark 3 we can show conservation of the two phase volumes in the case that Choosing χ = 1 in (20c) and ϕ = ( L d (Ω) ) ∈ P h (t) in (20b), we then obtain, on recalling (15b) and as U h ∈ U, that which is the discrete analogue of (6). Clearly, for fitted bulk meshes T h (t) with S h 0 (t) ⊂ P h (t) the condition (22) trivially holds. In the case of the unfitted approach, on the other hand, discrete pressure spaces P h (t) based on piecewise polynomials, such as S h 0 , S h 1 or S h 0 + S h 1 , will in general not satisfy the condition (22). However, the assumption (22) can now be satisfied with the help of the extended finite element method (XFEM), see e.g. [23, §7.9.2]. Here the pressure spaces P m need to be suitably extended, so that they satisfy the time-discrete analogue of (22), i.e.
X Ω m − ∈ P m , which means that then (13b) implies U m+1 , ν m Γ m = 0, which together with (13c) and (15b) then yields that Hence the obvious strategy to guarantee (24) in the context of unfitted bulk meshes is to add only a single new basis function to the basis of P m , namely X Ω m − . We remark that in practice (24) leads to excellent phase volume conservation properties for the fully discrete scheme (13a-d). Moreover, we note that the contributions to (13a,b) coming from L d (Ω) ∈ P m can be written in terms of integrals over Γ m , since (∇ . ξ, 1) = 0 and for all ξ ∈ U m . We will call this particular enrichment procedure the XFEM Γ approach. For example, P m may be given by one of with T m being independent of Γ m .
We note that the above XFEM approach is different to that in e.g. [23,2,34], where in the level set context a standard finite element pressure space is enriched with numerous discontinuous basis functions in the vicinity of the interface in order to improve its approximation of the pressure jump across the interface. Some of these basis functions may have support on a small fraction of a bulk element cut by the interface, and this can lead to ill-conditioning of the associated linear system. Of course, the support of the single additional basis function in the XFEM Γ approach is very non-local, and is given by Ω m − . However, its contributions are easily calculated with the help of surface integrals; recall (25).
Similarly to the other XFEM approaches, we are unable to prove that the elements (U m , P m ), where P m is given by one of (26), satisfy the LBB condition (9). This means that we cannot easily prove existence and uniqueness of the discrete pressure P m+1 ∈ P m for the system (13a-d) with P m given as in (26). However, the relevant result in Theorem 1, in the absence of an LBB condition, and Theorem 2 still hold.

Properties of discrete stationary solutions
We now consider stationary states, Γ m+1 = Γ m , of the fully discrete system. It follows from (15a) that a stationary solution to (13a-d) satisfies where we note (14). We recall from [5] that (27) in the case d = 2 implies that Γ m is equidistributed, with the possible exception of elements σ m j that are locally parallel to each other. Moreover, we recall from [6] that surfaces in R 3 that satisfy (27) are called conformal polyhedral surfaces.
Next we consider discrete stationary states when no outer forces act, i.e. when f = 0. Here it turns out that our stability results from §3.3 have an immediate consequence.
The above lemma guarantees, independently of the choice of µ ± , that no spurious velocities appear for discrete stationary solutions, Γ m+1 = Γ m . For the XFEM Γ approach we can even show that polyhedral surfaces with constant discrete mean curvature and zero velocity are stationary solutions.

Remark 6.
A stationary solution to the continuous problem with f = 0 is a circle (d = 2) or a sphere (d = 3) with zero velocity and a piecewise constant pressure with a discontinuity across the interface, see (36a,b) below. For d = 2, one can choose Γ m with equidistributed points on a circle as an approximation of this circle, i.e. a closed regular polygon. Such a Γ m has constant discrete curvature, i.e. there exists a κ ∈ R such that (28) is satisfied. Hence Lemma 5 yields that in this situation ( U m+1 , X m+1 , κ m+1 ) = ( 0, X m , κ) is the unique solution to the reduced system with f m+1 = 0.

Alternative curvature treatment
As mentioned in Section 2, there is an alternative way to approximate the curvature vector κ ν in (7). In contrast to the strategy employed in (13a-d), where κ and ν are discretized separately, it is also possible to discretize κ := κ ν directly, as proposed in the seminal paper [15]. We then obtain the following linear finite element approximation. For m = 0, . . . , M − 1, find U m+1 ∈ U m , P m+1 ∈ P m , X m+1 ∈ V(Γ m ) and κ m+1 ∈ V(Γ m ) such that and set Γ m+1 = X m+1 (Γ m ). A discretization based on (29a-d) has first been proposed by Bänsch in [3] for one-phase flow with a free capillary surface in the very special situation that Clearly, (30) requires the fitted approach and in that case can be satisfied e.g. for the lowest order Taylor-Hood element, P2-P1, and a piecewise quadratic variant of V(Γ m ), see [3], or for the MINI element, P1 bubble -P1, with the piecewise linear V(Γ m ). We note that if (30) holds, then (29a-d) can be equivalently rewritten as together with (29b,c). For a nonlinear variant of this scheme involving space-time finite elements, in the context of the Navier-Stokes equations with a free capillary surface, Bänsch proved existence, uniqueness and stability of discrete solutions, see [3]. It is not difficult to prove existence, uniqueness and stability results also for the linear scheme (29a-d) for two-phase Stokes flow. In particular, one can show that there exists a unique solution to (29a-d) that also satisfies the stability bounds (18) and (19). Similarly, the analogues of Lemma 4 and, if (30) is satisfied, of Lemma 5 hold. In the latter case we observe that the curvature part of the discrete stationary solution is given by the unique κ m+1 ∈ V(Γ m ) such that κ m+1 , η Γ m = κ ν m , η Γ m for all η ∈ V(Γ m ). We stress that if (30) does not hold, then it does not appear possible to show the analogue of Lemma 5, which means that it is not possible to prove the existence of discrete stationary solutions for (29a-d).
However, the crucial difference between (29a-d) and (13a-d) is that in (29c) the tangential velocity of the discrete interface is fixed by U m+1 , and this has two consequences. Firstly, there is no guarantee that the mesh quality of Γ m will be preserved. In fact, as mentioned in the Introduction, typically the mesh will deteriorate over time. And secondly, even for the case that X Ω m − ∈ P m , it is not possible to prove (24) for (29a-d), as χ = ν m is not a valid test function in (29c), and so true volume conservation in the semidiscrete setting, recall (23), cannot be shown. It is for these reasons that we prefer to use (13a-d).
We now return to the equivalent rewrite (31), (29b,c) of (29ad). This has the advantage that the explicit computation of the curvature vector κ m+1 is avoided, although in practice the gain in computational efficiency is negligible because the main computational task is to solve the Stokes equations in the bulk. We recall that for the equivalence between (31), (29b,c) and the original (29a-d) it was crucial to enforce the strong assumption (30), which relies on a fitted bulk mesh. A variant of this rewrite can be obtained for (29a-d), with ·, · Γ m replaced by ·, · h Γ m , also in the absence of the assumption (30). This then leads to (31) with ∇ s ξ replaced by ∇ s ( π m ξ). Being equivalent, this rewrite inherits all the theoretical properties of (29a-d), with ·, · Γ m replaced by ·, · h Γ m , namely the stability results (18) and (19) and the analogue of Lemma 4.
It is important to note that unconditional stability for (31), (29b,c), even in the presence of (30), can no longer be shown if the third term in (31) is changed to where we note that ∇ s X m , ∇ s ξ Γ m = 1, ∇ s . ξ Γ m . Proving existence and uniqueness for this simpler variant is trivial, but it is no longer possible to establish stability, or the analogues of Lemmas 4 and 5. Alternatively, the term in (32), in view of (29c), can also be replaced by which collapses to (31) if the condition (30) holds. The formulation (31), (29b,c) with the third term in (31) replaced by (33) has been exploited in [19] for an axisymmetric formulation of two-phase Navier-Stokes flow. Finally, we mention that the ideas presented in this section on how to discretize the curvature term arising from (2e) can also be applied in the level set approach. Here (29c) is replaced with an approximation of the level set transport equation and this is then combined with (29a,b), where the third term in (29a) is replaced by (32) with Γ m now being a suitable reconstruction of the discrete interface arising from the zero level set of a discretization of the level set function φ in (34). See e.g. [22,23,2] for some examples. We stress that the level set method is a convenient computational tool for the interface motion, but that it does not appear possible to derive stability results in the spirit of (18) and (19) for the level set method.

Solution methods
As is standard practice for the solution of linear systems arising from discretizations of Stokes and Navier-Stokes equations, we avoid the complications of the constrained pressure space P m in practice by considering an overdetermined linear system with P m instead. With a view towards some numerical test cases in Section 4, we also allow for an inhomogeneous Dirichlet boundary condition g on ∂Ω and for ease of exposition consider only piecewise quadratic velocity approximations. Then we reformulate (13a-d) as follows. Find If (U m , P m ) satisfy the LBB condition (9), then the existence and uniqueness proof for a solution to (13a,c,d), (35) is as before. In the absence of (9), the existence and uniqueness of a solution to the reduced system that is analogous to (13a,c,d), with U m replaced by U m 0 , hinges on the nonemptiness of the set U m 0 ( g) := { U ∈ U m ( g) : (∇ . U, ϕ) = 0 ∀ ϕ ∈ P m }. The  linear system (13a,c,d), (35) can be solved with the help of a Schur complement approach, which reduces the system to a standard saddle point problem arising from discretizations of Stokes problems.

Numerical results
For details on the assembly of the linear system arising at each time step, as well as details on the adaptive mesh refinement algorithm and the solution procedure, we refer to the forthcoming article [12]. In particular, we recall that the scheme in general uses an adaptive bulk mesh that has a fine mesh size h f around Γ m and a coarse mesh size h c further away from it. The special case h f = h c leads to a uniform bulk mesh which will be sufficient for some of the simple test problems considered in this section. For all the numerical results presented in this paper no refinement or remeshing procedure was applied to the discrete interface approximations Γ m . We remark that we implemented our scheme with the help of the finite element toolbox ALBERTA, see [35].
In order to test our finite element approximation (13a-d), we consider the trivial true solution of a stationary circle/sphere, as it has been considered in e.g. [18,22]. In particular, together with where λ(t) = λ(0) = γ (d − 1) [r(0)] −1 , is an exact solution to the problem (2a-g). A nontrivial divergence free and radially symmetric solution u can be constructed on a domain that does not contain the origin. To this end, consider e.g.
together with is an exact solution to the problem (2a-g) with the homogeneous right hand side in (2c) replaced by g, where g( z) = α | z| −d z.
From now on we fix Γ(0) = { z ∈ R d : | z| = 1 2 }. Throughout this section we use uniform time steps τ m = τ, m = 0, . . . , M−1. For later use, we define h m Γ := max j=1,...,J m Γ diam(σ m j ). We also define the errors In order to evaluate the errors in the pressure, we define When we use XFEM Γ , we also evaluate the following errors for the pressure, P c − p c L 2

Numerical results in 2d
For our first set of experiments we fix Ω = (−1, 1) 2 and use the true solution (36a,b) for the parameters This means that the true solution reduces to r(t) = 1 2 , u(·, t) = 0 and p(t) = 2 [X Ω − (0) − 1 4 L 2 (Ω − (0))] for all t ≥ 0. Some errors for our approximation (13a-d) for the P2-P1 element can be seen in Table 1. The same convergence test for the pressure spaces P0 and P1+P0 are shown in Tables 2 and 3, respectively. Here we always choose uniform spatial discretizations such that h c = h f = h and h m Γ ≈ h/8. We varied the time discretization parameter τ from 10 −2 to 10 −4 , but as the errors were nearly indistinguishable for these runs, we only report on the computations with τ = 10 −2 . It appears that the three errors Tables 1-3   2 ), respectively. The latter is to be expected for the approximation of the discontinuous pressure on uniform meshes.
We note that for the experiments in Tables 1-3 we choose an equidistributed approximation Γ 0 of the circle Γ(0). In this special case our approximation with XFEM Γ , i.e. (13a-d) with X Ω m − ∈ P m , yields the exact solution U m+1 = 0. In particular, on recalling Lemma 5, we have that the unique solution to the reduced system (13a,c,d), with U m replaced by U m 0 , for m = 0 is given by U m+1 = 0, X m+1 = X 0 and κ m+1 = −λ 0 ∈ R, where −λ 0 ≈ −2 approximates the curvature of Γ(0), and by induction for all m = 0, . . . , M − 1. This implies that U m+1 = 0, , X m+1 = X 0 and κ m+1 = −λ 0 is a solution to (13a-d), and it is this solution that is found by our solution method in practice, see Table 4. We remark that these results remain unchanged for nonconstant µ, e.g. when choosing µ + = 10 µ − = 1. We visualize the final pressures for the finest runs in Tables 1-4 in Figure 3. Here in the case of the enrichment XFEM Γ being used, we plot P M c , which is almost identically equal to a constant, and λ M X Ω M−1 − separately. We now demonstrate that this remarkable property is generic to our method, in the sense that the circular, equidistributed numerical steady state solution is recovered by our method even if we choose very noncircular or very nonuniform initial data Γ 0 . Of course, this is the discrete analogue of the fact that circles are the unique steady state solutions in the continuous case, recall Section 1. In particular, we choose Γ 0 to be a very nonuniform approximation of Γ(0), where we represent the upper half of the circle by a single vertex, while the lower half is properly resolved to resemble a semicircle. In total we use K 0 Γ = 64 vertices for Γ 0 , and we use an adaptive bulk mesh with h c = 8 h f = 2 − 1 2 . Choosing τ = 10 −4 we simulate the evolution with our scheme for the time interval [0,5]. In Figure 4 we show some snapshots of the evolution, while in Figure 5 a plot of U m L ∞ (Ω) over time can be seen. Here we use the P2-P1 element. As expected, the approximations Γ m converge towards an equidistributed circle, while U m converges to zero.
Here we always choose the spatial discretization parameters such that h c = 8 h f and h m Γ ≈ h f . The same convergence experiments now for the pressure spaces enriched with XFEM Γ can be found in Tables 14-16. We visualize the final pressures for the coarsest runs in Tables 11-16 in Figure 6. Here in the case of the enrichment XFEM Γ being used, we plot the pressure parts P M c and λ M X Ω M−1 − separately. We recall from (37b) that the jump λ(T ) in the pressure is made up from two components. A jump of 0.55 − 1 2 ≈ 1.35 due to the curvature, and a jump of 0.27/0.55 ≈ 0.49 due to the difference in µ. When the XFEM Γ enrichment of the pressure space is used, then the discretizations will treat these two jumps differently. As can be seen from (13a), the jump due to curvature can be absorbed by X Ω m − , while the jump in µ, recall (12), is left to the standard bulk pressure space. It appears from our plots in Figure 6 that this is indeed what happens in practice.
We end this section with a shear flow experiment that is motivated by similar simulations in [28]. Here we take Ω = (−1, 1) 3 and for u prescribe the inhomogeneous Dirichlet boundary condition g( z) = (z 3 , 0, 0) T on ∂Ω. We also let µ = 1 and γ = 3. For the discretization parameters we fix h c = 8 h f = 2 − 1 2 and τ = 10 −2 , and we use a triangulation of Γ 0 that is made up of 768 elements, and so h 0 Γ = 0.14. At first we present the numerical results for the alternative formulation (29a-d), see Figure 8, where we plot the discrete interface Γ m at times t = 0 and t = 4.2. We note the elongated elements and the high density of vertices in certain regions. In fact, the simulation breaks down soon after because of the coalescence of mesh points.
The same simulation for our approximation (13a-d) can be seen in Figure 9. Here the interface at the final time T = 5 is close to being a numerical steady state. The "sphericity" of Γ M , see [40], is given by π We observe an excellent mesh quality throughout the evolution, in contrast to the elongated elements that can be seen in Figure 8 and in e.g. [28,Fig. 14].

Conclusions
We have presented a novel front-tracking method for viscous incompressible two-phase flow which can be shown to be stable. The numerical method couples a parametric finite element approximation of the interface with a standard finite element approximation of the Stokes equations in the bulk. Here the bulk mesh may be chosen to be either fitted (or adapted) to the interface, or it can be totally independent of the interface mesh. In the latter case, we introduce an XFEM approach to guarantee that our scheme conserves the volumes of the two phases.
The stability of the proposed method implies that the wellknown static bubble test problem can be computed exactly. More generally, we can show that our scheme admits timeindependent discrete solutions, which all have the property that the velocity is zero. This means, in particular, that it is possible to eliminate spurious velocities for stationary solutions.
The second prominent feature of our numerical method is the excellent mesh quality of the interface approximation. This is induced by an inherent discrete tangential motion of the vertices that make up the discrete interface. In particular, for a semidiscrete continuous-in-time variant of our scheme it can be shown that the discrete interfaces are equidistributed polygonal curves (d = 2) and conformal polyhedral surfaces (d = 3), respectively.