A stable numerical method for the dynamics of fluidic membranes

We develop a finite element scheme to approximate the dynamics of two and three dimensional fluidic membranes in Navier–Stokes flow. Local inextensibility of the membrane is ensured by solving a tangential Navier–Stokes equation, taking surface viscosity effects of Boussinesq–Scriven type into account. In our approach the bulk and surface degrees of freedom are discretized independently, which leads to an unfitted finite element approximation of the underlying free boundary problem. Bending elastic forces resulting from an elastic membrane energy are discretized using an approximation introduced by Dziuk (Numer Math 111:55-80, 2008). The obtained numerical scheme can be shown to be stable and to have good mesh properties. Finally, the evolution of membrane shapes is studied numerically in different flow situations in two and three space dimensions. The numerical results demonstrate the robustness of the method, and it is observed that the conservation properties are fulfilled to a high precision.


Introduction
The evolution of lipid bilayer membranes is driven by the bending energy, which involves the curvature of the membrane, and hydrodynamics. Lipid membranes typically form vesicles, i.e. bag-like structures containing fluid, which are surrounded by a possibly different fluid. The omnipresence of membranes in biological systems has led to a growing interest in vesicles over the past decades. Much of the work on vesicles was motivated by the fact that their shape at rest resembles the biconcave forms of red blood cells. It is the goal of this paper to introduce, and analyze, a finite element method of a model for the evolution of lipid membranes, which was introduced by Arroyo and DeSimone [2].
Their model couples a tangential Navier-Stokes system on the membrane to a bulk Navier-Stokes system. On the membrane, forces stemming from a first variation of the curvature energy appear. In this paper, we introduce a finite element method discretizing the bulk and surface degrees of freedoms independently, which leads to an unfitted approximation of the two-phase flow problem. The forces resulting from the elastic membrane energy are discretized using an approach that was introduced by Dziuk [18] for Willmore flow. The two Navier-Stokes systems are coupled and subsequently discretized with the help of a suitable variational formulation, which allows us to show a stability bound for the discretization of the underlying complex free boundary problem. It will turn out that the surface finite element mesh has good mesh properties also for strongly deforming membrane evolutions. This fact results from a suitable discrete local incompressibility condition on the surface, which we will also analyse. Before we state the governing equations and the numerical method in more detail, we discuss the physical background and review approaches used by other authors to numerically solve similar problems.
The bending energy for a lipid membrane used in this paper is where the bilayer is modelled as a closed hypersurface Γ in R d , d = 2 or 3. By we denote the mean curvature (the sum of the principal curvatures) of Γ , α ∈ R >0 is the bending rigidity and dH d−1 indicates integration with respect to the (d−1)-dimensional surface measure. In the simplest energetical model for vesicles one minimizes the energy E α (Γ ) under the constraints that the area of Γ is fixed and that Γ encloses a fixed volume. The latter is due to the fact that the osmotic conditions of the fluids surrounding the membrane lead to a fixed volume. Furthermore, the vesicle can be considered as locally incompressible, which leads to a fixed total surface area. For a deeper physical discussion of these conditions we refer to the overview article [39], where also other aspects of fluidic membranes and vesicles are thoroughly discussed.
In the fluid regions, Ω − and Ω + , inside and outside of the membrane, one requires the incompressible Navier-Stokes equations, i.e.
At typical temperatures the membrane itself is in a fluidic state, which leads to the fact that on the membrane the incompressible surface Navier-Stokes equations have to hold. Here ∂ • t denotes the material time derivative on Γ , σ Γ is the surface stress tensor, ν is the unit normal on Γ and [σ ν] + − describes stresses acting on the membrane via the normal stresses σ ν from both sides of the membrane (see Sect. 2 for precise definitions). The operator ∇ s . is the surface divergence and ∇ s . u = 0 models the fact that the membrane is locally incompressible.
Interfacial fluid mechanics was first thoroughly discussed by Scriven [38], generalizing earlier ideas of Boussinesq. In this context the surface stress tensor σ Γ was first introduced, and is hence called the Boussinesq-Scriven tensor. In addition, α f Γ ν models forces acting on the membrane, which result from the curvature elasticity E α (Γ ). The forces act in a direction normal to the membrane and f Γ is given as minus the first variation of E(Γ ), i.e.
where Δ s is the surface Laplacian and ∇ s is the surface gradient.
In recent years, many papers have appeared which numerically approximate the where V is the normal velocity of the evolving membrane Γ . This geometric evolution equation is called Willmore flow, and we refer to [6,7,14,15,17,18,22,29] for different computational approaches to Willmore flow. Since the enclosed volume and the total surface area are preserved for lipid membranes, the volume and area preserving variant of (1.2), which is called Helfrich flow, is of particular interest. Helfrich flow has been considered numerically in e.g. [6,12]. Other authors included additional physical effects in the geometrical model, such as lateral inhomogeneity and line tension effects, see [20,31]. In [13] a fluid-membrane system, in which forces resulting from the Willmore energy act on an interior flow, is considered. In that model surface area is maintained with the help of a global Lagrange multiplier. As pointed out above, the membrane is locally incompressible and hence the condition ∇ s . u = 0 should be enforced on the flow. This condition has been dealt with in numerical simulations in [27,35,36] within a level set context, in [1,24] with the help of a phase field approach and in [23] by using an immersed boundary method. In these approaches the local incompressibility constraint on the membrane is enforced by a Lagrange multiplier leading to an inhomogeneous surface pressure. However, in the computations in the latter paper the constraint is relaxed by a spring-like elastic force. In addition, there exists work on the surface Stokes system without taking the bulk fluid flow into account. There the volume conservation is enforced by a global Lagrange multiplier. We refer to [33,34] for numerical results using this modelling variant. The only numerical work taking simultaneously surface and bulk viscosity effects in the fluidic membrane evolution into account is [3]. However, these results are restricted to the axisymmetric situation. In addition, their numerical method cannot be shown to be stable, as is the case for all of the above numerical methods.
We also refer to numerical work in [21,26,30,40] on the evolution of red blood cells, which study the influence of the elastic effects resulting from the cytoskeleton on the membrane evolution. Finally, we mention that analytical well-posedness issues for the model considered in this paper are currently being addressed in [25,28].
Building on earlier work by the present authors on two-phase flow and by Dziuk [18] on Willmore flow, it is the main goal of this paper to introduce and analyze a numerical method for the full membrane evolution problem. Our numerical approach has the following features: -The bulk and surface degrees of freedom are discretized with standard bulk and surface finite elements leading to an unfitted finite element method. -The effects of the bulk fluid and of the fluidic membrane are taken into account simultaneously. In particular, surface viscous effects are accounted for through the Boussinesq-Scriven law. -Local volume and local membrane area conservation result naturally from the volume and surface incompressibility conditions. Local area conservation can be shown for a continuous-in-time semidiscrete variant of our proposed scheme. In addition, for a simple modification of our scheme, which can be interpreted as a virtual element method, see e.g. [41], volume conservation properties can also be shown. -Elastic forcing from the curvature energy E(Γ ) is taken into account, and this is discretized with the help of a weak formulation due to [18]. -Stability of a semidiscrete version can be shown. To our knowledge, this is the first stability result in the literature for a numerical approximation of the dynamics of fluidic membranes. -The interface is advected with the help of the fluid velocity. In other fluid flow problems with a free boundary this typically leads to distortions of the parametric surface mesh, see the discussion in [4]. However, in our case the local surface area conservation ∇ s . u = 0 guarantees that the surface mesh quality remains good during the evolution, see Remark 4.1 and the numerical simulations in Sect. 7. We also refer to [32], where a strategy for the tangential redistribution of mesh points by conserving the relative surface area during the evolution was designed. -Fully three dimensional simulations, without making any symmetry assumptions, have been performed, and to the knowledge of the authors this paper presents the first such numerical computations for the full fluidic membrane problem, i.e. taking the bulk viscosity, the surface viscosity and the local incompressibility of the bulk and surface fluid into account.
The outline of the paper is as follows. After introducing the governing equations in Sect. 2, we present a weak formulation in Sect. 3. This weak formulation is the basis for our semidiscrete and fully discrete finite element approximations, which are formulated and analyzed in Sects. 4 and 5. After stating the solutions methods in Sect. 6, we present numerical simulations in Sect. 7.

