A Finite Element Method for Two-Phase Flow with Material Viscous Interface

Maxim Olshanskii; Annalisa Quaini; Qi Sun

doi:10.1515/cmam-2021-0185

Article Publicly Available

A Finite Element Method for Two-Phase Flow with Material Viscous Interface

Maxim Olshanskii , Annalisa Quaini and Qi Sun

Published/Copyright: November 3, 2021

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

From the journal Computational Methods in Applied Mathematics Volume 22 Issue 2

Abstract

This paper studies a model of two-phase flow with an immersed material viscous interface and a finite element method for the numerical solution of the resulting system of PDEs. The interaction between the bulk and surface media is characterized by no-penetration and slip with friction interface conditions. The system is shown to be dissipative, and a model stationary problem is proved to be well-posed. The finite element method applied in this paper belongs to a family of unfitted discretizations. The performance of the method when model and discretization parameters vary is assessed. Moreover, an iterative procedure based on the splitting of the system into bulk and surface problems is introduced and studied numerically.

Keywords: Two-Phase Flow; Surface Stokes Equations; Finite Elements; CutFEM; TraceFEM

MSC 2010: 65N12; 65N30; 76T30

1 Introduction

Motivated by applications in biomedical engineering [59] and cell biology [61], we are interested in the interaction between a rather viscous and dense interface material and a less viscous and dense bulk fluid, as is the case of, e.g., a lipid vesicle or cell membrane with its content and surrounding fluid. For this purpose, we consider two immiscible, viscous, and incompressible fluids separated by a viscous inextensible material interface modeled as a Boussinesq–Scriven surface fluid. The following coupling conditions are prescribed between the bulk two-phase flow and the surface fluid: (i) the immiscibility condition, i.e. the bulk fluid does not penetrate through the interface, (ii) slip with friction between the bulk fluid and the viscous interface, and (iii) the load exerted from the bulk fluid onto the surface fluid defined by the jump of the normal stress across the interface. To the best of our knowledge, it is the first time that the slip with friction conditions are imposed in this context. This extends the Navier boundary condition to the interface case and interpolates between no-slip of the bulk phase and embedded material layer, as considered in [2, 8, 44, 51], and uncoupled lateral dynamics of fluidic surfaces as studied, for example, in [52, 54, 57]. We remark that the classical Boussinesq–Scriven interface model as presented in, e.g., [8, 51] can be used to account for the viscous effects in fluidic membranes. Its extension is needed to address complex physical phenomena known to be relevant for lipid membranes, which include fluidity [18, 33], inextensibility [27, 41], inertia [29, 40], and slip with respect to the surrounding fluid [1]. By accounting for all of the above phenomena, our model represents another step towards a more realistic representation of lipid membranes interacting with bulk fluid.

We first state the problem in terms of incompressible two-phase Navier–Stokes flow in the bulk coupled to surface fluid equations and show the energy balance for this system. This balance demonstrates that the system with the proposed coupling conditions is thermodynamically consistent, i.e. in the absence of external forces and with no energy inflow through the boundary the system is dissipative. Then we introduce a simplifying assumption: the coupled system has reached a steady state and inertia terms can be neglected. Although the resulting surface-bulk Stokes problem is a strong simplification of the original problem, it still retains interesting features for the purpose of numerical analysis.

For the simplified model, we show well-posedness following the abstract framework for generalized saddle point problems proposed in, e.g., [3, 43]. Next, we present a sharp interface, geometrically unfitted finite element (FE) method. Unfitted methods allow the sharp interface to cut through the elements of a fixed background grid. Their main advantage is the relative ease in handling time-dependent domains, implicitly defined interfaces, and problems with strong geometric deformations [7]. For the bulk two-phase Stokes problem, we apply an isoparametric unfitted finite element approach [35] of the CutFEM (or Nitsche-XFEM) family [12, 15, 25, 26, 38, 58]. For a (bulk only) interface Stokes problem with slip between phases, the approach was studied in [46]. CutFEM uses overlapping fictitious domains in combination with ghost penalty stabilization [9] to enrich and stabilize the solution. In this paper, we consider the unfitted generalized Taylor–Hood finite element pair 𝐏 k + 1 - P k , k ≥ 1 . For more details on the isoparametric unfitted finite element, we refer to [36, 37]. For the discretization of the surface Stokes problem, we apply the trace finite element method [48, 31], which uses traces of the bulk finite element functions. In a set of numerical experiments, we address the dependence of the discretization error of the method on variations in the friction coefficient, ratio of bulk fluid viscosity, surface fluid viscosity, mesh size, and position of the interface relative to the fixed computational mesh.

One option to solve the coupled bulk-surface problem is to adopt a monolithic approach. However, such an approach would quickly lead to high computational costs as the mesh gets refined. Instead, we introduce a partitioned scheme based on fixed-point iterations [50]. With this algorithm, the bulk and surface flow problems are solved separately and sequentially, and the coupling conditions are enforced in an iterative fashion. We choose this algorithm for its simplicity of implementation and study its performance when model and discretization parameters vary. Our numerical experiments provide evidence of the robustness of the proposed approach with respect to the contrast in viscosity in the bulk fluid, surface fluid viscosity, and position of the interface relative to the background mesh. At the same time, convergence slowdown was observed for certain values of the slip coefficients.

The outline of the paper is as follows. In Section 2, we introduce the strong formulation of the coupled problem and the associated energy balance. Section 3 presents the strong and weak formulations of the simplified problem, together with the finite element discretization. In Section 4, we propose a partitioned algorithm for the numerical solution of the coupled problem. Numerical results in three dimensions are reported in Section 5. For all the simulations in this paper, we have used NGsolve [62, 17], a high performance multiphysics finite element software with a Python interface, and add-on library ngsxfem [63].

2 A Two-Phase Fluid with Material Viscous Interface

Consider a fixed volume Ω ⊂ ℝ 3 filled with two immiscible, viscous, and incompressible fluids separated by an interface Γ ⁢ ( t ) for all t ∈ [ 0 , T ] . We assume that Γ ⁢ ( t ) stays closed and sufficiently smooth (at least C 2 ) for all t ∈ [ 0 , T ] . Surface Γ ⁢ ( t ) separates Ω into two phases (subdomains) Ω + ⁢ ( t ) and Ω - ( t ) : = Ω ∖ Ω + ⁢ ( t ) ¯ . We assume Ω - ⁢ ( t ) to be completely internal, i.e. ∂ ⁡ Ω - ⁢ ( t ) ∩ ∂ ⁡ Ω = ∅ for all times. See Figure 1. Denote by 𝐧 ± the outward normals for Ω ± ⁢ ( t ) and 𝐧 the normal on Γ pointing from Ω - ⁢ ( t ) to Ω + ⁢ ( t ) ; it holds that 𝐧 - = 𝐧 and 𝐧 + = - 𝐧 at Γ. For ease of notation, from now on, we will drop the dependence on t for Γ, Ω + , and Ω - .

$Figure 1 Fluid domains and interface in ℝ 2 {\mathbb{R}^{2}} .$

Figure 1

Fluid domains and interface in ℝ 2 .

Consider the bulk fluid velocities 𝐮 ± ⁢ ( t ) : Ω ± → ℝ 3 and pressures p ± ⁢ ( t ) : Ω ± → ℝ to describe the fluid motion. The motion of the fluids occupying subdomains Ω ± is governed by the incompressible Navier–Stokes equations

(2.1) ρ ± ⁢ ∂ ˙ t ⁢ 𝐮 ± = div ⁡ 𝝈 ± + 𝐟 ± ⁢ in ⁢ Ω ± ,

(2.2) div ⁡ 𝐮 ± = 0 ⁢ in ⁢ Ω ± ,

for all t ∈ ( 0 , T ) . In (2.1), constants ρ ± represent the fluid density, ∂ ˙ t denotes the material derivative, 𝐟 ± are the external body forces, and 𝝈 ± are the Cauchy stress tensors. For Newtonian fluids in both phases, the Cauchy stress tensors have the following expression:

𝝈 ± = - p ± ⁢ 𝐈 + 2 ⁢ μ ± ⁢ 𝐃 ⁢ ( 𝐮 ± ) , 𝐃 ⁢ ( 𝐮 ± ) = 1 2 ⁢ ( ∇ ⁡ 𝐮 ± + ( ∇ ⁡ 𝐮 ± ) T ) ⁢ in ⁢ Ω ± ,

where constants μ ± represent the fluid dynamic viscosity.

We assume interface Γ to be a thin material layer with possibly different material properties from the bulk fluid. Motivated by applications in cell biology, we consider a viscous inextensible interface modeled as an “incompressible” surface fluid. The evolution of the material interface can be described in terms of the velocity of this surface fluid denoted by 𝐔 . Later, we will need the decomposition of 𝐔 into tangential and normal components: 𝐔 = 𝐔 T + U N ⁢ 𝐧 , with 𝐔 T ⋅ 𝐧 = 0 , U N = 𝐧 ⋅ 𝐔 . The surface Navier–Stokes equations governing the motion of a fluidic deformable layer appear in several works [60, 55, 30, 39]. Here, we adopt the formulation in terms of tangential differential operators from [30], where the equations have been derived from conservation principles. To introduce these equations, we need some further definitions. Let 𝐏 ( 𝐱 ) : = 𝐈 - 𝐧 ( 𝐱 ) 𝐧 ( 𝐱 ) T for 𝐱 ∈ Γ be the orthogonal projection onto the tangent plane. For a scalar function π : Γ → ℝ or a vector function 𝐔 : Γ → ℝ 3 , we define π e : 𝒪 ⁢ ( Γ ) → ℝ , 𝐔 e : 𝒪 ⁢ ( Γ ) → ℝ 3 , smooth extensions of π and 𝐔 from Γ to its neighborhood 𝒪 ⁢ ( Γ ) . The surface gradient and covariant derivatives on Γ are then defined as ∇ Γ ⁡ π = 𝐏 ⁢ ∇ ⁡ π e and ∇ Γ 𝐔 : = 𝐏 ∇ 𝐔 e 𝐏 . The definitions of ∇ Γ ⁡ π and ∇ Γ ⁡ 𝐔 are independent of the particular smooth extension of π and 𝐔 off Γ. On Γ, we consider the surface rate-of-strain tensor [22] given by

𝐃 Γ ( 𝐔 ) : = 1 2 𝐏 ( ∇ 𝐔 + ( ∇ 𝐔 ) T ) 𝐏 = 1 2 ( ∇ Γ 𝐔 + ( ∇ Γ 𝐔 ) T ) .

The surface divergence for a vector 𝐠 : Γ → ℝ 3 and a tensor 𝐀 : Γ → ℝ 3 × 3 are defined as: div Γ 𝐠 : = tr ( ∇ Γ 𝐠 ) , div Γ 𝐀 : = ( div Γ ( 𝐞 1 T 𝐀 ) , div Γ ( 𝐞 2 T 𝐀 ) , div Γ ( 𝐞 3 T 𝐀 ) ) T , where 𝐞 i is the ith standard basis vector. With this notation, the surface Navier–Stokes equations take the form

(2.3) ρ Γ ⁢ ∂ ˙ t ⁢ 𝐔 = - ∇ Γ ⁡ π + 2 ⁢ μ Γ ⁢ div Γ ⁡ 𝐃 Γ ⁢ ( 𝐔 ) + 𝐟 Γ + 𝐛 e + π ⁢ κ ⁢ 𝐧 ⁢ on ⁢ Γ ,

(2.4) div Γ ⁡ 𝐔 = 0 ⁢ on ⁢ Γ ,

where ρ Γ is the surface fluid density, μ Γ is the surface fluid dynamic viscosity, κ denotes pointwise doubled mean curvature on Γ, and π is the surface fluid pressure. The material derivative in (2.3) is taken with respect to surface fluid trajectories, i.e. ∂ ˙ t ⁢ 𝐔 = ∂ ⁡ 𝐔 ∂ ⁡ t + ( 𝐔 ⋅ ∇ ) ⁢ 𝐔 . Note that ∂ ˙ t ⁢ 𝐔 is an intrinsic surface quantity, although both terms ∂ ⁡ 𝐔 ∂ ⁡ t and ( 𝐔 ⋅ ∇ ) ⁢ 𝐔 depend on extension of 𝐔 in the bulk. On the right-hand side of (2.3), 𝐟 Γ denotes the external area force acting on the surface as a result of the interaction with the bulk fluids (specified below), while 𝐛 e denotes other possible area force (such as elastic bending forces) and is not specified further in the present paper.

Next, we turn to the coupling conditions between equations (2.1)–(2.2) posed in the bulk and equations (2.3)–(2.4) posed on Γ. First, the immiscibility condition means that the bulk fluid does not penetrate through Γ, which implies that

(2.5) 𝐮 + ⋅ 𝐧 = U N = 𝐮 - ⋅ 𝐧 on ⁢ Γ .

Normal velocity U N determines radial deformations of Γ ⁢ ( t ) , and so it governs the geometric evolution of the interface, which can be defined through the Lagrangian mapping Ψ ⁢ ( t , ⋅ ) from Γ ⁢ ( 0 ) to Γ ⁢ ( t ) : for 𝐱 ∈ Γ ⁢ ( 0 ) , Ψ ⁢ ( t , 𝐱 ) solves the ODE system

(2.6) Ψ ⁢ ( 0 , 𝐱 ) = 𝐱 , ∂ ⁡ Ψ ⁢ ( t , 𝐱 ) ∂ ⁡ t = U N ⁢ ( t , Ψ ⁢ ( t , 𝐱 ) ) , t ∈ [ 0 , T ] .

In fluid vesicles and cells, typically a viscous and dense lipid membrane, represented by Γ, is surrounded by a less viscous and less dense liquid. We are interested in modeling slip with friction between the bulk fluid and the viscous membrane. Thus, we consider Navier-type conditions