Governing equations
In this section we state the equations governing the evolution of fluidic membranes, as introduced in [2 , occupied by the outer phase, and a domain Ω − (t) := Ω \Ω + (t), which is occupied by the inner phase, see Fig. 1 for an illustration. For later use, we assume that (Γ (t)) t∈[0,T ] is a sufficiently smooth evolving hypersurface without boundary that is parameterized by defines the velocity of Γ (t), and V := V . ν is the normal velocity of the evolving hypersurface Γ (t), where ν(t) is the unit normal on Γ (t) pointing into Ω + (t). Moreover, we define the space-time surface , with ρ ± ∈ R ≥0 , denote the fluid densities, where here and throughout X A defines the characteristic function for a set A . Denoting by u : Ω × [0, T ] → R d the fluid velocity, by σ : Ω × [0, T ] → R d×d the stress tensor, and by f : Ω × [0, T ] → R d a possible volume force, the incompressible Navier-Stokes equations in the two phases are given by denotes the boundary of Ω with outer unit normal n. Hence (2.2c) prescribes a possibly inhomogeneous Dirichlet condition for the velocity on ∂ 1 Ω, which collapses to the standard no-slip condition when g = 0, while (2.2d) prescribes a stress-free condition on ∂ 2 Ω. Throughout this paper we assume that H d−1 (∂ 1 Ω) > 0. We will also assume w.l.o.g. that g is extended so that g : Ω → R d . In addition, the stress tensor in (2.2a) is defined by where Id ∈ R d×d denotes the identity matrix and D( u) : , with μ ± ∈ R >0 , denotes the dynamic viscosities in the two phases. On the free surface Γ (t), the following conditions need to hold: where ρ Γ ∈ R ≥0 denotes the surface material density, α ∈ R >0 is the bending rigidity and f Γ := f Γ ν is defined by (1.1). In addition, ∇ s . denotes the surface divergence on Γ (t), and the surface stress tensor is given by where μ Γ ∈ R ≥0 is the interfacial shear viscosity and p Γ denotes the surface pressure, which acts as a Lagrange multiplier for the incompressibility condition (2.4c). Here is the projection onto the tangent space of Γ (t), and Here and throughout, we employ the shorthand notation b ± := b | Ω ± (t) for a function b : Ω × [0, T ] → R d ; and similarly for scalar and matrix-valued functions. In addition, denotes the material time derivative of ζ on Γ (t). We compute ∂ • t ζ with the help of an extension of ζ to a neighborhood of G T . Here we stress that the derivative in (2.7) is well-defined, and depends only on the values of ζ on G T , even though ζ t and ∇ ζ do not make sense separately for a function defined on G T ; see e.g. [19, p. 324 where Γ 0 ⊂ Ω and u 0 : Ω → R d are given initial data satisfying ρ ∇ . u 0 = 0 in Ω, ρ Γ ∇ s . u 0 = 0 on Γ 0 and ρ + u 0 = ρ + g on ∂ 1 Ω. Of course, in the case ρ − = ρ + = ρ Γ = 0 the initial data u 0 is not needed. Similarly, in the case ρ − = ρ + = 0 and ρ Γ > 0 the initial data u 0 is only needed on Γ 0 . However, for ease of exposition, and in view of the unfitted nature of our numerical method, we will always assume that u 0 , if required, is given on all of Ω. It is not difficult to show that the conditions (2.2b) enforce volume preservation for the phases, while (2.4c) leads to the conservation of the total surface area H d−1 (Γ (t)), see Sect. 3 below for the relevant proofs. As an immediate consequence we obtain that spheres remain spheres, and that spheres with a zero bulk velocity are stationary solutions. In addition, in the case d = 2 the condition (2.4c) immediately implies that D s ( u) = 0 on Γ (t), which means that in two space dimensions the problem is independent of the value of μ Γ .
Furthermore, we note that Here denotes the mean curvature of Γ (t), i.e. the sum of the principal curvatures where we have adopted the sign convention that is negative where Ω − (t) is locally convex. In particular, it holds that Γ . Thereforep ± is constant on Ω ± (t). In addition, since ∇ spΓ is tangential, we obtain that ∇ spΓ = 0, and hencep Γ is a constant. Moreover, (2.11b) implies that p Γ =p + −p − . So if is not constant, which is the case if Γ (t) is not a sphere, thenp Γ = 0 andp + =p − . Hence p Γ in this case is unique, and p is unique in Ω up to an additive constant. If is constant, however, i.e. if Γ (t) is a sphere, then p Γ is only unique up to an additive constant, which is fixed by the two additive constants in the bulk phases. For more details see the discussion around (3.13) below.
Finally, we recall that the source term where ·, · Γ (t) denotes the L 2 -inner product on Γ (t). It does not appear possible to derive a stable discretization of the system (2.2a-d), (2.3), (2.4a-d), (2.5) based on the formulation (1.1). Hence in this paper we will make use of the stable approximation of Willmore flow introduced in [18], which is based on a discretization of the curvature vector = ν, and on the identity , (2.12) where we note that our notation is such that ∇ s χ = (∇ Γ χ) T , with ∇ Γ χ = ∂ s i χ j d i, j=1 defined as in [18].

Weak formulation
Before introducing our finite element approximation, we will state an appropriate weak formulation. With this in mind, we introduce the following function spaces for a given b ∈ [H 1 (Ω)] d : In addition, we let P := L 2 (Ω) and define Letting (·, ·) and ·, · ∂ 2 Ω denote the L 2 -inner products on Ω and ∂ 2 Ω, respectively, we recall from [11] that it follows from (2.2a-d), (2.4a, d) and (2.3) that and where we have also noted for symmetric matrices A ∈ R d×d that A : B = A : Only slip or free-slip conditions were considered in [11], and so the boundary integral over ∂ 2 Ω did not appear there. But it is easily established that the more general (3.1) also holds, on noting [11, (3.2)]. We define, similarly to (2.7), the following time derivative that follows the parameterization x(·, t) of Γ (t), rather than u. In particular, we let where we stress once again that this definition is well-defined, even though ζ t and ∇ ζ do not make sense separately for a function ζ ∈ H 1 (G T ). On recalling (2.7) we obtain that Moreover, for later use we note that (3.5) see Definition 2.11 and Lemma 5.2 in [19], respectively.
In what follows we would like to derive an energy bound for a solution of (3.7a-f), where for ease of exposition we consider only the case g = 0. All of the following considerations are formal, in the sense that we make the appropriate assumptions about the existence, boundedness and regularity of a solution to (3.7a-f). Firstly, it follows from (3.5), (3.7d) and (3.7c . (3.9) (3.10) Moreover, we note that it immediately follows from (3.5) and (3.7c, d) that In addition, the volume of Ω − (t) is preserved in time, i.e. the mass of each phase is conserved. To see this, choose Recalling the argument on the uniqueness of the pressures p and p Γ below (2.11a, b), we note the following LBB-type condition: 13) which we also refer to as the LBB Γ condition. Here we have defined the space In the case that the smooth hypersurface Γ (t) is not a sphere, then (3.13) is shown to hold if ∂ 1 Ω = ∂Ω is a smooth boundary in [28, p. 15].

Semidiscrete finite element approximation
For simplicity we consider Ω to be a polyhedral domain. Then let T h be a regular where P k (o) denotes the space of polynomials of degree k on o. We also introduce denotes the coordinates of the degrees of freedom of S h k (Ω), k ≥ 1. In addition we define the standard projection operator Our approximation to the velocity and pressure on T h will be finite element spaces For the finite element spaces (U h ( g), P h ) we may choose, for example, the lowest order Taylor-Hood element P2-P1, the P2-P0 element or the P2-(P1+P0) element on setting , respectively. The lowest order Taylor-Hood element satisfies the standard LBB condition in the bulk for d = 2 and d = 3, while the other two choices satisfy it for d = 2.
For the numerical approximation of the evolution of fluidic membranes it is desirable to maintain the surface area of the interface, recall (3.11), as well as the volume of the two phases, recall (3.12). In [8,11] the present authors augmented the pressure space by the characteristic function of the inner phase in order to obtain discretizations that maintain the volume of the two phases. This enrichment of the pressure space is an example of an XFEM approach, and we refer to this particular approach as XFEM Γ . Unfortunately, it does not appear possible to prove a discrete analogue of (3.11) for the XFEM Γ approach from [8,11]. Hence in this paper we will modify the XFEM Γ approach so that we obtain numerical approximations that satisfy discrete analogues of both (3.12) and (3.11). From a practical point of view, this approach is very close to the procedure in [8,11]. But the introduced modifications mean that the adjustments to the finite element approximations no longer have an interpretation within the XFEM framework. However, the new approach may be interpreted as an example of the recently proposed virtual element method, see below for further details.
The parametric finite element spaces in order to approximate e.g. and p Γ are defined as follows, see also [5] e. a union of non-degenerate (d − 1)-simplices with no hanging vertices (see [16, p. 164 We also introduce S h 0 (Γ h (t)), the space of piecewise constant functions on {σ h j (t)} J Γ j=1 . For ease of presentation, we introduce the following notations for the spaces of piecewise linear functions on For later purposes, we also introduce the standard interpolation operators π h k (t) : If v, w are piecewise continuous, with possible jumps across the edges of Following [19, (5.23)], we define the discrete material velocity for z ∈ Γ h (t) by Then we define, similarly to (3.3), For later use, we also introduce the finite element spaces as well as On differentiating (4.1) with respect to t, it immediately follows that see [19,Lem. 5.5]. It follows directly from (4.4) that Similarly, we recall from [9, Lem. 3.1] that We then partition the elements of the bulk mesh T h into interior, exterior and interfacial elements as follows. Let Clearly . Moreover, we introduce the discrete density ρ h (t) ∈ S h 0 (Ω) and the discrete viscosity We introduce, similarly to (2.6a, b), In what follows we will introduce a finite element approximation for the free boundary problem . When designing such a finite element approximation, a careful decision has to be made about the discrete tangential velocity of Γ h (t). The most natural choice is to select the velocity of the fluid, i.e. (3.7d) is appropriately discretized, and that is the approach we adopt in this paper. Overall, we then obtain the following semidiscrete continuous-in-time finite element approximation, which is the semidiscrete analogue of the weak formulation (3.7a-f). Given where we recall (4.2). Here we have defined f h (·, t) := I h 2 f (·, t), where here and throughout we assume that f ∈ L 2 (0, T ; [C(Ω)] d ). We observe that (4.9d) collapses , which on recalling (4.3) turns out to be crucial for the stability analysis for (4.9a-f). It is for this reason that we use mass lumping in (4.9d). The superscript · (h) in (4.9e, f) means that we can consider the corresponding terms either with or without mass lumping. Here we note that the scheme (4.9d-f), with true integration used throughout, and with U h in (4.9d) replaced by F h Γ , is the stable approximation of Willmore flow from [18], see also [15] for the case d = 2. In fact, for d = 2 we observe that where · s in the last term denotes differentiation with respect to arclength; compare also [7, (3.12 a, b)].
In the following theorem we derive discrete analogues of (3.10) and (3.11) for the scheme (4.9a-f).
In the case = 1, it immediately follows from (4.6) and (4.4), similarly to (3.11), on choosing η = χ h k in (4.9c), and on recalling from (4.9d) that While this property may appear desirable at first, our numerical experience for the fully discrete variant of this modified (4.9a-f) indicates that in the case d = 3 the constraint (4.17) is too severe. For example, the evolution for the fully discrete analogues of Γ h (t) may lag behind the observed evolution for the equivalent simulation with S h 0 (Γ h (t)) replaced by W (Γ h (t)). This can even lead to locking, where the linear solvers are no longer able to find a discrete solution at a certain time step. Here we note that for typical triangulations of Γ h (t) it holds that J Γ ≈ 2 K Γ . It is for this reason that we prefer the scheme (4.9a-f) as stated.

Remark 4.2
For the proofs of (4.11) and (4.12) it is crucial that (4.9c) holds with linear interpolation, i.e. = 1. Similarly, for the proof of (4.10) with ρ Γ > 0 it is necessary to choose = 1. For the surface Navier-Stokes system, this means that we use linear velocity approximations with linear pressures, something that is unlikely to satisfy a discrete surface LBB condition. In fact, in practice this can lead to an oscillatory approximation of the surface pressure, as can be seen in Fig. 5 below.
This can be avoided by using a quadratic interpolation of the bulk velocities, i.e. choosing = 2. This gives better behaved surface pressure approximations in practice, but the mesh quality is no longer maintained. Moreover, it can no longer be shown that (4.12) holds. However, in practice the evolutions of the fully discrete analogues of Γ h (t) are nearly identical for the two cases = 1 and = 2, at least when d = 2. When = 2 that means that in practice the surface area is maintained well, while for = 1 it means that despite the oscillatory surface pressures, the approximation of the velocity is well-behaved. Note that in three dimensional simulations we observe that the scheme with = 2 does not conserve the overall surface area well. Hence it does not appear to be a viable scheme in practice.
We observe that it does not appear possible to prove a discrete analogue of (3.12) for the scheme (4.9a-f). The reason is that χ = ν h is not a valid test function in (4.9d). However, a procedure similar to the XFEM Γ approach introduced by the authors in [8,11] ensures that a modified variant of (4.9a-f) conserves the enclosed volumes. To this end, we introduce the vertex normal function For later use we note that We are now in a position to propose the following adaptation of (4.9a-f).
and (4.9c-f) hold. Of course, χ = ω h is a valid test function in (4.9d), and so combining with (4.19b) yields a discrete volume preservation property, as is shown in the following theorem.
In order to interpret the adaptation in (4.19a, b) physically, we note the following. Of course, P h and P h sing act as Lagrange multipliers for the conditions in (4.19b). Moreover, on noting that we observe that replacing corresponds to augmenting P h in (4.9a, b) with the single additional basis function . This is the XFEM Γ approach introduced by the authors in [8,11], which for the schemes introduced there naturally leads to the conservation of the volume of the two phases. Such an XFEM interpretation is no longer possible for the modifications (4.21), as one cannot identify the corresponding additional basis function in the bulk. Therefore this can be viewed as an example of the recently proposed framework of virtual element methods, see e.g. [41]. In addition, on recalling (4.18) we have for a fixed time t ∈ [0, T ] that . It follows from (4.22) that we can interpret P h (·, t) + P h as the natural approximation to the pressure p(·, t) arising from (4.19a, b), (4.9c-f).

Fully discrete finite element approximation
In this section we consider a fully discrete variant of the scheme (4.19a, b), (4.9c-f) from §4. Here we will choose the time discretization such that existence and uniqueness of the discrete solutions can be guaranteed, and such that we inherit as much of the structure of the stable schemes in [8,11] We also introduce the standard interpolation operators 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 ).
Given Γ m , we let Ω m + denote the exterior of Γ m and let Ω m − denote the interior of In addition, we define the piecewise constant unit normal ν m to Γ m such that ν m points into Ω m + . We then partition the elements of the bulk mesh T m into interior, exterior and interfacial elements as before, and we introduce ρ m , μ m ∈ S m 0 (Ω), for m ≥ 0, as We also introduce the L 2 -inner product ·, · Γ m over the current polyhedral surface Γ m , as well as the the mass lumped inner product ·, · h Γ m . We introduce, similarly to (4.8a, b), We introduce the following pushforward operator for the discrete interfaces Γ m and Γ m−1 , for m = 0, . . . , M. Here we set Γ −1 : for m = 1, . . . , M, and set Π 0 −1 := π 0 1 . Analogously to (5.1) we also introduce Π m m−1 : C(Γ m−1 ) → W (Γ m ). We note, similarly to (4.18), that where ω m := K Γ k=1 χ m k ω m k ∈ V (Γ m ), and where for k = 1, . . . , K Γ we let Θ m k := { j : q m k ∈ σ m j } and set Λ m k := ∪ j∈Θ m k σ m j and ω m k := For the approximation to the velocity and pressure on T m we use the finite element spaces U m ( g) and P m , which are the direct time discrete analogues of U h ( g) and P h (t m ), as well as P m ⊂ P. Analogously to (3.13), we also say that (U m ( 0), P m , W (Γ m ), π m ) satisfy the discrete LBB Γ inf-sup condition if there exists a C 0 ∈ R >0 , independent of T m and {σ m j } J Γ j=1 , such that Unfortunately, it does not appear possible to prove that (5.2) holds for e.g. (U m ( 0), , because T m and Γ m are totally independent. Recall that also in the much simpler situation of the XFEM Γ approach from [8,11], which corresponds to setting η = 0 in (5.2) and replacing ω m , ξ h Γ m with ν m , ξ Γ m , the authors were unable to show that an LBB condition holds. In fact, in this simpler situation it is easily possible to construct a counterexample, e.g. when P m = S m 0 (Ω). Then, if Ω m − is a union of bulk elements, i.e. Γ m happens to be a union of bulk faces, clearly (5.2) does not hold, as X Ω m − ∈ P m . Hence for the choice Of course, in practice this situation never occurs, because the totally independent Γ m is never exactly aligned with the bulk mesh. Our proposed fully discrete equivalent of (4.19a, b), (4.9c-f), for a fixed ∈ {1, 2}, is then given as follows. Let Γ 0 , an approximation to Γ (0), as well as κ 0 ∈ V (Γ 0 ) and U 0 ∈ U 0 ( g) be given. For m = 0, . . . , M − 1, find U m+1 ∈ U m ( g), P m+1 ∈ P m , P m+1 sing ∈ R, P m+1 Here we have defined f m+1 := I m 2 f (·, t m+1 ). We observe that (5.3a-f) is a linear scheme in that it leads to a linear system of equations for the unknowns ( U m+1 , P m+1 , P m+1 sing , P m+1 at each time level. In the absence of the LBB Γ condition (5.2) we need to consider the reduced system In order to prove the existence of a unique solution to (5.3a-f) we make the following very mild well-posedness assumption.

Remark 5.1
We always choose U m ( 0) and P m so that the standard LBB condition in the bulk holds. Therefore U m 0 ( g), recall (5.4), is non-empty in the absence of the two Γ m constraints in (5.4). Clearly, there is the possibility of U m 0 ( g) being empty if the number of vertices on Γ m , K Γ , is too large compared to the number of bulk mesh vertices in the vicinity of Γ m . In practice we refine our bulk meshes in the neighbourhood of the interface, which lessens the likelihood of this occurring. In fact, in practice we do not experience problems for our choices of bulk and surface meshes. For example, in the case d = 2, we recall that the P2-(P1+P0) element satisfies the standard LBB condition in the bulk. This means that for the lowest order Taylor-Hood element P2-P1, which we employ in practice, there are some additional degrees of freedom in the velocity space, which prevent U m 0 ( g) from being empty in practice.

Remark 5.2
The scheme (5.3a-f) clearly leads to a coupled system of linear equations. On replacing F m+1 Γ in (5.3a) with F m Γ the system decouples into (5.3a-c) for ( U m+1 , P m+1 , P m+1 sing , P m+1 Γ ) and into (5.3d-f) for ( X m+1 , κ m+1 , F m+1 Γ ). Of course, the subsystem (5.3d-f) itself decouples into three equations for the three unknowns. While the decoupled system offers the advantage of being easier to solve, we found in practice that the coupled scheme (5.3a-f) preserved the surface area better than the decoupled scheme. An additional drawback of the decoupled scheme is that it is less stable and so in general needs smaller time steps than the coupled scheme (5.3a-f). The latter fact can partly be explained with the following observation.
On replacing F m+1 Γ with κ m+1 , and in the case ∂ 1 Ω = ∂Ω, g = 0 and ρ Γ = 0, we obtain an unconditionally stable approximation for two-phase flow in the spirit of [11], but with the additional side constraint (2.4c). In particular, for fixed bulk meshes in time one can show that see [11,Theorem 4.1] and [10, Theorem 4.2] for more details. It is for these reasons that we prefer the coupled scheme (5.3a-f).

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. Introducing the obvious abuse of notation, the linear system (5.3a-c) for α = 0, with P m replaced by P m , can be formulated as: Find ( U m+1 , P m+1 , P m+1 sing , P m+1 where ( U m+1 , P m+1 , P m+1 sing , P m+1 here denote the coefficients of these finite element functions with respect to the their standard bases. The definitions of the matrices and vectors in (6.1) directly follow from (5.3a-c), but we state them here for completeness in the case g = 0. Let i, j = 1, . . . , K m U , n, q = 1, . . . , K m P and k, l = 1, . . . , K Γ . Then where { e r } d r =1 denotes the standard basis in R d , and where we have used the convention that the subscripts in the matrix notations refer to the test and trial domains, respectively. A single subscript is used where the two domains are the same. The entries of D Ω , for i = 1, . . . , K m U , are given by The only new term compared to previous works by the authors on two-phase flows, see [10,11], is S Γ,Ω . Here we note that In order to provide a matrix-vector formulation for the full system (5.3a-f), and in particular in view of (5.3f), we recall from [18, p. 64 Hence, in addition to (6.2), we introduce the following matrices, where q = 1, . . . , K m U , and k, l = 1, . . . , K Γ Here we have made use of the fact that Moreover, it clearly holds that Denoting the system matrix in (6.1) as B Ω C C T 0 , and letting P m+1 = (P m+1 , P m+1 sing , P m+1 Γ ) T , then the linear system (5.3a-f), with numerical integration in (5.3e, f), can be written as ⎛ For the solution of (6.4) a Schur complement approach similar to [11] can be used. In particular, the Schur approach for eliminating ( κ m+1 , δ X m+1 , F m+1 Γ ) from (6.4) can be obtained as follows. Let Then (6.4) can be reduced to In (6.5a) we have used the definitions For the linear system (6.5a) well-known solution methods for finite element discretizations for the standard Navier-Stokes equations may be employed. We refer to [11, § 5], where we describe such solution methods in detail for a very similar situation.

Numerical results
For the bulk mesh adaptation we use the strategy from [11], which results in a fine mesh size h f around Γ m and a coarse mesh size h c further away from it. Here h f = To summarize the discretization parameters we use the shorthand notation n adapt k,l from [11]. The subscripts refer to the fineness of the spatial discretizations, i.e. for the set n adapt k,l it holds that N f = 2 k and N c = 2 l . For the case d = 2, in this paper, we have in addition that K Γ = J Γ = 2 k + 1. Finally, the uniform time step size for the set n adapt k,l is given by τ = 10 −3 /n, and if n = 1 we write adapt k,l . We remark that we implemented the scheme (5.3a-f) with the help of the finite element toolbox ALBERTA, see [37].
In all the numerical simulations we employ the scheme with numerical integration in (5.3e, f), i.e. we choose the superscript · h in the two brackets · (h) . Moreover, unless otherwise stated, we use the scheme (5.3a-f) with = 1. The initial data κ 0 ∈ V (Γ 0 ) is always computed as the solution to In addition, we employ the lowest order Taylor-Hood element P2-P1 in all computations and set U 0 = I 0 2 u 0 . Unless stated otherwise we fix ∂ 1 Ω = ∂Ω, g = 0 and u 0 = 0. The volume force is always set to f = 0. Moreover, we set all physical parameters to unity, i.e. ρ ± = μ ± = ρ Γ = μ Γ = α = 1, unless stated otherwise.
At times we will discuss the discrete energy of the numerical solutions. On recalling Theorem 4.1 the discrete energy is defined by represents the kinetic part of the discrete energy. For the simulation of vesicles the reduced volume is often mentioned as a characteristic number. In the case d = 3, and for the initial discrete interface Γ 0 , this is defined as 3 2 , see e.g. [42]. In a similar fashion, for the case d = 2 we define the reduced area as

Numerical simulations in 2D
For all our two-dimensional simulations we choose the discretization parameters 2 adapt 9,4 . In all the simulations presented here the areas of the two phases, as well as the length of the interface, are well preserved, with the relative differences over time in each case being less than 0.2 %. Moreover, the ratio of the largest and smallest elements' lengths was always bounded by 1.005. Here we note that we always choose the initial polygon Γ 0 to be equidistributed. We conducted the following shearing experiments on the domain Ω = (−2, 2) 2 for an initial interface in the form of an ellipse, centred at the origin, with axis lengths 1 and 2.5, so that a r = 0.745. In particular, we prescribe the inhomogeneous Dirichlet boundary condition g( z) = (z 2 , 0) T on ∂ 1 Ω = [−2, 2] × {±2}. For the initial data u 0 we choose the function u 0 ( z) = η(z 2 ) e 1 , where η : [−2, 2] → R is a continuous piecewise linear function with η(±2) = ±2 and η(s) = 0 if |s| ≤ 1.5. Hence u 0 satisfies the required conditions ∇ . u 0 = 0 in Ω and ∇ s . u 0 = 0 on Γ (0), recall (2.8), and is such that u 0 = g on ∂ 1 Ω. The remaining parameters are given by α = 0.05, ρ ± = ρ Γ = μ Γ = 1 and either (a) μ + = 1, μ − = 1 , or (b) μ + = 1, μ − = 10, (7.1) where we recall that the continuous problem, as d = 2, is independent of the value of μ Γ . The results can be seen in Figs. 2 and 3. In the first case we observe that the evolution reaches a steady state in which the interfacial fluid rotates along the interface. This motion is often called tank treading, see e.g. [36]. The second example, on the other hand, leads to a rotation of the whole vesicle, and this is called tumbling. The numerical simulation of a vesicle flowing through a constriction can be seen in Fig. 4. This example shows that membranes can drastically deform in order to pass through a constriction. This resembles the remarkable properties of red blood cells, which show a similar behaviour when flowing through capillaries. Here we choose the initial shape of the interface to be an elongated tube of total dimension 0.2 × 1.5. This gives a reduced area of a r = 0.351. As the computational domain we we prescribe the inhomogeneous boundary conditions g( z) = (1 − z 2 2 , 0) T in order to model Poiseuille flow. In this computation we consider the quasi-static variant with ρ ± = 0. Moreover, we set α = 0.1 and let ρ Γ = 0 or ρ Γ = 15. In the latter case the effect of inertia on the evolution is clearly visible. We note that for the experiment with ρ Γ = 0, the ratio r a increases at the very first time step from 1 to 1.0023, which explains the graph in Fig. 4. Here we recall from Remark 4.1 that equidistribution is only maintained for the semidiscrete variant of our scheme.
In order to highlight the local interface length preservation, we also show in Fig. 5 how a scalar field, initialized as Ψ 0 = 1, is transported along the interface by the fluid. To this end, at each time step, we find Ψ m+1 ∈ W (Γ m+1 ) such that recall (4.4e) in [9] without diffusion. We compare the results to a simulation for the scheme (5.3a-f) with = 2, when no local interface length preservation can be expected, recall Remark 4.2. Clearly, using the scheme with = 2 leads to a nonuniform distribution of the scalar field Ψ M , which coincides with a nonuniform distribution of the vertices along the discrete interface Γ M . Despite this difference in the approximation of Γ (t M ), the two simulations produce interfaces Γ M that are nearly identical. This suggests that the oscillatory surface pressure exhibited by the scheme (5.3a-f) with = 1 compared to = 2, see the bottom row in Fig. 5, does not have a detrimental effect on the velocity approximation. A very pronounced difference between ρ Γ = 0 and ρ Γ > 0 can be observed in our next simulation, where we start with an initial shape in the form of a smooth letter "C" with reduced area a r = 0.326. As the computational domain we choose Ω = (−1, 1) 2 , and we let ρ ± = α = 1. See Figs. 6 and 7 for the evolutions in the cases ρ Γ = 0 and ρ Γ = 1, respectively. In the latter case, the two arms of the vesicle swing up and down due to inertia, which is clearly visible in the plot of the kinetic energy as well.

Numerical simulations in 3D
Unless otherwise stated, for the uniform time step size we choose τ = 10 −3 in this subsection. In all the simulations presented here the volumes of the two phases, as well as the area of the interface, are almost exactly preserved, with the relative differences over time in each case being less than 0.2 %. As a first example for a three-dimensional simulation, we consider the evolution for an initially flat plate of total dimension 4 × 4 × 1, similarly to [6, Fig. 15]. As the computational domain we choose Ω = (−2.5, 2.5) 3 , and we set ρ Γ = 0. We note that the reduced volume for this shape is given by v r = 0.569. As discretization parameters we choose adapt 5,2 , and the initial triangulation Γ 0 satisfies (K Γ , J Γ ) = (1538, 3072) and r a = 1.898. The results for the scheme (5.3a-f) can be seen in Fig. 8, where we note that the interface assumes the shape of a red blood cell. Plots of the discrete energies and of the ratio r a are also shown in Fig. 8. We note that the discrete energy is monotonically decreasing, while the ratio r a always remains bounded below 2.
The numerical simulation of a vesicle flowing through a constriction can be seen in Fig. 9. Here we choose the initial shape of the interface to be a biconcave surface resembling a human red blood cell, with a reduced volume of v r = 0.568. As the computational domain we choose Ω = (−2, −1)×(−1, 1) 2 ∪[−1, 1]×(−0.5, 0.5) 2 ∪ (1, 2) × (−1, 1) 2 . We define ∂ 2 Ω = {2} × (−1, 1) 2 and on ∂ 1 Ω we set no-slip Similarly to in Fig. 5, we consider the transport of a scalar field Ψ 0 = 1 on Γ 0 for the simulation in Fig. 9. As can be seen from the plot in Fig. 10, the local surface area is maintained almost exactly throughout the evolution. The same cannot be said for the scheme (5.3a-f) with = 2. Here the final distribution of Ψ M is very uneven, and the overall surface area decreases by 14 %. As a consequence, the final shape Γ M differs dramatically from the run with = 1.
We also show the final surface pressure P M Γ in Fig. 10, and also compare the results with a run for the scheme (5.3a-f) with = 2. Similarly to the results in Fig. 5, we observe oscillatory surface pressures when = 1, while = 2 yields well behaved surface pressure approximations. The same simulation for piecewise constant surface pressures P m+1 ∈ S m 0 (Γ m ), and for both = 1 and = 2, yields the plots in Fig. 11. Here we replace (7.2) with the suitable formulation for the basis functions of S m 0 (Γ m ) and S m+1 0 (Γ m+1 ). We observe that this evolution can be computed also for piecewise constant surface pressure approximations. However, for more complex evolutions, such as three dimensional analogues of Figs. 2 and 3, we observe locking. Here the iterative linear solver is unable to find a discrete solution, even at the first time step.