(2.7) 𝐏 ⁢ 𝝈 + ⁢ 𝐧 = f + ⁢ ( 𝐏𝐮 + - 𝐔 T ) ⁢ on ⁢ Γ ,

(2.8) 𝐏 ⁢ 𝝈 - ⁢ 𝐧 = - f - ⁢ ( 𝐏𝐮 - - 𝐔 T ) ⁢ on ⁢ Γ ,

where f - and f + are friction coefficients at Γ on the Ω - and Ω + side, respectively. Conditions (2.7)–(2.8) model an incomplete adhesion of a bulk fluid to the material surface with 1 / f ± often referred to as a “slip length” [42]; see, e.g., [4, 34] for the modern description of experimental and theoretical validations. In particular, the acceptance of nonzero slip length resolves the well-known “no-collision paradox” [13, 28, 19], thus suggesting (2.7)–(2.8) to be an important modeling assumption in the simulation of a lipid vesicle – cell membrane contact (and fusion). We finally note that the Navier conditions are not an uncommon choice in numerical models if the flow in boundary region is under-resolved [32].

The area force in (2.3) coming from the bulk fluid is defined by the jump of the normal stress on Γ,

(2.9) 𝐟 Γ = [ 𝝈 ⁢ 𝐧 ] - + = 𝝈 + ⁢ 𝐧 - 𝝈 - ⁢ 𝐧 on ⁢ Γ .

We summarize the complete system of equations and coupling conditions below:

(2.10) { ρ ± ⁢ ∂ ˙ t ⁢ 𝐮 ± - μ ± ⁢ Δ ⁢ 𝐮 ± + ∇ ⁡ p ± = 𝐟 ± in ⁢ Ω ± , div ⁡ 𝐮 ± = 0 in ⁢ Ω ± , ρ Γ ⁢ ∂ ˙ t ⁢ 𝐔 - 2 ⁢ μ Γ ⁢ div Γ ⁡ 𝐃 Γ ⁢ ( 𝐔 ) + ∇ Γ ⁡ π - π ⁢ κ ⁢ 𝐧 = [ 𝝈 ⁢ 𝐧 ] - + + 𝐛 e on ⁢ Γ , div Γ ⁡ 𝐔 = 0 on ⁢ Γ , 𝐮 + ⋅ 𝐧 = 𝐮 - ⋅ 𝐧 = U N on ⁢ Γ , 𝐏 ⁢ 𝝈 ± ⁢ 𝐧 = ± f ± ⁢ ( 𝐏𝐮 ± - 𝐔 T ) on ⁢ Γ ,

with the evolution of Ω ± and Γ defined by the velocity solving the system through (2.6). On ∂ ⁡ Ω , the system is endowed with boundary conditions either for the bulk velocity or for the bulk normal stress,

(2.11) 𝐮 + = 𝐠 ⁢ on ⁢ ∂ ⁡ Ω D ,

(2.12) 𝝈 + ⁢ 𝐧 + = 𝐟 N ⁢ on ⁢ ∂ ⁡ Ω N .

Here ∂ ⁡ Ω D ¯ ∪ ∂ ⁡ Ω N ¯ = ∂ ⁡ Ω ¯ and ∂ ⁡ Ω D ∩ ∂ ⁡ Ω N = ∅ . See Figure 1. At t = 0 , initial velocity is given by 𝐮 ± = 𝐮 0 ± in Ω ± ⁢ ( 0 ) and 𝐔 = 𝐔 0 on Γ ⁢ ( 0 ) .

2.1 Energy Balance of the Continuous Coupled Problem

We look for the energy balance of the coupled system (2.10)–(2.12). We make use of the following identities for time-dependent domains Ω ± ⁢ ( t ) , whose only moving part of the boundary is Γ ⁢ ( t ) :

(2.13) 1 2 ⁢ d d ⁢ t ⁢ ∫ Ω ± ⁢ ( t ) | 𝐮 ± | 2 ⁢ 𝑑 V = ∫ Ω ± ⁢ ( t ) 𝐮 ± ⋅ ∂ ⁡ 𝐮 ± ∂ ⁡ t ⁢ 𝑑 V + 1 2 ⁢ ∫ Γ ⁢ ( t ) | 𝐮 ± | 2 ⁢ 𝐮 ± ⋅ 𝐧 ± ⁢ 𝑑 S ,

(2.14) ∫ Ω ± ⁢ ( t ) ( 𝐮 ± ⋅ ∇ ) ⁢ 𝐮 ± ⋅ 𝐮 ± ⁢ 𝑑 V = 1 2 ⁢ ∫ ∂ ⁡ Ω ± ⁢ ( t ) | 𝐮 ± | 2 ⁢ 𝐮 ± ⋅ 𝐧 ± ⁢ 𝑑 S .

Identity (2.13) is the Reynolds transport theorem, while identity (2.14) is obtained from integration by parts.

Let us start from the kinetic energy of the fluid in Ω - ,

(2.15) d d ⁢ t ⁢ E - = 1 2 ⁢ d d ⁢ t ⁢ ∫ Ω - ρ - ⁢ | 𝐮 - | 2 ⁢ 𝑑 V = ∫ Ω - ρ - ⁢ 𝐮 - ⋅ ∂ ⁡ 𝐮 - ∂ ⁡ t ⁢ 𝑑 V + 1 2 ⁢ ∫ Γ ρ - ⁢ | 𝐮 - | 2 ⁢ 𝐮 - ⋅ 𝐧 ⁢ 𝑑 S = ∫ Ω - 𝐮 - ⋅ ( div ⁡ 𝝈 - + 𝐟 - - ρ - ⁢ ( 𝐮 - ⋅ ∇ ) ⁢ 𝐮 - ) ⁢ 𝑑 V + 1 2 ⁢ ∫ Γ ρ - ⁢ | 𝐮 - | 2 ⁢ 𝐮 - ⋅ 𝐧 ⁢ 𝑑 S = ∫ Ω - ∇ ⋅ ( 𝝈 - ⁢ 𝐮 - ) ⁢ 𝑑 V - ∫ Ω - 𝝈 - : ∇ ⁡ 𝐮 - ⁢ d ⁢ V + ∫ Ω - 𝐮 - ⋅ 𝐟 - ⁢ 𝑑 V = ∫ Γ 𝐮 - ⋅ ( 𝝈 - ⁢ 𝐧 ) ⁢ 𝑑 S - ∫ Ω - 𝝈 - : 𝐃 ⁢ ( 𝐮 - ) ⁢ d ⁢ V + ∫ Ω - 𝐮 - ⋅ 𝐟 - ⁢ 𝑑 V = ∫ Γ 𝐮 - ⋅ ( 𝝈 - ⁢ 𝐧 ) ⁢ 𝑑 S - 2 ⁢ μ - ⁢ ∫ Ω - ∥ 𝐃 ⁢ ( 𝐮 - ) ∥ 2 ⁢ 𝑑 V + ∫ Ω - 𝐮 - ⋅ 𝐟 - ⁢ 𝑑 V .

Above, we have used (2.13), (2.1), (2.2), (2.14), and integration by parts.

We repeat similar steps for the kinetic energy of the fluid in Ω + , the main difference being that

∂ ⁡ Ω + ¯ = Γ ¯ ∪ ∂ ⁡ Ω D ¯ ∪ ∂ ⁡ Ω N ¯ ,

while ∂ ⁡ Ω - ¯ = Γ ¯ . We obtain

(2.16) d d ⁢ t ⁢ E + = - ∫ Γ 𝐮 + ⋅ ( 𝝈 + ⁢ 𝐧 ) ⁢ 𝑑 S - 2 ⁢ μ + ⁢ ∫ Ω + ∥ 𝐃 ⁢ ( 𝐮 + ) ∥ 2 ⁢ 𝑑 V + ∫ Ω + 𝐮 + ⋅ 𝐟 + ⁢ 𝑑 V + B ,

where

(2.17) B = ∫ ∂ ⁡ Ω D ( 𝐠 ⋅ ( 𝝈 + ⁢ 𝐧 ) - 1 2 ⁢ ρ + ⁢ | 𝐠 | 2 ⁢ 𝐠 ⋅ 𝐧 + ) ⁢ 𝑑 S + ∫ ∂ ⁡ Ω N ( 𝐮 + ⋅ 𝐟 - 1 2 ⁢ ρ + ⁢ | 𝐮 + | 2 ⁢ 𝐮 + ⋅ 𝐧 + ) ⁢ 𝑑 S .

Putting together (2.15) and (2.16), we obtain the total kinetic energy for the bulk flow,

(2.18) d ⁢ E d ⁢ t = d d ⁢ t ⁢ ( E + + E - ) = - 2 ⁢ μ - ⁢ ∫ Ω - ∥ 𝐃 ⁢ ( 𝐮 - ) ∥ 2 ⁢ 𝑑 V - 2 ⁢ μ + ⁢ ∫ Ω + ∥ 𝐃 ⁢ ( 𝐮 + ) ∥ 2 ⁢ 𝑑 V - ∫ Γ U N ⁢ [ 𝐧 T ⁢ σ ⁢ 𝐧 ] - + ⁢ 𝑑 S - ∫ Γ f - ⁢ 𝐏𝐮 - ⋅ ( 𝐏𝐮 - - 𝐔 T ) ⁢ 𝑑 S - ∫ Γ f + ⁢ 𝐏𝐮 + ⋅ ( 𝐏𝐮 + - 𝐔 T ) ⁢ 𝑑 S + ∫ Ω - 𝐮 - ⋅ 𝐟 - ⁢ 𝑑 V + ∫ Ω + 𝐮 + ⋅ 𝐟 + ⁢ 𝑑 V + B ,

where we have used (2.5), (2.7), and (2.8). For the kinetic energy of the surface flow, we will use the surface analogue of the Reynolds transport theorem (see, e.g., [14, Lemma 2.1]),

(2.19) d d ⁢ t ⁢ ∫ Γ f ⁢ 𝑑 S = ∫ Γ ( ∂ ˙ t ⁢ f + f ⁢ div Γ ⁡ 𝐔 ) ⁢ 𝑑 S .

We also need the surface integration by part identity

∫ Γ ( g ⁢ div Γ ⁡ 𝐟 + 𝐟 ⋅ ∇ Γ ⁡ g ) ⁢ 𝑑 S = ∫ Γ κ ⁢ g ⁢ ( 𝐟 ⋅ 𝐧 ) ⁢ 𝑑 S ,

which is valid for a smooth closed surface with smooth scalar field g and vector field 𝐟 . A componentwise application of this equality yields the identity

(2.20) ∫ Γ 𝐠 ⋅ div Γ ⁡ 𝐆 ⁢ d ⁢ S = - ∫ Γ 𝐆 : ∇ Γ ⁡ 𝐠 ⁢ d ⁢ S

for a vector field 𝐠 ∈ ( C 1 ⁢ ( Γ ) ) 3 and a matrix function 𝐆 ∈ ( C 1 ⁢ ( Γ ) ) 3 × 3 such that 𝐆 = 𝐏𝐆𝐏 .

The energy balance on Γ is given by

(2.21) d ⁢ E Γ d ⁢ t = 1 2 ⁢ d d ⁢ t ⁢ ∫ Γ ρ Γ ⁢ | 𝐔 | 2 ⁢ 𝑑 S = ∫ Γ ( ρ Γ ⁢ 𝐔 ⋅ ∂ ˙ t ⁢ 𝐔 + 1 2 ⁢ ρ Γ ⁢ | 𝐔 | 2 ⁢ div Γ ⁡ 𝐔 ) ⁢ 𝑑 S = ∫ Γ 𝐔 ⋅ ( - ∇ Γ ⁡ π + 2 ⁢ μ Γ ⁢ div Γ ⁡ 𝐃 Γ ⁢ ( 𝐔 ) + 𝐟 Γ + 𝐛 e ) ⁢ 𝑑 S = - 2 ⁢ μ Γ ⁢ ∫ Γ 𝐃 Γ ⁢ ( 𝐔 ) : ∇ Γ ⁡ 𝐔 ⁢ d ⁢ S + ∫ Γ 𝐔 ⋅ [ 𝝈 ⁢ 𝐧 ] - + ⁢ 𝑑 S + ∫ Γ 𝐔 T ⋅ 𝐛 e ⁢ 𝑑 S = - 2 ⁢ μ Γ ⁢ ∫ Γ ∥ 𝐃 Γ ⁢ ( 𝐔 ) ∥ 2 ⁢ 𝑑 S + ∫ Γ U N ⁢ [ 𝐧 T ⁢ 𝝈 ⁢ 𝐧 ] - + ⁢ 𝑑 S + ∫ Γ 𝐔 T ⋅ ( f + ⁢ ( 𝐏𝐮 + - 𝐔 T ) + f - ⁢ ( 𝐏𝐮 - - 𝐔 T ) ) ⁢ 𝑑 S + ∫ Γ 𝐔 ⋅ 𝐛 e ⁢ 𝑑 S ,

where we have applied (2.20) and used (2.19), (2.3)–(2.5), (2.7)–(2.9), and integration by parts.

Combining (2.18) and (2.21), we get the kinetic energy for the bulk and surface flows,

(2.22) d ⁢ E d ⁢ t + d ⁢ E Γ d ⁢ t = - 2 ⁢ μ - ⁢ ∫ Ω - ∥ 𝐃 ⁢ ( 𝐮 - ) ∥ 2 ⁢ 𝑑 V - 2 ⁢ μ + ⁢ ∫ Ω + ∥ 𝐃 ⁢ ( 𝐮 + ) ∥ 2 ⁢ 𝑑 V ⏟ bulk fluid viscous dissipation ⁢ - 2 ⁢ μ Γ ⁢ ∫ Γ ∥ 𝐃 Γ ⁢ ( 𝐔 T ) ∥ 2 ⁢ 𝑑 S ⏟ surface fluid viscous dissipation - ∫ Γ f - ⁢ ∥ 𝐏𝐮 - - 𝐔 T ∥ 2 ⁢ 𝑑 S - ∫ Γ f + ⁢ ∥ 𝐏𝐮 + - 𝐔 T ∥ 2 ⁢ d ⏟ frictional energy dissipation + ∫ Γ 𝐮 - ⋅ 𝐟 - ⁢ 𝑑 V + ∫ Γ 𝐮 + ⋅ 𝐟 + ⁢ 𝑑 V + ∫ Γ 𝐔 ⋅ 𝐛 e ⁢ 𝑑 S ⏟ work of external forces ⁢ + B ⏟ work of b.c. ,

where B, i.e. the work of the boundary conditions, is defined in (2.17).

From (2.22), we see that, in the absence of external forces and with no energy inflow through the boundary, the system is dissipative, i.e. thermodynamically consistent.

3 A Simplified Steady Problem

In this section, we consider a (strongly) simplified version of the problem presented in Section 2. Our main assumption is that the coupled bulk and surface fluid system has reached a steady state and inertia terms can be neglected. Since the steady state implies Γ ⁢ ( t ) = Γ ⁢ ( 0 ) , we have U N = 0 and hence 𝐔 = 𝐔 T . This simplified surface-bulk Stokes problem models a viscosity dominated two-phase flow with the viscous interface in a dynamical equilibrium; see also Remark 1. It is an interesting model problem for the purpose of numerical analysis. With these simplifications, equations (2.1)–(2.2) become

(3.1) - μ ± ⁢ Δ ⁢ 𝐮 ± + ∇ ⁡ p = 𝐟 ± ⁢ in ⁢ Ω ± ,

(3.2) div ⁡ 𝐮 ± = 0 ⁢ in ⁢ Ω ± .

We impose a non-homogeneous Dirichlet condition on the entire outer boundary of Ω, i.e. problem (3.1)–(3.2) is supplemented with boundary condition (2.11) with 𝐠 ∈ [ H 1 / 2 ⁢ ( ∂ ⁡ Ω D ) ] 3 and ∂ ⁡ Ω D = ∂ ⁡ Ω . Under our assumption, the momentum equation for the surface fluid simplifies to - 2 ⁢ μ Γ ⁢ div Γ ⁡ 𝐃 Γ ⁢ ( 𝐔 T ) + ∇ Γ ⁡ π - π ⁢ κ ⁢ 𝐧 = [ 𝝈 ⁢ 𝐧 ] - + + 𝐛 e . The tangential part of the above momentum equation together with the inextensibility condition (2.4) leads to the surface Stokes problem

(3.3) - 2 ⁢ μ Γ ⁢ 𝐏 ⁢ div Γ ⁡ 𝐃 Γ ⁢ ( 𝐔 T ) + ∇ Γ ⁡ π = [ 𝐏 ⁢ 𝝈 ⁢ 𝐧 ] - + + 𝐏𝐛 e ⁢ on ⁢ Γ ,

(3.4) div Γ ⁡ 𝐔 T = 0 ⁢ on ⁢ Γ ,

while the normal part simplifies to

(3.5) [ 𝐧 T ⁢ 𝝈 ⁢ 𝐧 ] + - = π ⁢ κ on ⁢ Γ .

The interface condition above is standard in many models of two-phase flows, where π has the meaning of the surface tension coefficient.

Coupling condition (2.5) is replaced by

(3.6) 𝐮 + ⋅ 𝐧 = 𝐮 - ⋅ 𝐧 on ⁢ Γ ,

while conditions (2.7) and (2.8) still hold,

(3.7) 𝐏 ⁢ 𝝈 ± ⁢ 𝐧 = ± f ± ⁢ ( 𝐏𝐮 ± - 𝐔 T ) on ⁢ Γ .

Finally, we will see that the (weak formulation of the) problem is well-posed under two mean conditions for the bulk pressure,

∫ Ω ± p ± ⁢ 𝑑 x = 0 .

Remark 1.

Since 𝐔 ⋅ 𝐧 = 0 , condition (3.6) allows the flow through the steady interface Γ. This is inconsistent with (2.5), which assumes immiscibility of fluids. For a physically consistent formulation that describes the true equilibrium, one has to set 𝐮 + ⋅ 𝐧 = 𝐮 - ⋅ 𝐧 = 0 on Γ, but allow the shape of Γ to be the unknown, i.e. to be determined as a part of the problem. For such equilibrium to exist, external forces and boundary conditions may have to satisfy additional constraints. Finding such constraints and solving the resulting nonlinear problem is outside the scope of this paper. We rather follow a common convention in the analysis of models for steady two-phase problems and allow (3.6) for the steady interface; see, e.g., [16, 47, 21, 6].

3.1 Variational Formulation

The purpose of this section is to derive the variational formulation of coupled problem (3.1)–(3.7). Let us introduce some standard notation. The space of functions whose square is integrable in a domain ω is denoted by L 2 ⁢ ( ω ) . The space of functions whose distributional derivatives of order up to m ≥ 0 (integer) belong to L 2 ⁢ ( ω ) is denoted by H m ⁢ ( ω ) . The space of vector-valued functions with components in L 2 ⁢ ( ω ) is denoted by L 2 ⁢ ( ω ) 3 , and H 1 ⁢ ( div , ω ) is the space of functions in L 2 ⁢ ( ω ) with divergence in L 2 ⁢ ( ω ) . Moreover, we introduce the following functional spaces:

V - = H 1 ( Ω - ) 3 , V + = { 𝐮 ∈ H 1 ( Ω + ) 3 , 𝐮 | ∂ ⁡ Ω D = 𝐠 } , V 0 + = { 𝐮 ∈ H 1 ( Ω + ) 3 , 𝐮 | ∂ ⁡ Ω D = 𝟎 } ,

V ± = { 𝐮 = ( 𝐮 - , 𝐮 + ) ∈ V - × V + , 𝐮 - ⋅ 𝐧 = 𝐮 + ⋅ 𝐧 on Γ } ,

V 0 ± = { 𝐮 = ( 𝐮 - , 𝐮 + ) ∈ V - × V 0 + , 𝐮 - ⋅ 𝐧 = 𝐮 + ⋅ 𝐧 on Γ } ,

L ± 2 = { p = ( p - , p + ) ∈ L 2 ( Ω - ) × L 2 ( Ω + ) : ∫ Ω ± p ± d x = 0 } ,

V Γ = { 𝐔 ∈ H 1 ⁢ ( Γ ) 3 : 𝐔 ⋅ 𝐧 = 0 } .

The space V ± can be also characterized as ( V - × V + ) ∩ H 1 ⁢ ( div , Ω ) . We use ( ⋅ , ⋅ ) ω and 〈 ⋅ , ⋅ 〉 ω to denote the L 2 product and the duality pairing, respectively.

Multiplying (3.1) by 𝐯 ∈ V 0 ± and (3.2) by q ∈ L 0 2 ⁢ ( Ω ) and integrating over each subdomain, we see that smooth bulk velocity and pressure satisfy integral identity,

(3.8) - ( p - , div ⁡ 𝐯 - ) Ω - - ( p + , div ⁡ 𝐯 + ) Ω + + 2 ⁢ ( μ - ⁢ 𝐃 ⁢ ( 𝐮 - ) , 𝐃 ⁢ ( 𝐯 - ) ) Ω - + 2 ⁢ ( μ + ⁢ 𝐃 ⁢ ( 𝐮 + ) , 𝐃 ⁢ ( 𝐯 + ) ) Ω + - 〈 π ⁢ κ , 𝐯 - ⋅ 𝐧 〉 Γ + 〈 f - ⁢ ( 𝐏𝐮 - - 𝐔 ) , 𝐏𝐯 - 〉 Γ + 〈 f + ⁢ ( 𝐏𝐮 + - 𝐔 ) , 𝐏𝐯 + 〉 Γ = ( 𝐟 - , 𝐯 - ) Ω - + ( 𝐟 + , 𝐯 + ) Ω +

(3.9) ( div ⁡ 𝐮 - , q - ) Ω - + ( div ⁡ 𝐮 + , q + ) Ω + = 0

for all ( 𝐯 , q ) ∈ V 0 ± × L 0 2 ⁢ ( Ω ) . The interface terms in (3.8) have been obtained using coupling conditions (3.7) and (3.5) as follows:

- 〈 𝝈 - ⁢ 𝐧 , 𝐯 - 〉 Γ + 〈 𝝈 + ⁢ 𝐧 , 𝐯 + 〉 Γ = - 〈 𝐏 ⁢ 𝝈 - ⁢ 𝐧 , 𝐏𝐯 - 〉 Γ + 〈 𝐏 ⁢ 𝝈 + ⁢ 𝐧 , 𝐏𝐯 + 〉 Γ - 〈 [ 𝐧 T ⁢ 𝝈 ⁢ 𝐧 ] + - , 𝐯 - ⋅ 𝐧 〉 Γ = 〈 f - ⁢ ( 𝐏𝐮 - - 𝐔 ) , 𝐏𝐯 - 〉 Γ + 〈 f + ⁢ ( 𝐏𝐮 + - 𝐔 ) , 𝐏𝐯 + 〉 Γ - 〈 π ⁢ κ , 𝐯 - ⋅ 𝐧 〉 Γ .

Likewise, we find that the surface velocity and pressure satisfy the following integral identities:

(3.10) - ( π , div Γ ⁡ 𝐕 ) Γ + 2 ⁢ ( μ Γ ⁢ 𝐃 Γ ⁢ ( 𝐔 ) , 𝐃 Γ ⁢ ( 𝐕 ) ) Γ - 〈 f - ⁢ ( 𝐏𝐮 - - 𝐔 ) , 𝐕 〉 Γ - 〈 f + ⁢ ( 𝐏𝐮 + - 𝐔 ) , 𝐕 〉 Γ = ( 𝐏𝐛 e , 𝐕 ) Γ ,

(3.11) ( div Γ ⁡ 𝐔 , τ ) Γ = 0

for all ( 𝐕 , τ ) ∈ V Γ × L 0 2 ⁢ ( Γ ) .

The weak formulation of the coupled problem (3.1)–(3.7) follows by combining (3.8)–(3.9) and (3.10)–(3.11). In order to write it, we introduce the following forms for all 𝐮 ∈ V ± , 𝐯 ∈ 𝐕 0 ± , 𝐔 , 𝐕 ∈ V Γ , p ∈ L 2 ⁢ ( Ω ) , π ∈ L 2 ⁢ ( Γ ) :

a ⁢ ( { 𝐮 , 𝐔 } , { 𝐯 , 𝐕 } ) = 2 ⁢ ( μ - ⁢ 𝐃 ⁢ ( 𝐮 - ) , 𝐃 ⁢ ( 𝐯 - ) ) Ω - + 2 ⁢ ( μ + ⁢ 𝐃 ⁢ ( 𝐮 + ) , 𝐃 ⁢ ( 𝐯 + ) ) Ω + + 2 ⁢ ( μ Γ ⁢ 𝐃 Γ ⁢ ( 𝐔 ) , 𝐃 Γ ⁢ ( 𝐕 ) ) Γ + 〈 f - ⁢ ( 𝐏𝐮 - - 𝐔 ) , 𝐏𝐯 - - 𝐕 〉 Γ + 〈 f + ⁢ ( 𝐏𝐮 + - 𝐔 ) , 𝐏𝐯 + - 𝐕 〉 Γ ,

b ⁢ ( { 𝐯 , 𝐕 } , { p , π } ) = - ( p - , div ⁡ 𝐯 - ) Ω - - ( p + , div ⁡ 𝐯 + ) Ω + - ( π , div Γ ⁡ 𝐕 ) Γ ,

s ⁢ ( 𝐯 , π ) = - 〈 π ⁢ κ , 𝐯 ⋅ 𝐧 〉 Γ ,

r ⁢ ( 𝐯 , 𝐕 ) = ( 𝐟 - , 𝐯 - ) Ω - + ( 𝐟 + , 𝐯 + ) Ω + + ( 𝐏𝐛 e , 𝐕 ) Γ .

Then the weak formulation reads: find ( 𝐮 , p ) ∈ V ± × L ± 2 and ( 𝐔 , π ) ∈ V Γ × L 2 ⁢ ( Γ ) such that

(3.12) { a ⁢ ( { 𝐮 , 𝐔 } , { 𝐯 , 𝐕 } ) + b ⁢ ( { 𝐯 , 𝐕 } , { p , π } ) + s ⁢ ( 𝐯 , π ) = r ⁢ ( 𝐯 , 𝐕 ) b ⁢ ( { 𝐮 , 𝐔 } , { q , τ } ) = 0

for all ( 𝐯 , q ) ∈ V 0 ± × L 0 2 ⁢ ( Ω ) and ( 𝐕 , τ ) ∈ V Γ × L 0 2 ⁢ ( Γ ) . Note that test and trial pressure spaces both involve two (different) gauge conditions.

3.2 Well-Posedness

With the goal of proving the well-posedness of the stationary problem, we start by showing that a ⁢ ( { ⋅ , ⋅ } , { ⋅ , ⋅ } ) is coercive. Let ∥ 𝐮 ∥ H 1 ⁢ ( Ω ± ) 2 = ∥ 𝐮 + ∥ H 1 ⁢ ( Ω + ) 2 + ∥ 𝐮 - ∥ H 1 ⁢ ( Ω - ) 2 . Let ∥ p ∥ L 2 ⁢ ( Ω ± ) 2 = ∥ p ∥ Ω + 2 + ∥ p ∥ Ω - 2 . We define the following additional norms:

| | | { 𝐯 , 𝐕 } | | | 2 = ∥ 𝐯 ∥ H 1 ⁢ ( Ω ± ) 2 + ∥ 𝐕 ∥ H 1 ⁢ ( Γ ) 2 , | | | { p , π } | | | 2 = ∥ p ∥ Ω ± 2 + ∥ π ∥ Γ 2 .

The coercivity result is formulated in the form of a lemma.

Lemma 1.

For any v ∈ V 0 ± and V ∈ V Γ , it holds

(3.13) a ⁢ ( { 𝐯 , 𝐕 } , { 𝐯 , 𝐕 } ) ≥ C ⁢ | | | { 𝐯 , 𝐕 } | | | 2

with a positive constant C, which may depend on the viscosity values and Ω ± .

Proof.

One readily computes that

(3.14) a ⁢ ( { 𝐯 , 𝐕 } , { 𝐯 , 𝐕 } ) = 2 ⁢ ( μ - ⁢ 𝐃 ⁢ ( 𝐯 - ) , 𝐃 ⁢ ( 𝐯 - ) ) Ω - + 2 ⁢ ( μ + ⁢ 𝐃 ⁢ ( 𝐯 + ) , 𝐃 ⁢ ( 𝐯 + ) ) Ω + + 2 ⁢ ( μ Γ ⁢ 𝐃 Γ ⁢ ( 𝐕 ) , 𝐃 Γ ⁢ ( 𝐕 ) ) Γ + f - ⁢ ∥ 𝐏𝐯 - - 𝐕 ∥ Γ 2 + f + ⁢ ∥ 𝐏𝐯 + - 𝐕 ∥ Γ 2 .

Since function 𝐯 + satisfies homogeneous Dirichlet boundary condition on ∂ ⁡ Ω + ∖ Γ , we apply the following Korn inequality in Ω + :

(3.15) ∥ 𝐯 + ∥ H 1 ⁢ ( Ω + ) ≤ C ⁢ ∥ 𝐃 ⁢ ( 𝐯 + ) ∥ Ω + .

By the triangle and trace inequalities in Γ, we get

(3.16) ∥ 𝐕 ∥ Γ ≤ ∥ 𝐏𝐯 + - 𝐕 ∥ Γ + ∥ 𝐏𝐯 + ∥ Γ ≤ ∥ 𝐏𝐯 + - 𝐕 ∥ Γ + C ⁢ ∥ 𝐯 + ∥ H 1 ⁢ ( Ω + ) .

We further apply Korn’s inequality on Γ (see [30]),

(3.17) ∥ 𝐕 ∥ H 1 ⁢ ( Γ ) ≤ C ⁢ ( ∥ 𝐕 ∥ Γ + ∥ 𝐃 Γ ⁢ ( 𝐕 ) ∥ Γ ) .

Next, we can estimate the trace of 𝐯 - on Γ through the triangle inequality,

(3.18) ∥ 𝐯 - ∥ Γ ≤ ∥ 𝐏𝐯 - - 𝐕 ∥ Γ + ∥ 𝐕 ∥ Γ ≤ ∥ 𝐏𝐯 - - 𝐕 ∥ Γ + ∥ 𝐕 ∥ H 1 ⁢ ( Γ ) .

We finally apply the following Korn inequality in Ω - :

(3.19) ∥ 𝐯 ∥ H 1 ⁢ ( Ω - ) ≤ C ⁢ ( ∥ 𝐃 ⁢ ( 𝐯 - ) ∥ Ω - + ∥ 𝐯 - ∥ Γ ) .

Identity (3.14) and inequalities (3.15)–(3.19) lead to (3.13) after easy computations. ∎

The continuity of the bilinear forms a ⁢ ( { ⋅ , ⋅ } , { ⋅ , ⋅ } ) , b ⁢ ( { ⋅ , ⋅ } , { ⋅ , ⋅ } ) , and s ⁢ ( ⋅ , ⋅ ) follows from standard arguments based on the Cauchy–Schwarz and triangle inequalities,

(3.20)

a ⁢ ( { 𝐮 , 𝐔 } , { 𝐯 , 𝐕 } ) ≤ C ⁢ | | | { 𝐮 , 𝐔 } | | | ⁢ | | | { 𝐯 , 𝐕 } | | | for all ⁢ 𝐮 , 𝐯 ∈ V 0 ± , 𝐔 , 𝐕 ∈ V Γ ,

b ⁢ ( { 𝐮 , 𝐔 } , { p , π } ) ≤ C ⁢ | | | { 𝐮 , 𝐔 } | | | ⁢ | | | { p , π } | | | for all ⁢ 𝐮 ∈ V 0 ± , 𝐔 ∈ V Γ , p ∈ L 2 ⁢ ( Ω ) , π ∈ L 2 ⁢ ( Γ ) ,

s ⁢ ( 𝐯 , π ) ≤ C ⁢ | | | { 𝐯 , 0 } | | | ⁢ | | | { 0 , π } | | | for all ⁢ 𝐯 ∈ V ± , π ∈ L 2 ⁢ ( Γ ) .

Problem (3.12) falls into the class of so-called generalized saddle point problems. An abstract well-posedness result for such problems can be found, e.g., in [3, 43], which extend the Babuşka–Brezzi theory. Applied to (3.12), this well-posedness result requires coercivity (3.13), continuity (3.20) and two inf-sup conditions formulated in the following lemma.

Lemma 2.

The following inf-sup conditions hold with positive constants γ 1 and γ 2 :

(3.21) sup 𝐯 ∈ V 0 ± , 𝐕 ∈ V Γ ⁡ b ⁢ ( { 𝐯 , 𝐕 } , { p , π } ) + s ⁢ ( 𝐯 , π ) | | | { 𝐯 , 𝐕 } | | | ≥ γ 1 ⁢ | | | { p , π } | | | ⁢ for all ⁢ p ∈ L ± 2 , π ∈ L 2 ⁢ ( Γ ) ,

(3.22) sup 𝐯 ∈ V 0 ± , 𝐕 ∈ V Γ ⁡ b ⁢ ( { 𝐯 , 𝐕 } , { p , π } ) | | | { 𝐯 , 𝐕 } | | | ≥ γ 2 ⁢ | | | { p , π } | | | ⁢ for all ⁢ p ∈ L 0 2 ⁢ ( Ω ) , π ∈ L 0 2 ⁢ ( Γ ) .

Proof.

The proof follows by combining well-known results about the existence of a continuous right inverse of the divergence operator in H 0 1 ⁢ ( Ω ) 3 (see [5]) and V Γ (see [30]): for arbitrary p ∈ L 0 2 ⁢ ( Ω ) and π ∈ L 0 2 ⁢ ( Γ ) , there exist 𝐯 ∈ H 0 1 ⁢ ( Ω ) 3 and 𝐕 ∈ V Γ such that

(3.23)

p = div ⁡ 𝐯 ⁢ in ⁢ Ω and ∥ 𝐯 ∥ H 1 ⁢ ( Ω ) ≤ c Ω ∥ p ∥ L 2 ⁢ ( Ω ) ,

π = div Γ ⁡ 𝐕 ⁢ on ⁢ Γ and ∥ 𝐕 ∥ H 1 ⁢ ( Γ ) ≤ c Γ ∥ π ∥ Γ .

Letting 𝐯 ± = 𝐯 | Ω ± , ( 𝐯 - , 𝐯 + ) T ∈ V 0 ± , and adding estimates in (3.23), we get

(3.24) | | | { p , π } | | | 2 ≤ b ⁢ ( { 𝐯 , 𝐕 } , { p , π } ) ,

(3.25) | | | { 𝐯 , 𝐕 } | | | ≤ c Ω ⁢ ∥ p ∥ Ω + c Γ ⁢ ∥ π ∥ Γ ≤ ( c Ω + c Γ ) ⁢ | | | { p , π } | | | .

This proves (3.22) with γ 2 = 1 / ( c Ω + c Γ ) .

To show (3.21), we split π = π 0 + π ⊥ with π 0 ∈ L 0 2 ⁢ ( Γ ) and π ⊥ = | Γ | - 1 ⁢ ∫ Γ π ⁢ 𝑑 s . For the π 0 part of π, we use again (3.23) as above, while for p ± ∈ L 0 2 ⁢ ( Ω ± ) , we use the existence of a continuous right inverse of div in H 0 1 ⁢ ( Ω ± ) 3 to claim the existence of 𝐯 ∈ H 0 1 ⁢ ( Ω - ) 3 × H 0 1 ⁢ ( Ω + ) 3 ⊂ V 0 ± and 𝐕 ∈ V Γ such that

(3.26) | | | { p , π 0 } | | | ≤ b ⁢ ( { 𝐯 , 𝐕 } , { p , π 0 } ) + s ⁢ ( 𝐯 , π ) , | | | { 𝐯 , 𝐕 } | | | ≤ ( c Ω + c Γ ) ⁢ | | | { p , π 0 } | | | ,

with some positive c Ω , c Γ depending only on Γ and Ω. We also used that 𝐯 = 𝟎 on Γ implies s ⁢ ( 𝐯 , π ) = 0 .

Let C ± = ± | Ω ± | - 1 ⁢ | Γ | ⁢ ∫ Γ κ ⁢ 𝑑 s . To control ∥ π ⊥ ∥ Γ , we need 𝐯 1 ∈ V 0 ± such that

(3.27) div ⁡ 𝐯 1 = - C ± ⁢ in ⁢ Ω ± , 𝐯 1 ⋅ 𝐧 = κ ⁢ on ⁢ Γ and ∥ 𝐯 1 ∥ H 1 ⁢ ( Ω ± ) ≤ C .

Such 𝐯 1 can be built, for example, as follows. Let 𝐯 1 - = ∇ ⁡ ψ , where ψ ∈ H 2 ⁢ ( Ω - ) solves the Neumann problem - Δ ⁢ ψ = C - in Ω - , 𝐧 ⋅ ∇ ⁡ ψ = κ on Γ. Since Γ = ∂ ⁡ Ω - is smooth, by the H 2 -regularity of the Neumann problem, we have that ∥ 𝐯 1 - ∥ H 1 ⁢ ( Ω - ) ≤ ∥ ψ ∥ H 2 ⁢ ( Ω - ) ≤ C . The boundary ∂ ⁡ Ω is only Lipschitz, and so the Neumann problem in Ω + is not necessarily H 2 -regular. To handle this, we first extend 𝐯 1 - from Ω - to a function 𝐯 ~ 1 in H 0 1 ⁢ ( Ω ) 3 such that ∥ 𝐯 ~ 1 ∥ H 1 ⁢ ( Ω + ) ≤ c ⁢ ∥ 𝐯 1 - ∥ H 1 ⁢ ( Ω - ) (see [56]). Next, we consider 𝐰 ∈ H 0 1 ⁢ ( Ω + ) 3 such that div ⁡ 𝐰 = C + - div ⁡ 𝐯 ~ 1 ∈ L 0 2 ⁢ ( Ω + ) , and ∥ 𝐰 ∥ H 1 ⁢ ( Ω + ) ≤ c Ω + ⁢ ∥ div ⁡ 𝐰 ∥ L 2 ⁢ ( Ω + ) ≤ C (see [5]). The desired 𝐯 1 + is given in Ω + by 𝐯 1 + = 𝐯 ~ 1 + 𝐰 . Since div ⁡ 𝐯 1 = - C ± , for p ∈ L ± 2 ⁢ ( Ω ) and π ∈ L 2 ⁢ ( Γ ) , we have identities

(3.28) b ⁢ ( { 𝐯 1 , 0 } , { p , π } ) = 0 = ( div Γ ⁡ 𝐕 , π ⊥ ) Γ .

We also note the equality ∥ π ⊥ ∥ Γ 2 = c ^ ⁢ s ⁢ ( 𝐯 1 , π ⊥ ) , with

c ^ = π ⊥ ⁢ | Γ | / ∫ Γ κ 2 ⁢ 𝑑 s .

The denominator above is positive since Γ is closed, and so κ cannot be zero everywhere on Γ. We use (3.24)–(3.28) to estimate, for some β > 0 ,

| | | { p , π } | | | 2 = | | | { p , π 0 } | | | 2 + β ⁢ ∥ π ⊥ ∥ Γ 2 ≤ b ⁢ ( { 𝐯 , 𝐕 } , { p , π 0 } ) + s ⁢ ( 𝐯 , π ) + s ⁢ ( β ⁢ c ^ ⁢ 𝐯 1 , π ⊥ ) = b ⁢ ( { 𝐯 + β ⁢ c ^ ⁢ 𝐯 1 , 𝐕 } , { p , π } ) + s ⁢ ( 𝐯 + β ⁢ c ^ ⁢ 𝐯 1 , π ) - ( β ⁢ c ^ ⁢ 𝐯 1 , π 0 ) Γ ≤ b ⁢ ( { 𝐯 + β ⁢ c ^ ⁢ 𝐯 1 , 𝐕 } , { p , π } ) + s ⁢ ( 𝐯 + β ⁢ c ^ ⁢ 𝐯 1 , π ) + β 2 ⁢ c ^ 2 2 ⁢ ∥ 𝐯 1 ∥ Γ 2 + 1 2 ⁢ ∥ π 0 ∥ Γ 2 ≤ b ⁢ ( { 𝐯 + β ⁢ c ^ ⁢ 𝐯 1 , 𝐕 } , { p , π } ) + s ⁢ ( 𝐯 + β ⁢ c ^ ⁢ 𝐯 1 , π ) + c 3 ⁢ β 2 ⁢ ∥ π ⊥ ∥ Γ 2 + 1 2 ⁢ ∥ π 0 ∥ Γ 2 ,

with some c 3 > 0 depending only on Γ and Ω. For β > 0 sufficiently small such that β 2 - c 3 ⁢ β 2 ≥ 0 , we get

(3.29) c ⁢ | | | { p , π } | | | 2 ≤ b ⁢ ( { 𝐯 + β ⁢ c ^ ⁢ 𝐯 1 , 𝐕 } , { p , π } ) + s ⁢ ( 𝐯 + β ⁢ c ^ ⁢ 𝐯 1 , π ) ,

with c > 0 depending only on Γ and Ω. Thanks to the triangle inequality, the second estimate in (3.26), and the definition of c ^ and 𝐯 1 , we find the bound

| | | 𝐯 + β ⁢ c ^ ⁢ 𝐯 1 , 𝐕 | | | ≤ | | | 𝐯 , 𝐕 | | | + β ⁢ c ^ ⁢ | | | β ⁢ c ^ ⁢ 𝐯 1 , 0 | | | ≤ ( c Ω + c Γ ) ⁢ | | | { p , π 0 } | | | + C ⁢ ∥ π ⊥ ∥ Γ ≤ C ⁢ | | | { p , π } | | | ,

with C > 0 depending only on Γ and Ω. The combination of the above bound and (3.29) completes the proof of the lemma. ∎

3.3 Finite Element Discretization

Let Ω ⊂ ℝ 3 be a fixed polygonal domain that strictly contains Γ. We consider a family of shape regular tetrahedral triangulations { 𝒯 h } h > 0 of Ω. We adopt the convention that the elements T and edges e are open sets and use the over-line symbol to refer to their closure. Let h T = diam ⁡ ( T ) for T ∈ 𝒯 h . The set of elements intersecting Ω ± and the set of elements having a nonzero intersection with Γ are

𝒯 h ± = { T ∈ 𝒯 h : T ∩ Ω ± ≠ ∅ } , 𝒯 h Γ = { T ∈ 𝒯 h : T ¯ ∩ Γ ≠ ∅ } ,

respectively. We assume { 𝒯 h Γ } h > 0 to be quasi-uniform. The domain formed by all tetrahedra in 𝒯 h Γ is denoted by Ω h Γ : = int ( ⋃ T ∈ 𝒯 h Γ T ¯ ) . We define the h-dependent domains

Ω h ± = int ⁡ ( ⋃ T ∈ 𝒯 h ± T ¯ )

and the set of faces of 𝒯 h Γ restricted to the interior of Ω h ± ,

ℰ h Γ , ± = { e = int ⁡ ( ∂ ⁡ T 1 ∩ ∂ ⁡ T 2 ) : T 1 , T 2 ∈ 𝒯 h ± ⁢ and ⁢ T 1 ∩ Γ ≠ ∅ ⁢ or ⁢ T 2 ∩ Γ ≠ ∅ } .

For the space discretization of the bulk fluid problems, we restrict our attention to inf-sup stable finite element pair 𝐏 k + 1 - P k , k ≥ 1 , i.e. Taylor–Hood elements. Specifically, we consider the spaces of continuous finite element pressures given by

Q h - = { p ∈ C ⁢ ( Ω h - ) : q | T ∈ P k ⁢ ( T ) ⁢ for all ⁢ T ∈ 𝒯 h - } .

Space Q h + is defined analogously. The trial FE pressure space is given by:

L ± 2 ⁢ ( Ω ) h = { p = ( p - , p + ) ∈ Q h - × Q h + : ∫ Ω - p - = ∫ Ω + p + = 0 }

and the test space by Q h ± = Q h - × Q h + ∩ L 0 2 ⁢ ( Ω ) . Let

V h - = { 𝐮 ∈ C ⁢ ( Ω h - ) 3 : 𝐮 | T ∈ 𝐏 k + 1 ⁢ ( T ) ⁢ for all ⁢ T ∈ 𝒯 h - } ,

with the analogous definition for V h + . Our FE velocity space is given by

V h ± = { 𝐮 = ( 𝐮 - , 𝐮 + ) ∈ ( V h - × V h + ) } .

Functions in L ± 2 ⁢ ( Ω ) and V h ± and their derivatives are multivalued in Ω h Γ , the overlap of Ω h - and Ω h + . The jump of a multivalued function over the interface is defined as the difference of components coming from Ω h - and Ω h + , i.e. [ 𝐮 ] = 𝐮 - - 𝐮 + on Γ. Note that this is the jump that we have previously denoted with [ ⋅ ] + - . Moreover, we define the following averages:

(3.30) { 𝐮 } = α ⁢ 𝐮 + + β ⁢ 𝐮 - , 〈 𝐮 〉 = β ⁢ 𝐮 + + α ⁢ 𝐮 - ,

where α and β are weights to be chosen such that α + β = 1 , 0 ≤ α , β ≤ 1 . For example, in [12], the setting α = μ - / ( μ + + μ - ) and β = μ + / ( μ + + μ - ) is suggested. In [11], the authors choose α = 0 , β = 1 if μ - ≤ μ + and α = 1 , β = 0 otherwise. Below, in (3.34) and (3.39), we will use the relationship

(3.31) [ a ⁢ b ] = [ b ] ⁢ { a } + 〈 b 〉 ⁢ [ a ] .

For the discretization of the surface Stokes problem, we first consider the generalized Taylor–Hood bulk spaces in the strip Ω h Γ ,

V Γ , h = { 𝐔 ∈ C ⁢ ( Ω h Γ ) 3 : 𝐔 | T ∈ 𝐏 k + 1 ⁢ ( T ) ⁢ for all ⁢ T ∈ 𝒯 h Γ } ,

Q Γ , h = { π ∈ C ⁢ ( Ω h Γ ) : π | T ∈ P k ⁢ ( T ) ⁢ for all ⁢ T ∈ 𝒯 h Γ } ,

Q Γ , h 0 = Q Γ , h ∩ L 0 2 ⁢ ( Γ ) . In the trace finite element method, we use the traces of functions from V Γ , h and Q Γ , h on Γ. The inf-sup stability of the resulting trace FEM was analyzed in [48] for k = 1 and extended to higher order isoparametric trace elements in [31].

In the treatment of the surface Stokes problem, one has to enforce the tangentiality condition 𝐔 ⋅ 𝐧 = 0 on Γ. In order to enforce it while avoiding locking, we follow [24, 23, 30, 53, 45] and add a penalty term to the weak formulation.

A discrete variational analogue of problem (3.12) reads: find

( 𝐮 h , p h ) ∈ V h ± × L ± 2 ⁢ ( Ω ) h and ( 𝐔 h , π h ) ∈ V Γ , h × Q Γ , h

such that

(3.32) { a h ⁢ ( { 𝐮 h , 𝐔 h } , { 𝐯 h , 𝐕 h } ) + b h ⁢ ( { 𝐯 h , 𝐕 h } , { p h , π h } ) + s h ⁢ ( 𝐯 h , π h ) = r h ⁢ ( 𝐯 h , 𝐕 h ) , b h ⁢ ( { 𝐮 h , 𝐔 h } , { q h , τ h } ) - b p ⁢ ( p h , q h ) - b s ⁢ ( π h , τ h ) = 0

for all ( 𝐯 h , q h ) ∈ V 0 , h ± × Q h ± and ( 𝐕 h , τ h ) ∈ V Γ , h × Q Γ , h 0 . We define all the bilinear forms in (3.32) for all 𝐮 h ∈ V h ± , 𝐯 h ∈ V 0 , h ± , 𝐔 , 𝐕 ∈ V Γ , h , p ∈ L 2 ⁢ ( Ω ) , π ∈ L 2 ⁢ ( Γ ) . Let us start from form a h ⁢ ( { ⋅ , ⋅ } , { ⋅ , ⋅ } ) ,

(3.33) a h ⁢ ( { 𝐮 h , 𝐔 h } , { 𝐯 h , 𝐕 h } ) = a i ⁢ ( { 𝐮 h , 𝐔 h } , { 𝐯 h , 𝐕 h } ) + a n ⁢ ( 𝐮 h , 𝐯 h ) + a p ⁢ ( { 𝐮 h , 𝐔 h } , { 𝐯 h , 𝐕 h } ) + a s ⁢ ( 𝐔 h , 𝐕 h ) ,

where we group together the terms that arise from the integration by parts of the divergence of the stress tensors,

(3.34) a i ⁢ ( { 𝐮 h , 𝐔 h } , { 𝐯 h , 𝐕 h } ) = 2 ⁢ ( μ - ⁢ 𝐃 ⁢ ( 𝐮 h - ) , 𝐃 ⁢ ( 𝐯 h - ) ) Ω - + 2 ⁢ ( μ + ⁢ 𝐃 ⁢ ( 𝐮 h + ) , 𝐃 ⁢ ( 𝐯 h + ) ) Ω + + 〈 f - ⁢ ( 𝐏𝐮 h - - 𝐔 h ) , 𝐏𝐯 h - - 𝐕 h 〉 Γ + 〈 f + ⁢ ( 𝐏𝐮 h + - 𝐔 h ) , 𝐏𝐯 h + - 𝐕 h 〉 Γ - 2 ⁢ 〈 { μ ⁢ 𝐧 T ⁢ 𝐃 ⁢ ( 𝐮 h ) ⁢ 𝐧 } , [ 𝐯 h ⋅ 𝐧 ] 〉 Γ + 2 ⁢ ( μ Γ ⁢ 𝐃 Γ ⁢ ( 𝐔 h ) , 𝐃 Γ ⁢ ( 𝐕 h ) ) Γ ,

the terms that enforce condition (3.6) weakly using Nitsche’s method,

(3.35) a n ⁢ ( 𝐮 h , 𝐯 h ) = ∑ T ∈ 𝒯 h Γ γ h T ⁢ { μ } ⁢ ( [ 𝐮 h ⋅ 𝐧 ] , [ 𝐯 h ⋅ 𝐧 ] ) Γ - 2 ⁢ 〈 { μ ⁢ 𝐧 T ⁢ 𝐃 ⁢ ( 𝐯 ) ⁢ 𝐧 } , [ 𝐮 h ⋅ 𝐧 ] 〉 Γ ,

and the stabilization and penalty terms

(3.36) a p ⁢ ( { 𝐮 h , 𝐔 h } , { 𝐯 h , 𝐕 h } ) = 𝐉 h - ⁢ ( 𝐮 h , 𝐯 h ) + 𝐉 h + ⁢ ( 𝐮 h , 𝐯 h ) + τ s ⁢ ( 𝐔 h ⋅ 𝐧 , 𝐕 h ⋅ 𝐧 ) Γ ,

(3.37) 𝐉 h ± ⁢ ( 𝐮 h , 𝐯 h ) = ∑ ℓ = 1 k + 1 | e | 2 ⁢ ℓ - 1 ⁢ ∑ e ∈ ℰ h Γ , ± γ 𝐮 ± ⁢ μ ± ⁢ ( [ ∂ n ℓ ⁡ 𝐮 h ± ] , [ ∂ n ℓ ⁡ 𝐮 h ± ] ) e .

In (3.37), ∂ n ℓ ⁡ 𝐮 h - denotes the derivative of order ℓ of 𝐮 h - in the direction of 𝐧 . The 𝐉 h terms in (3.36) are so-called ghost-penalty stabilization [9, 10] included to avoid poorly conditioned algebraic systems due to possible small cuts of tetrahedra from 𝒯 h Γ by the interface. The terms in (3.38) and (3.40) have the same role for the surface bilinear forms.

The last form in (3.33) is related to the algebraic stability of the surface Stokes problem,

(3.38) a s ⁢ ( 𝐔 h , 𝐕 h ) = ρ 𝐮 ⁢ ( ∇ ⁡ 𝐮 h ⁢ 𝐧 , ∇ ⁡ 𝐯 h ⁢ 𝐧 ) Ω h Γ .

Similarly, the terms coming from the integration by parts of the divergence of the stress tensors are contained in

(3.39) b h ⁢ ( { 𝐯 h , 𝐕 h } , { p h , π h } ) = - ( p h - , div ⁡ 𝐯 h - ) Ω - - ( p h + , div ⁡ 𝐯 h + ) Ω + + 〈 { p h } , [ 𝐯 h ⋅ 𝐧 ] 〉 Γ + ( ∇ Γ ⁡ π h , 𝐕 h ) Γ ,

the penalty terms are grouped together in

b p ⁢ ( p h , q h ) = J h - ⁢ ( p h , q h ) + J h + ⁢ ( p h , q h ) , J h ± ⁢ ( p h , q h ) = γ p ± μ ± ⁢ ∑ e ∈ ℰ h Γ , ± ∑ ℓ = 1 k | e | 2 ⁢ ℓ + 1 ⁢ ( [ ∂ n ℓ ⁡ p h ± ] , [ ∂ n ℓ ⁡ q h ± ] ) e ,

and we have a term related to algebraic stability of the surface Stokes problem in

(3.40) b s ⁢ ( π h , τ h ) = ρ p ⁢ ( ∇ ⁡ p h ⋅ 𝐧 , ∇ ⁡ p h ⋅ 𝐧 ) Ω h Γ .

Finally,

s h ⁢ ( 𝐯 h , π h ) = - 〈 π h ⁢ κ , 〈 𝐯 h ⋅ 𝐧 〉〉 Γ ,

r h ⁢ ( 𝐯 h , 𝐕 h ) = ( 𝐟 h - , 𝐯 h - ) Ω - + ( 𝐟 h + , 𝐯 h + ) Ω + + ( 𝐏𝐛 h e , 𝐕 h ) Γ .

We recall that some of the interface terms in a i ⁢ ( { ⋅ , ⋅ } , { ⋅ , ⋅ } ) and b h ⁢ ( { ⋅ , ⋅ } , { ⋅ , ⋅ } ) have been obtained using relationship (3.31).

Parameters γ 𝐮 ± , γ p ± , and γ are all assumed to be independent of μ ± , h, and the position of Γ against the underlying mesh. Parameter γ in (3.35) needs to be large enough to provide the bilinear form a h ⁢ ( { ⋅ , ⋅ } , { ⋅ , ⋅ } ) with coercivity. Parameters γ 𝐮 ± and γ p ± can be tuned to improve the numerical performance of the method. As for the parameters required by the discretization of the surface Stokes problem, we allow

τ s = c τ ⁢ h - 2 , ρ p = c p ⁢ h , ρ 𝐮 ∈ [ c 𝐮 ⁢ h , C 𝐮 ⁢ h - 1 ] ,

where c τ , c p , c 𝐮 , and C 𝐮 are positive constants independent of h and how Γ cuts the bulk mesh.

The definition of bilinear forms requires integration over Γ ∩ T and T ∩ Ω ± for T from Ω h Γ . In general, there are no exact quadrature formulas to accomplish this task [49]. In practice, approximations should be made which introduce geometric errors. To keep these geometric errors of the order consistent with the approximation properties of the finite element spaces, we use isoparametric variants of the above spaces introduced in [35]; see also [37, 20].

We expect that the stability of the finite element formulation can be analyzed largely following the same steps of the well-posedness analysis for the weak formulation in Section 3.2, with a special treatment of cut elements, Nitsche terms, and surface elements as available in the literature for bulk Stokes interface and surface Stokes problems. For the sake of brevity, we do not work out these details here, but will present them in a follow-up paper.

4 A Partitioned Method for the Coupled Bulk-Surface Flow

For the solution of the coupled problem described in Section 3, we intend to use a partitioned strategy, i.e. each sub-problem is solved separately and the coupling conditions are enforced in an iterative fashion. Partitioned methods are appealing for solving coupled problems because they allow to reuse existing solvers with minimal modifications. In order to devise such a method for the simplified problem in Section 3, let us take a step back and look at the original problem (2.10).

Discretize problem (2.10) in time with, e.g., the backward Euler method, and consider the coupled problem at a particular time t = t n + 1 . Let S b be the map that associates the jump in the normal stress across the interface to any given surface flow velocity 𝐔 = 𝐔 T + U N ⁢ 𝐧 ,

S b ⁢ ( 𝐔 ) = [ 𝝈 ⁢ 𝐧 ] - + = 𝝈 + ⁢ ( 𝐮 + , p + ) ⁢ 𝐧 - 𝝈 - ⁢ ( 𝐮 - , p - ) ⁢ 𝐧 on ⁢ Γ ,

where ( 𝐮 + , p + ) and ( 𝐮 - , p - ) represent the solution of the two-phase time-discrete Navier–Stokes problem at time t associated to (2.1)–(2.2) endowed with interface conditions (2.5), (2.7), and (2.8). Moreover, let S s be the operator associated to the surface flow such that, to any given surface flow velocity 𝐔 , it associates the load 𝐟 Γ , S s ⁢ ( 𝐔 ) = 𝐟 Γ on Γ, through the time-discrete surface Navier–Stokes problem at time t associated to (2.3)–(2.4). Note that S b and S s are nonlinear and their definitions can involve also forcing terms and, in the case of the bulk fluid problem, terms due to the boundary conditions. For the surface operator, we can define S s - 1 as the map that associates the surface flow velocity 𝐔 to any given load 𝐟 Γ on Γ.

With the above definitions, we can express the time-discrete version of coupled problem (2.10) in terms of the solution 𝐔 of a nonlinear equation defined only on Γ. This interface equation is usually presented in one of three formulations that are equivalent from the mathematical point of view, but give rise to different iterative algorithms. The first and perhaps most used formulation is the fixed-point one: find 𝐔 such that

(4.1) S s - 1 ⁢ ( S b ⁢ ( 𝐔 ) ) = 𝐔 on ⁢ Γ .

The second formulation is a slight modification of (4.1), which lends itself to a Newton iterative method: find 𝐔 such that

S s - 1 ⁢ ( S b ⁢ ( 𝐔 ) ) - 𝐔 = 𝟎 on ⁢ Γ .

The third approach is given by the Steklov–Poincaré equation: find 𝐔 such that

S b ⁢ ( 𝐔 ) - S s ⁢ ( 𝐔 ) = 𝟎 on ⁢ Γ .

See, e.g., [50] for more details on these three formulations.

A standard algorithm for equation (4.1) uses fixed-point iterations: given 𝐔 k , compute

(4.2) 𝐔 k + 1 = 𝐔 k + ω k ⁢ ( 𝐔 ¯ k - 𝐔 k ) with ⁢ 𝐔 ¯ k = S s - 1 ⁢ ( S b ⁢ ( 𝐔 k ) ) .

The choice of the relaxation parameter ω k determines the efficiency of the algorithm or it might be crucial for convergence in certain ranges of the physical parameters. An effective strategy for setting ω k is Aitken’s acceleration method.

For simplicity, we present algorithm (4.2) applied to the time-discrete version of coupled problem (2.10) with ω k = 1 for all k (i.e., no relaxation). At time t = t n + 1 , assuming that 𝐔 k is known, perform the following steps.

Solve the two-phase time-discrete Navier–Stokes problem at time t associated to (2.1)–(2.2) for the bulk flow variables ( 𝐮 k + 1 - , p k + 1 - ) and ( 𝐮 k + 1 + , p k + 1 + ) with interface conditions

𝐮 k + 1 + ⋅ 𝐧 = U N k = 𝐮 k + 1 - ⋅ 𝐧 on ⁢ Γ ,

𝐏 ⁢ 𝝈 k + 1 + ⁢ 𝐧 = f + ⁢ ( 𝐏𝐮 k + 1 + - 𝐔 T k ) on ⁢ Γ ,

𝐏 ⁢ 𝝈 k + 1 - ⁢ 𝐧 = - f - ⁢ ( 𝐏𝐮 k + 1 - - 𝐔 T k ) on ⁢ Γ .
Solve the time-discrete surface Navier–Stokes problem at time t associated to (2.3)–(2.4) for variables ( 𝐔 k + 1 , π k + 1 ) with interface condition

𝐟 Γ k + 1 = [ 𝝈 k + 1 ⁢ 𝐧 ] - + on ⁢ Γ .
Check the stopping criterion

∥ 𝐔 k + 1 - 𝐔 k ∥ Γ < ϵ ⁢ ∥ 𝐔 k ∥ Γ ,

where ϵ is a given stopping tolerance.

Notice that the bulk and surface flow problems are solved separately and sequentially. In general, this algorithm is easy to implement, but convergence could be slow in certain ranges of the physical parameters and require relaxation for speed-up.

The above algorithm adapted to the simplified problem (3.1)–(3.7) reads as follows. At iteration k + 1 , assuming that ( 𝐔 T k , π k ) are known, perform the following steps.

Solve two-phase problem (3.1)–(3.2) for the bulk flow variables ( 𝐮 k + 1 - , p k + 1 - ) and ( 𝐮 k + 1 + , p k + 1 + ) with interface conditions

𝐮 k + 1 + ⋅ 𝐧 = 𝐮 k + 1 - ⋅ 𝐧 ⁢ on ⁢ Γ ,

(4.3) 𝐏 ⁢ 𝝈 k + 1 + ⁢ 𝐧 = f + ⁢ ( 𝐏𝐮 k + 1 + - 𝐔 T k ) ⁢ on ⁢ Γ ,

(4.4) 𝐏 ⁢ 𝝈 k + 1 - ⁢ 𝐧 = - f - ⁢ ( 𝐏𝐮 k + 1 - - 𝐔 T k ) ⁢ on ⁢ Γ ,

(4.5) = + - π k κ ⁢ on ⁢ Γ .
Solve surface flow problem (3.3)–(3.4) for variables ( 𝐔 T k + 1 , π k + 1 ) with interface condition

(4.6) 𝐏𝐟 Γ k + 1 = [ 𝐏 ⁢ 𝝈 k + 1 ⁢ 𝐧 ] - + on ⁢ Γ .
Check the stopping criterion

(4.7) ∥ 𝐔 T k + 1 - 𝐔 T k ∥ Γ < ϵ ⁢ ∥ 𝐔 T k ∥ Γ .

Notice that only interface conditions (4.3)–(4.6) are coupling conditions for bulk and surface flows. If one was to compute the load exerted on the surface fluid in (4.6) directly from the solution of the problem at Step 1, the overall accuracy of the method would be spoiled. Instead, one can compute 𝐏𝐟 Γ k + 1 by plugging (4.3)–(4.4) into (4.6),

𝐏𝐟 Γ k + 1 = f + ⁢ 𝐏𝐮 k + 1 + + f - ⁢ 𝐏𝐮 k + 1 - - ( f + + f - ) ⁢ 𝐔 T k on ⁢ Γ .

However, we prefer to use a more implicit version of the above condition,

𝐏𝐟 Γ k + 1 = f + ⁢ 𝐏𝐮 k + 1 + + f - ⁢ 𝐏𝐮 k + 1 - - ( f + + f - ) ⁢ 𝐔 T k + 1 on ⁢ Γ

since it could help have a better control of approximate rigid rotations (Killing vector fields).

5 Numerical Results

The aim of the numerical results collected in this section is to provide evidence of the robustness of the proposed finite element approach with respect to the contrast in viscosity in the bulk fluid, surface fluid viscosity, value of the slip coefficients, and position of the interface relative to the fixed computational mesh.

For the averages in (3.30), we set α = 0 and β = 1 for all the numerical experiments since we have μ - ≤ μ + . In addition, we set γ 𝐮 ± = 0.05 , γ p ± = 0.05 , and γ = 80 . The value of all other parameters will depend on the specific test. The stopping tolerance for criterion (4.7) is set to ϵ = 10 - 6 . For all the simulations, we choose to use finite element pair 𝐏 2 - P 1 for both the bulk and surface fluid problems.

For all the results presented below, we will report the L 2 error and a weighted H 1 error for the bulk velocity defined as

(5.1) ( 2 ⁢ μ - ⁢ ∥ D ⁢ ( 𝐮 - 𝐮 h - ) ∥ Ω - 2 + 2 ⁢ μ + ⁢ ∥ D ⁢ ( 𝐮 - 𝐮 h + ) ∥ Ω + 2 ) 1 2 ,

and a weighted L 2 error for the bulk pressure defined as

(5.2) ( μ - - 1 ⁢ ∥ p - p h - ∥ Ω - 2 + μ + - 1 ⁢ ∥ p - p h + ∥ Ω + 2 ) 1 2 .

Such weighted norms naturally arise in the error analysis of the Stokes interface problem [47]. In addition, we will report the L 2 and H 1 errors for the surface velocity and L 2 error for the surface pressure.

5.1 Sphere Embedded in a Cube

We perform a series of tests where domain Ω is the cube [ - 1.5 , 1.5 ] 3 and interface Γ is the unit sphere centered at the origin. Let 𝐱 = ( x , y , z ) ∈ Ω . Surface Γ is characterized as the zero-level set of function ϕ ⁢ ( 𝐱 ) = ∥ 𝐱 ∥ 2 2 - 1 . We consider the following solution for the bulk flow:

(5.3) p - = 3 ⁢ x ⁢ x 2 + y 2 + z 2 - 2 ⁢ x ⁢ ( x 2 + y 2 + z 2 ) , ⁢ 𝐮 - = 2 ⁢ f - f - - μ - ⁢ 𝐚 ⁢ ( x , y , z ) ,

(5.4) p + = 6 ⁢ x ⁢ x 2 + y 2 + z 2 - 4 ⁢ x ⁢ ( x 2 + y 2 + z 2 ) , ⁢ 𝐮 + = 2 ⁢ f + f + + μ + ⁢ 𝐚 ⁢ ( x , y , z ) ,

where

𝐚 ⁢ ( x , y , z ) = ( 3 2 - x 2 + y 2 + z 2 ) ⁢ [ ( - y - z ) ⁢ x + y 2 + z 2 ( - x - z ) ⁢ y + x 2 + z 2 ( - x - y ) ⁢ z + y 2 + x 2 ] ,

coupled to the following exact solution for the surface flow:

(5.5) π = x , 𝐔 = [ ( - y - z ) ⁢ x + y 2 + z 2 ( - x - z ) ⁢ y + x 2 + z 2 ( - x - y ) ⁢ z + y 2 + x 2 ] .

The forcing terms 𝐟 - and 𝐟 + are found by plugging solution (5.3)–(5.4) in (3.1). We impose a Dirichlet condition (2.11) on the faces x = 1.5 , y = - 1.5 , z = - 1.5 , where function 𝐠 is found from 𝐮 + in (5.4). On the remaining part of the boundary, we impose a Neumann condition (2.12), where 𝐟 N is found from p + in (5.3) and 𝐮 + in (5.4).

The value of the physical parameters will be specified for each test.

Spatial Convergence

To check the spatial accuracy of the finite element method described in Section 3.3, we consider exact solution (5.3)–(5.5) with viscosities μ - = 1 , μ + = 10 , and μ Γ = 1 , and friction coefficients f - = 2 and f + = 10 . Notice that the fluid outside the sphere has a larger viscosity than the fluid inside the sphere, which has the same viscosity as the surface fluid. We consider structured meshes of tetrahedra with five levels of refinement, the coarsest mesh having mesh size h = 0.5 , while the finest mesh has h = 0.05 . All the meshes feature a local one-level refinement near the corners of Ω. Table 1 reports the number of DOFs for each mesh. Figure 2 (left) shows the L 2 error and weighted H 1 error (5.1) for the bulk velocity, weighted L 2 error (5.2) for the bulk pressure, L 2 and H 1 errors for the surface velocity, and L 2 error for the surface pressure against the mesh size h. We observe optimal convergence rates for all the norms under consideration. Figure 2 (right) shows the number of bulk-surface iterations to satisfy stopping criterion (4.7) as h varies. As we can see, the number of iterations is fairly insensitive to a mesh refinement or coarsening.

Table 1

Sphere: DOFs for bulk and surface variables for all the meshes under consideration in the spatial convergence test.

h	0.5	0.25	0.125	0.0625	0.05
# bulk velocity DOFs	1.1e4	7.4e4	5.2e5	3.6e6	6.4e6
# bulk pressure DOFs	6.2e2	3.7e3	2.3e4	1.6e5	2.8e5
# surface velocity DOFs	2.4e3	1.0e4	4.0e4	1.5e5	2.2e5
# surface pressure DOFs	1.4e2	5.9e2	2.3e3	8.5e3	1.3e4

Figure 2

Sphere: (left) bulk and surface FE errors against the mesh size h. (right) Number of bulk-surface iterations of the partitioned method as h varies.

Robustness with Respect to the Viscosity Contrast

It is known that the case of high contrast for the viscosities in a two-phase problem is especially challenging from the numerical point of view. To test the robustness of our approach with respect to the viscosity contrast in the bulk, we consider exact solution (5.3)–(5.5) and fix μ - = 1 , while we let μ + vary from 1 to 256. We set μ Γ = 1 and friction coefficients f - = 2 and f + = 10 .

We consider one of the meshes adopted for the previous sets of simulations (with h = 0.125 ). Figure 3 (left) shows the L 2 error and weighted H 1 error (5.1) for the bulk velocity, weighted L 2 error (5.2) for the bulk pressure, L 2 and H 1 errors for the surface velocity, and L 2 error for the surface pressure against the value of μ + . We see that the errors remain mostly unchanged as μ + varies, with the exception of the weighted L 2 error for the bulk pressure, which decreases as μ + increases. In [46], which focuses only on two-phase bulk flow, we found that such error reaches a plateau as μ + is further increased. Figure 3 (left) shows that our approach is substantially robust with respect to the viscosity contrast μ + / μ - .

$Figure 3 Sphere: (left) bulk and surface FE errors against the value of μ + {\mu^{+}} . (right) Number of bulk-surface iterations of the partitioned method as μ + {\mu^{+}} varies.$

Figure 3

Sphere: (left) bulk and surface FE errors against the value of μ + . (right) Number of bulk-surface iterations of the partitioned method as μ + varies.

Figure 3 (right) reports the number of bulk-surface iterations to satisfy stopping criterion (4.7) as μ + varies. We observe that the number of iterations increases as the μ + / μ - ratio decreases, indicating that the coupled bulk-surface problem becomes more stiff as μ + decreases to match μ - and μ Γ .

Robustness with Respect to the Value of the Surface Viscosity

We now let μ Γ vary from 1 to 256 and keep all the other physical parameters fixed to the following values: μ - = 1 , μ + = 10 , f - = 2 , and f + = 10 . Again, we consider exact solution (5.3)–(5.5) and the mesh with mesh size h = 0.125 . Figure 4 (left) shows all the errors we have considered so far against the value of μ Γ . We notice that all the bulk errors stay constant as μ Γ varies. The L 2 errors for the surface velocity and pressure increase as μ Γ increases, while the H 1 error for the surface velocity slightly decreases as μ Γ increases. This experiment suggests that more viscous embedded layer is less controlled by the bulk fluid which effects the numerical stability of the complete system. In a water-lipid membrane system, the ratio of lateral dynamic viscosities of the embedded bi-layer and bulk water is typically 1–10 μm (depending on the temperature and composition) with the size of a vesicle being generally between 0.1 and 10 μm. Hence the observed increase of the numerical error does not look critical for this application.

$Figure 4 Sphere: (left) bulk and surface FE errors against the value of μ Γ {\mu_{\Gamma}} . (right) Number of bulk-surface iterations of the partitioned method as μ Γ {\mu_{\Gamma}} varies.$

Figure 4

Sphere: (left) bulk and surface FE errors against the value of μ Γ . (right) Number of bulk-surface iterations of the partitioned method as μ Γ varies.

Figure 4 (right) shows the number of bulk-surface iterations to satisfy stopping criterion (4.7) as μ Γ varies. Our partitioned method seems to be insensitive to a variation in the value of μ Γ . In particular, for the range of μ Γ under consideration the number of iterations stays constant and equal to 12.

Robustness with Respect to the Slip Coefficients

To check the sensitivity of the errors and partitioned method to the value of the slip coefficients, we run two sets of experiments, both involving exact solution (5.3)–(5.5). In the first set, we fix f + = 2 and vary f - from 1 to 256, while in the second set, we take f + = f - and let them both vary from 1 to 256. The viscosities are set as follows: μ - = 1 , μ + = 10 , and μ Γ = 1 . We consider again the mesh with mesh size h = 0.125 . Figures 5 (left) and 6 (left) show all the errors under consideration against the value of the slip coefficient(s) for both sets of tests. The only error that shows a substantial variation is the weighted H 1 error of the bulk velocity, which increases as the slip coefficient(s) increase. However, such error seems to reach a plateau in both cases.

$Figure 5 Sphere: (left) bulk and surface FE errors against the value of f + {f^{+}} . (right) Number of bulk-surface iterations of the partitioned method as f + {f^{+}} varies.$

Figure 5

Sphere: (left) bulk and surface FE errors against the value of f + . (right) Number of bulk-surface iterations of the partitioned method as f + varies.

Figures 5 (right) and 6 (right) report the number of bulk-surface iterations to satisfy stopping criterion (4.7) as the value of the coefficient(s) varies for both sets of tests. In Figure 5 (right), we see a rather sharp increase in the number of iterations as f - increases. This is even more true when both slip coefficients are increased together, as we can see from Figure 6 (right). Figure 7 reports the relative difference of the surface velocity between subsequent iterations in L 2 norm until stopping criterion (4.7) is met for f + = f - = 2 2 and f + = f - = 2 8 . We see that such relative difference decreases regularly for f + = f - = 2 2 , while for f + = f - = 2 8 , it decreases quickly for the first few iterations and then it slows down. A heuristic explanation we have for this is that, as the two friction coefficients increase, interface conditions (2.7)–(2.8) become close to Dirichlet conditions, making the surface flow more “passive”. Thus, separating the surface flow from the bulk flow as in the partitioned algorithm might not make much sense.

$Figure 6 Sphere: (left) bulk and surface FE errors against the value of f + = f - {f^{+}=f^{-}} . (right) Number of bulk-surface iterations of the partitioned method as the value of f + {f^{+}} and f - {f^{-}} (with f + = f - {f^{+}=f^{-}} ) varies.$

Figure 6

Sphere: (left) bulk and surface FE errors against the value of f + = f - . (right) Number of bulk-surface iterations of the partitioned method as the value of f + and f - (with f + = f - ) varies.

$Figure 7 Sphere: relative difference of the surface velocity between subsequent iterations in L 2 {L^{2}} norm until stopping criterion (4.7) is met.$

Figure 7

Sphere: relative difference of the surface velocity between subsequent iterations in L 2 norm until stopping criterion (4.7) is met.

5.2 Torus Embedded in a Cube

The domain Ω is the cube [ - 2 , 2 ] 2 and surface Γ is a torus centered at 𝐜 = ( c 1 , c 2 , c 3 ) . Let

( x , y ) = ( x ~ - c 1 , y ~ - c 2 , z ~ - c 3 ) , ( x ~ , y ~ , z ~ ) ∈ Ω .

We can characterize Γ as the zero-level set of function ϕ ⁢ ( 𝐱 ) = z 2 + ( x 2 + y 2 - 1 ) 2 - 1 2 . Finding an exact solution to problem (3.1)–(3.5), (2.7), and (2.8) with this more complicated surface is highly non-trivial. To simplify the task, we relax interface conditions (2.7), (2.8), and (3.5) as follows:

𝐏 ⁢ 𝝈 ± ⁢ 𝐧 = ± f ± ⁢ ( 𝐏𝐮 ± - 𝐔 ) + g ± on ⁢ Γ ,

[ 𝐧 T ⁢ 𝝈 ⁢ 𝐧 ] + - = π ⁢ κ + g n on ⁢ Γ ,

where g + , g - , and g n are computed such that the exact solution given below satisfies these relaxed interface conditions. The solution is given by

(5.6) p - = ( 1 2 - 2 - 4 ⁢ x 2 + y 2 x 2 + y 2 ) ⁢ ( x 3 + x ) , p + = 1 2 ⁢ ( x 3 + x ) , 𝐮 - = 𝐮 + = [ x 2 ⁢ y 5 - x ⁢ y 2 + z 2 - x ⁢ y ]

for the bulk and

(5.7) π = x 3 + x , 𝐔 = [ - z ⁢ x x 2 + y 2 , - z ⁢ y x 2 + y 2 , x 2 + y 2 - 1 ] T

for the surface. The forcing terms 𝐟 - and 𝐟 + are found by plugging solution (5.6)–(5.7) in (3.1). We impose a Dirichlet condition (2.11) on the faces x = 2 , y = - 2 , z = - 2 , where function 𝐠 is found from 𝐮 + in (5.6). On the remaining part of the boundary, we impose a Neumann condition (2.12), where 𝐟 n is found from p + and 𝐮 + in (5.6).

Spatial Convergence

Once again, we start by checking spatial accuracy. To this end, we consider exact solution (5.6)–(5.7) with 𝐜 = ( 0 , 0 , 0 ) , viscosities μ - = 1 , μ + = 10 , μ Γ = 1 , and friction coefficients f - = 2 and f + = 10 . Just like in the case of the sphere, we consider structured meshes of tetrahedra that feature a local one-level refinement near the corners of Ω. The details of the meshes under consideration are reported in Table 2. Figure 8 shows the L 2 error and weighted H 1 error (5.1) for the bulk velocity, weighted L 2 error (5.2) for the bulk pressure, L 2 and H 1 errors for the surface velocity, and L 2 error for the surface pressure against the mesh size h. Also for this second convergence test, we observe optimal convergence rates for all the norms.

Table 2

Torus: DOFs for bulk and surface variables for all the meshes under consideration in the spatial convergence test.

h	0.25	0.125	0.0625	0.05
# bulk velocity DOFs	1.6e5	1.2e6	8.5e6	1.5e7
# bulk pressure DOFs	7.6e3	5.4e4	3.7e5	6.7e5
# surface velocity DOFs	1.6e4	6.0e4	2.3e5	3.4e5
# surface pressure DOFs	9.0e2	3.4e3	1.3e4	2.0e4

Figure 8

Torus: bulk and surface FE errors against the mesh size h.

Robustness with Respect to the Position of the Interface

We conclude our series of numerical results with a set of simulations aimed at checking that our approach is not sensitive to the position of the interface with respect to the background mesh. We vary the center 𝐜 = ( c 1 , c 2 , c 3 ) of the torus that represents Γ,

(5.8) c 1 = h ⁢ k 20 ⁢ sin ⁡ ( k ⁢ π 10 ) , c 2 = h ⁢ k ⁢ 2 40 ⁢ cos ⁡ ( k ⁢ π 10 ) , c 3 = h ⁢ k ⁢ 2 40 ⁢ cos ⁡ ( k ⁢ π 10 ) ,

where h is the mesh size. The physical parameters are set like in the convergence test. We consider the mesh in Table 2 with h = 0.125 . Figure 9 shows all the errors against the value of k in (5.8). We see that all the errors are fairly insensitive to the position of Γ with respect to the background mesh, indicating robustness.

Figure 9

Torus: bulk and surface FE errors against the value of k in (5.8).

Funding source: National Science Foundation

Award Identifier / Grant number: DMS-1620384

Award Identifier / Grant number: DMS-1953535

Award Identifier / Grant number: DMS-2011444

Funding statement: This work has been partially supported by US National Science Foundation (NSF) through grants DMS-1620384, DMS-1953535, DMS-2011444.

Acknowledgements

We are grateful to Dr. Christoph Lehrenfeld for providing us with an ngsxfem implementation of isoparametric unfitted finite elements.

References

[1] G. J. Amador, D. van Dijk, R. Kieffer, M.-E. Aubin-Tam and D. Tam, Hydrodynamic shear dissipation and transmission in lipid bilayers, Proc. Natl. Acad. Sci. USA 118 (2021), 10.1073/pnas.2100156118. 10.1073/pnas.2100156118Search in Google Scholar PubMed PubMed Central

[2] J. W. Barrett, H. Garcke and R. Nürnberg, A stable numerical method for the dynamics of fluidic membranes, Numer. Math. 134 (2016), no. 4, 783–822. 10.1007/s00211-015-0787-5Search in Google Scholar PubMed PubMed Central

[3] C. Bernardi, C. Canuto and Y. Maday, Generalized inf-sup conditions for Chebyshev spectral approximation of the Stokes problem, SIAM J. Numer. Anal. 25 (1988), no. 6, 1237–1271. 10.1137/0725070Search in Google Scholar

[4] L. Bocquet and J.-L. Barrat, Flow boundary conditions from nano-to micro-scales, Soft Matter 3 (2007), 685–693. 10.1039/b616490kSearch in Google Scholar PubMed

[5] M. E. Bogovskiĭ, Solution of the first boundary value problem for an equation of continuity of an incompressible medium, Dokl. Akad. Nauk SSSR 248 (1979), no. 5, 1037–1040. Search in Google Scholar

[6] N. Bootland, A. Bentley, C. Kees and A. Wathen, Preconditioners for two-phase incompressible Navier–Stokes flow, SIAM J. Sci. Comput. 41 (2019), no. 4, B843–B869. 10.1137/17M1153674Search in Google Scholar

[7] S. P. A. Bordas, E. Burman, M. G. Larson and M. A. Olshanskii, Geometrically Unfitted Finite Element Methods and Applications, Lect. Notes Comput. Sci. Eng. 121. Springer, Cham, 2018. 10.1007/978-3-319-71431-8Search in Google Scholar

[8] D. Bothe and J. Prüss, On the two-phase Navier–Stokes equations with Boussinesq–Scriven surface fluid, J. Math. Fluid Mech. 12 (2010), no. 1, 133–150. 10.1007/s00021-008-0278-xSearch in Google Scholar

[9] E. Burman, Ghost penalty, C. R. Math. Acad. Sci. Paris 348 (2010), no. 21–22, 1217–1220. 10.1016/j.crma.2010.10.006Search in Google Scholar

[10] E. Burman, S. Claus, P. Hansbo, M. G. Larson and A. Massing, CutFEM: Discretizing geometry and partial differential equations, Internat. J. Numer. Methods Engrg. 104 (2015), no. 7, 472–501. 10.1002/nme.4823Search in Google Scholar

[11] E. Cáceres, J. Guzmán and M. Olshanskii, New stability estimates for an unfitted finite element method for two-phase Stokes problem, SIAM J. Numer. Anal. 58 (2020), no. 4, 2165–2192. 10.1137/19M1266897Search in Google Scholar

[12] S. Claus and P. Kerfriden, A CutFEM method for two-phase flow problems, Comput. Methods Appl. Mech. Engrg. 348 (2019), 185–206. 10.1016/j.cma.2019.01.009Search in Google Scholar

[13] M. Cooley and M. O’neill, On the slow motion generated in a viscous fluid by the approach of a sphere to a plane wall or stationary sphere, Mathematika 16 (1969), 37–49. 10.1112/S0025579300004599Search in Google Scholar

[14] G. Dziuk and C. M. Elliott, L 2 -estimates for the evolving surface finite element method, Math. Comp. 82 (2013), no. 281, 1–24. 10.1090/S0025-5718-2012-02601-9Search in Google Scholar

[15] T. Frachon and S. Zahedi, A cut finite element method for incompressible two-phase Navier–Stokes flows, J. Comput. Phys. 384 (2019), 77–98. 10.1016/j.jcp.2019.01.028Search in Google Scholar

[16] S. Ganesan, G. Matthies and L. Tobiska, On spurious velocities in incompressible flow problems with interfaces, Comput. Methods Appl. Mech. Engrg. 196 (2007), no. 7, 1193–1202. 10.1016/j.cma.2006.08.018Search in Google Scholar

[17] P. Gangl, K. Sturm, M. Neunteufel and J. Schöberl, Fully and semi-automated shape differentiation in NGSolve, Struct. Multidiscip. Optim. 63 (2021), no. 3, 1579–1607. 10.1007/s00158-020-02742-wSearch in Google Scholar PubMed PubMed Central

[18] R. B. Gennis, Biomembranes: Molecular Structure and Function, Springer, New York, 1989. 10.1007/978-1-4757-2065-5Search in Google Scholar

[19] D. Gérard-Varet, M. Hillairet and C. Wang, The influence of boundary conditions on the contact problem in a 3D Navier–Stokes flow, J. Math. Pures Appl. (9) 103 (2015), no. 1, 1–38. 10.1016/j.matpur.2014.03.005Search in Google Scholar

[20] J. Grande, C. Lehrenfeld and A. Reusken, Analysis of a high-order trace finite element method for PDEs on level set surfaces, SIAM J. Numer. Anal. 56 (2018), no. 1, 228–255. 10.1137/16M1102203Search in Google Scholar

[21] S. Gross and A. Reusken, Numerical Methods for Two-Phase Incompressible Flows, Springer Ser. Comput. Math. 40, Springer, Berlin, 2011. 10.1007/978-3-642-19686-7Search in Google Scholar

[22] M. E. Gurtin and A. I. Murdoch, A continuum theory of elastic material surfaces, Arch. Ration. Mech. Anal. 57 (1975), 291–323. 10.1007/BF00261375Search in Google Scholar

[23] P. Hansbo, M. G. Larson and K. Larsson, Analysis of finite element methods for vector Laplacians on surfaces, IMA J. Numer. Anal. 40 (2020), no. 3, 1652–1701. 10.1093/imanum/drz018Search in Google Scholar

[24] P. Hansbo, M. G. Larson and A. Massing, A stabilized cut finite element method for the Darcy problem on surfaces, Comput. Methods Appl. Mech. Engrg. 326 (2017), 298–318. 10.1016/j.cma.2017.08.007Search in Google Scholar

[25] P. Hansbo, M. G. Larson and S. Zahedi, A cut finite element method for a Stokes interface problem, Appl. Numer. Math. 85 (2014), 90–114. 10.1016/j.apnum.2014.06.009Search in Google Scholar

[26] X. He, F. Song and W. Deng, Stabilized nonconforming Nitsche’s extended finite element method for Stokes interface problems, preprint (2019), https://arxiv.org/abs/1905.04844. 10.3934/dcdsb.2021163Search in Google Scholar

[27] W. Helfrich, Elastic properties of lipid bilayers: Theory and possible experiments, Z. Naturforschung 28 (1973), 693–703. 10.1515/znc-1973-11-1209Search in Google Scholar PubMed

[28] L. M. Hocking, The effect of slip on the motion of a sphere close to a wall and of two adjacent spheres, J. Engrg. Math. 7 (1973), 207–221. 10.1007/BF01535282Search in Google Scholar

[29] M. Hömberg and M. Müller, The role of inertia and coarse-graining on the transverse modes of lipid bilayers, EPL 97 (2012), Article ID 68010. 10.1209/0295-5075/97/68010Search in Google Scholar

[30] T. Jankuhn, M. A. Olshanskii and A. Reusken, Incompressible fluid problems on embedded surfaces: modeling and variational formulations, Interfaces Free Bound. 20 (2018), no. 3, 353–377. 10.4171/IFB/405Search in Google Scholar

[31] T. Jankuhn, M. A. Olshanskii, A. Reusken and A. Zhiliakov, Error analysis of higher order trace finite element methods for the surface Stokes equation, J. Numer. Math. 29 (2021), no. 3, 245–267. 10.1515/jnma-2020-0017Search in Google Scholar

[32] V. John, Slip with friction and penetration with resistance boundary conditions for the Navier–Stokes equations—numerical tests and aspects of the implementation, J. Comput. Appl. Math. 147 (2002), no. 2, 287–300. 10.1016/S0377-0427(02)00437-5Search in Google Scholar

[33] K. Kawano, E. Onose, Y. Hattori and Y. Maitani, Higher liposomal membrane fluidity enhances the in vitro antitumor activity of folate-targeted liposomal mitoxantrone, Molecular Pharmaceutics 6 (2009), 98–104. 10.1021/mp800069cSearch in Google Scholar

[34] E. Lauga, M. Brenner and H. Stone, Microfluidics: The no-slip boundary condition, Springer Handbooks, Springer, Berlin (2007), 1219–1240. 10.1007/978-3-540-30299-5_19Search in Google Scholar

[35] C. Lehrenfeld, High order unfitted finite element methods on level set domains using isoparametric mappings, Comput. Methods Appl. Mech. Engrg. 300 (2016), 716–733. 10.1016/j.cma.2015.12.005Search in Google Scholar

[36] C. Lehrenfeld, A higher order isoparametric fictitious domain method for level set domains, Geometrically Unfitted Finite Element Methods and Applications, Springer, Cham (2017), 65–92. 10.1007/978-3-319-71431-8_3Search in Google Scholar

[37] C. Lehrenfeld and A. Reusken, Analysis of a high-order unfitted finite element method for elliptic interface problems, IMA J. Numer. Anal. 38 (2018), no. 3, 1351–1387. 10.1093/imanum/drx041Search in Google Scholar

[38] A. Massing, M. G. Larson, A. Logg and M. E. Rognes, A stabilized Nitsche overlapping mesh method for the Stokes problem, Numer. Math. 128 (2014), no. 1, 73–101. 10.1007/s00211-013-0603-zSearch in Google Scholar

[39] T.-H. Miura, On singular limit equations for incompressible fluids in moving thin domains, Quart. Appl. Math. 76 (2018), no. 2, 215–251. 10.1090/qam/1495Search in Google Scholar

[40] P. B. Moore, C. F. Lopez and M. L. Klein, Dynamical properties of a hydrated lipid bilayer from a multinanosecond molecular dynamics simulation, Biophys. J. 81 (2001), 2484–2494. 10.1016/S0006-3495(01)75894-8Search in Google Scholar

[41] C. Morris and U. Homann, Cell surface area regulation and membrane tension, J. Membrane Biol. 179 (2001), Paper No. 79. 10.1007/s002320010040Search in Google Scholar PubMed

[42] C. Navier, Mémoire sur les lois du mouvement des fluides, Mém. Acad. Roy. Sci. Inst. France 6 (1823), 389–440. Search in Google Scholar

[43] R. A. Nicolaides, Existence, uniqueness and approximation for generalized saddle point problems, SIAM J. Numer. Anal. 19 (1982), no. 2, 349–357. 10.1137/0719021Search in Google Scholar

[44] I. Nitschke, A. Voigt and J. Wensch, A finite element approach to incompressible two-phase flow on manifolds, J. Fluid Mech. 708 (2012), 418–438. 10.1017/jfm.2012.317Search in Google Scholar

[45] M. A. Olshanskii, A. Quaini, A. Reusken and V. Yushutin, A finite element method for the surface Stokes problem, SIAM J. Sci. Comput. 40 (2018), no. 4, A2492–A2518. 10.1137/18M1166183Search in Google Scholar

[46] M. A. Olshanskii, A. Quaini and Q. Sun, An unfitted finite element method for two-phase Stokes problems with slip between phases, preprint (2021), https://arxiv.org/abs/2101.09627. 10.1007/s10915-021-01658-xSearch in Google Scholar

[47] M. A. Olshanskii and A. Reusken, Analysis of a Stokes interface problem, Numer. Math. 103 (2006), no. 1, 129–149. 10.1007/s00211-005-0646-xSearch in Google Scholar

[48] M. A. Olshanskii, A. Reusken and A. Zhiliakov, Inf-sup stability of the trace 𝐏 2 – P 1 Taylor–Hood elements for surface PDEs, Math. Comp. 90 (2021), no. 330, 1527–1555. 10.1090/mcom/3551Search in Google Scholar

[49] M. A. Olshanskii and D. Safin, Numerical integration over implicitly defined domains for higher order unfitted finite element methods, Lobachevskii J. Math. 37 (2016), no. 5, 582–596. 10.1134/S1995080216050103Search in Google Scholar

[50] A. Quarteroni and A. Valli, Domain Decomposition Methods for Partial Differential Equations, Numer. Math. Sci. Comput., The Clarendon, New York, 1999. 10.1093/oso/9780198501787.001.0001Search in Google Scholar

[51] A. Reusken and Y. Zhang, Numerical simulation of incompressible two-phase flows with a Boussinesq–Scriven interface stress tensor, Internat. J. Numer. Methods Fluids 73 (2013), no. 12, 1042–1058. 10.1002/fld.3835Search in Google Scholar

[52] S. Reuther and A. Voigt, The interplay of curvature and vortices in flow on curved surfaces, Multiscale Model. Simul. 13 (2015), no. 2, 632–643. 10.1137/140971798Search in Google Scholar

[53] S. Reuther and A. Voigt, Solving the incompressible surface Navier–Stokes equation by surface finite elements, Phys. Fluids 30 (2018), Article ID 012107. 10.1063/1.5005142Search in Google Scholar

[54] D. S. Rodrigues, R. F. Ausas, F. Mut and G. C. Buscaglia, A semi-implicit finite element method for viscous lipid membranes, J. Comput. Phys. 298 (2015), 565–584. 10.1016/j.jcp.2015.06.010Search in Google Scholar

[55] G. Salbreux and F. Jülicher, Mechanics of active surfaces, Phys. Rev. E 96 (2017), Article ID 032404. 10.1103/PhysRevE.96.032404Search in Google Scholar PubMed

[56] E. M. Stein, Singular Integrals and Differentiability Properties of Functions, Princeton Math. Ser. 30, Princeton University, Princeton, 1970. 10.1515/9781400883882Search in Google Scholar

[57] A. Torres-Sánchez, D. Millán and M. Arroyo, Modelling fluid deformable surfaces with an emphasis on biological interfaces, J. Fluid Mech. 872 (2019), 218–271. 10.1017/jfm.2019.341Search in Google Scholar

[58] N. Wang and J. Chen, A nonconforming Nitsche’s extended finite element method for Stokes interface problems, J. Sci. Comput. 81 (2019), no. 1, 342–374. 10.1007/s10915-019-01019-9Search in Google Scholar

[59] Y. Wang, A. Zhiliakov, A. Quaini, M. Olshanskii and S. Majd, Lipid domain formation and dynamics in multicomponent membranes: Experimental validation of a phase-field model, Biophys. J. 120 (2021), Article ID 225. 10.1016/j.bpj.2020.11.1503Search in Google Scholar

[60] A. Yavari, A. Ozakin and S. Sadik, Nonlinear elasticity in a deforming ambient space, J. Nonlinear Sci. 26 (2016), no. 6, 1651–1692. 10.1007/s00332-016-9315-8Search in Google Scholar

[61] V. Yushutin, A. Quaini, S. Majd and M. Olshanskii, A computational study of lateral phase separation in biological membranes, Int. J. Numer. Methods Biomed. Eng. 35 (2019), no. 3, Article ID e3181. 10.1002/cnm.3181Search in Google Scholar PubMed

[62] Netgen/NGSolve, https://ngsolve.org/. Search in Google Scholar

[63] ngsxfem, https://github.com/ngsxfem/ngsxfem/tree/49205a1ae637771a0ed56d4993ce99008f3a00e0. Search in Google Scholar

Received: 2021-09-29

Accepted: 2021-10-03

Published Online: 2021-11-03

Published in Print: 2022-04-01

Articles in the same Issue

https://doi.org/10.1515/cmam-2021-0185

Keywords for this article

Two-Phase Flow; Surface Stokes Equations; Finite Elements; CutFEM; TraceFEM