{'fraction_non_alphanumeric': 0.041275947743848684, 'fraction_numerical': 0.013785365071731987, 'mean_word_length': 4.140061791967044, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A cut finite element method for coupled bulk-surface problems on time-dependent domains

In this contribution we present a new computational method for coupled bulk-surface problems on time-dependent domains. The method is based on a space-time formulation using discontinuous piecewise linear elements in time and continuous piecewise linear elements in space on a fixed background mesh. The domain is represented using a piecewise linear level set function on the background mesh and a cut finite element method is used to discretize the bulk and surface problems. In the cut finite element method the bilinear forms associated with the weak formulation of the problem are directly evaluated on the bulk domain and the surface defined by the level set, essentially using the restrictions of the piecewise linear functions to the computational domain. In addition a stabilization term is added to stabilize convection as well as the resulting algebraic system that is solved in each time step. We show in numerical examples that the resulting method is accurate and stable and results in well conditioned algebraic systems independent of the position of the interface relative to the background mesh.


Introduction
Problems involving phenomena that take place both on surfaces (or interfaces) and in bulk domains occur in a variety of applications in fluid dynamics and biology. In this paper, we consider two-phase flow problems with soluble surfactants. A soluble surfactant is dissolved in the bulk fluid but also exists in adsorbed form on the interface separating two immiscible fluids. Surfactants have a large influence on the dynamics in multiphase flow systems in that they may cause drop-breakup or coalescence due to their ability to reduce the surface tension. They were for example used to lower the surface tension of oil droplets in the 2010 Deepwater Horizon oil spill so that the oil became more soluble in the water. Other examples of applications where the effects of surfactants are important include drug delivery, treatment of lung diseases, and polymer blending [1]. Simulations can be a powerful tool to better understand such complex systems.
Soluble surfactants are modeled by a coupled system of time-dependent convection-diffusion equations describing the concentration of surfactants in the bulk fluid and on the interface. The interface is moving with the fluid velocity. The equations governing the fluid flow are typically the time-dependent Stokes or the Navier-Stokes equations in the subdomains, each occupied by a fluid, with interface conditions due to the surface tension force. The effect of surfactants on the flow dynamics is modeled through the surface tension. From a computational point of view, the main challenge is that the differential equations are defined on domains that are evolving with time and that these domains may undergo strong deformations.
A common strategy is to let the mesh conform to the time-dependent domain, see, e.g., [2]. This technique can be made accurate but requires remeshing and interpolation as an interface evolves with time and leads to significant complications when topological changes such as dropbreakup or coalescence occur, especially in three space dimensions. Therefore, computational methods that allow the interface to be arbitrarily located with respect to a fixed background mesh, so called fixed grid methods, are preferable and have become highly attractive and significant efforts have been directed to their development, see, e.g., [3,4,5,6]. In general, strategies for solving quantities such as the surfactant concentration on evolving domains are developed on the basis of the interface representation technique. Techniques to represent the interface can be roughly divided into two classes: explicit representation, e.g., by marker particles [7], and implicit representation, e.g., by the level set of a higher dimensional function [8]. In consequence, existing fixed grid methods are usually tightly coupled to the interface representation. Several methods have been proposed both based on explicit [9,10,11,12,13,14] and implicit [15,16,17,18,19] representations. Most work has been done on insoluble surfactants, i.e., surfactants that are only present at the interface without surfactant mass transfer between the interface and the bulk. In most existing methods using implicit representation techniques, the surfactant concentration is extended to a region embedding the interface, and instead of a surface Partial Differential Equation (PDE), a PDE on a higher dimensional domain must be solved for the interfacial surfactant.
In this paper, we present a new computational method for the coupled bulk-surface problem on time-dependent domains modeling the evolution of soluble surfactants. The method can be used both with explicit and implicit interface representation techniques. The surface PDE is solved on the interface which can be arbitrarily located with respect to the fixed background mesh. The method is accurate and stable and results in well conditioned linear systems independently of how the interface cuts through the background mesh, and the total mass of surfactants is accurately conserved.
Our strategy is to embed the time-dependent domain where the PDE has to be solved in a fixed background grid, equipped with a standard finite element space, and then take the restriction of the finite element functions to the time-dependent domain. This idea was first proposed for an elliptic problem with a stationary fictitious boundary in [6] and for the Laplace-Beltrami operator on a stationary interface in [20]. It has been extended to other equations with error analysis for example the Stokes equations involving two immiscible incompressible fluids with different viscosities and with surface tension [21] and to PDE:s on time-dependent surfaces in, e.g., [22,23]. We call this type of methods cut finite element methods (CutFEM) since the interface cuts through the background grid in an arbitrary fashion.
In this paper, we present an optimal order accurate cut finite element method based on a space-time approach with continuous linear elements in space and discontinuous piecewise linear elements in time. To ensure well posedness of the resulting algebraic systems of equations, independent of the position of the interface in the background grid and also to ensure a stable discretization in case of convection dominated problems, we add a consistent stabilization term [24,25]. In addition, the total mass of surfactants is accurately conserved using a Lagrange multiplier. Numerical results indicate that the method is optimal order accurate. We have proven optimal order of accuracy for a stationary coupled bulk-surface problem with a linear coupling term in [26].
The proposed finite element method is very convenient for multiphase flow simulations since the same finite element spaces can be used for solving problems both in the bulk domain and on the interfaces separating immiscible fluids. In this paper, we have used the standard level set method to represent the interface but the method can be used with other interface representation techniques as well. We have chosen to concentrate on the challenging task of solving the coupled bulk-surface problem on time-dependent domains and will throughout the paper assume that the velocity field is given. However, the method can be coupled to a flow solver solving for example the Stokes or the Navier-Stokes equations.
The remainder of the paper is outlined as follows. In Section 2 we formulate the transport equations for soluble surfactants. In Section 3 we introduce a discrete approximation of the interface and state our assumptions on the geometry. The computational method for the coupled bulk-surface problem modeling soluble surfactants is presented in Section 4. In Section 5, we show numerical examples. Finally, in Section 6 we summarize our results.

Transport equations for soluble surfactants
Let Ω be an open bounded domain in R d , d = 2, 3, with convex polygonal boundary ∂Ω and let I = [0, T ] be a time interval. We consider two immiscible incompressible fluids that occupy time dependent subdomains Ω i (t) ⊂ Ω, i = 1, 2, with t ∈ I, such that Ω = Ω 1 (t) ∪ Ω 2 (t) and Ω 1 (t) ∩ Ω 2 (t) = ∅, and are separated by a smooth interface defined by Γ(t) = ∂Ω 1 (t) ∩ ∂Ω 2 (t). See Fig. 1 for an illustration in two dimensions. We assume that Γ(t) is a simple closed curve (2D) or surface (3D) moving with a given divergence free velocity field β : I × Ω → R d in such a way that it does not intersect the boundary ∂Ω (Γ ∩ ∂Ω = ∅) or itself for any t ∈ I. In applications the velocity field β is typically obtained from the incompressible Stokes or the incompressible Navier-Stokes equations. For simplicity we further assume that the surfactant is soluble only in the outer fluid phase Ω 1 (t). The model for soluble surfactants is given by a time-dependent convection-diffusion equation on the interface coupled with a time-dependent convection-diffusion equation in the bulk. The concentration of surfactants in the bulk and the concentration of surfactants on the surface are coupled through a nonlinear term which enters as a source term for the interfacial surfactant and in a Neumann boundary condition for the bulk concentration.
More precisely, we consider the following time dependent coupled bulk-surface problem: find u B : I × Ω 1 → R and u S : with initial condition u B (0) = u 0 B and u S (0) = u 0 S on Ω 1 (0) and Γ(0). Here ∂ t = ∂ ∂t , ∇ is the usual R d gradient, ∇ Γ is the tangent gradient associated with Γ defined by n is the unit normal vector to Γ, outward-directed with respect to Ω 1 , n ∂Ω is the outward directed unit normal vector on ∂Ω, I is the identity matrix, ⊗ denotes outer product ((a ⊗ b) ij = a i b j for any two vectors a and b), and k B and k S are the bulk diffusion and the interfacial diffusion coefficients, respectively.
The exchange of surfactants between the interface and the bulk is modeled with the term f coupling . We consider, in particular, the Langmuir model where Here u ∞ S is the maximum surfactant concentration on Γ and k a and k d are adsorption and desorption coefficients, respectively. See [27] for other models for the coupling term f coupling . Using the fact that ∇ · β = 0 we also have the conservation law for the total amount of surfactants on the surface and in the bulk.

Non-dimensional form
In the following we use the same non-dimensional variables as in [2]. Let L, β ∞ , u ∞ S , and u ∞ B be the characteristic values for length, velocity, and surface and bulk surfactant concentration.
The non-dimensional form of the surfactant concentration equations is given by Here the non-dimensional numbers Pe and Pe s are the bulk and surface (interfacial) Peclet numbers, defined with respect to the bulk and surface diffusivity k B and k S , Da is the Damköhler number and Bi is the Biot number Due to the non-dimensionalization, the conserved total amount of surfactants is

Weak form
First note that we can write f coupling in the form In case f coupling is given by equation (2.7) we have b B = k a u ∞ S , b S = k d , and b BS = k a while in the non-dimensional form of the equations f coupling is given by equation (2.13) and b B = b BS = α and b S = Bi. We assume that b B , b S , and b BS are positive constants.

)). After integration by parts and incorporating the boundary conditions (2.2) -(2.3) we obtain the variational problem: find
Note that a BS (t, u, v) is the linear part of the coupling term f coupling . For the non-dimensional form of the equations we can multiply the bulk PDE (2.9) with a test function b B Da v B ∈ H 1 (Ω 1 (t)) and the surface PDE (2.12) as before with a test function b S v S ∈ H 1 (Γ(t)). This yields after integration by parts and incorporating the boundary conditions the following weak form and a B , a S , and a BS as before. Note that there is only a change in the coefficient in front of a B and (∂ t u B , v B ) Ω 1 (t) (terms coming from the bulk PDE).

The level set representation of the interface
For t ∈ I, let U(Γ(t)) be an open neighborhood of Γ(t) such that for each x ∈ U(Γ(t)) there is a uniquely determined closest point in Γ(t). We let ρ(t) : R d → R be the signed distance function. The exterior unit normal n = n(t, x) is the spatial gradient of the signed distance function, Given a vector field β the evolution of the surface Γ(t) is governed by the following problem for the level set function: find ρ : Let K 0,h be a quasiuniform partition of Ω into shape regular triangles for d = 2 and tetrahedra for d = 3 of diameter h. We approximate the level set function ρ by ρ h ∈ V 0,h/2 where V 0,h/2 is the space of piecewise linear continuous functions defined on the mesh K 0,h/2 obtained by refining K 0,h uniformly once. Γ h (t) is the zero level set of ρ h (t). We consider a continuous piecewise linear approximation Γ h of Γ such that Γ h ∩ K is a linear segment for d = 2 and is a subset of a hyperplane in R 3 , for each K ∈ K 0,h/2 . We assume that for every t ∈ I, Γ h (t) ⊂ U(Γ(t)) and that the following approximation assumptions hold: and for all t ∈ I. Here denotes less or equal up to a positive constant, n h denotes the piecewise constant exterior unit normal to Γ h with respect to Ω 1 and n e is the extension of the exact normal to Γ h by the closest point mapping. These assumptions are consistent with the piecewise linear nature of the discrete interface. We define Ω h,1 as the domain enclosed by Γ h ∪ ∂Ω.

The space-time cut finite element method
In a space-time formulation the variational formulation is written over the space-time domain which is divided into space-time slabs corresponding to each time step. We employ piecewise linear trial and test functions that are continuous in space but allowed to be discontinuous from one space-time slab to another. For a given time t n (where n is the time step number) we thus have two distinct solutions, at times t ± n := lim →0 t n ± , see, e.g. [28]. Using a suitable weak enforcement of continuity at t n the discrete equations can then be solved one space-time slab at a time. The discontinuous Galerkin method in time can be compared with implicit finite difference methods and good stability is expected (the method is related to the fisct subdiagonal Padé approximation, cf. Thomée [29]). In space we will use a CutFEM, which means that we will use restrictions of standard continuous finite element functions defined on the background grid to our time dependent domain. Since the interface can cut through the fixed background grid arbitrarily there is a lack of shape regularity which results in ill-conditioned system matrices.
To avoid this problem, we add stabilization terms similar to [26]. The same stabilization terms also stabilize the proposed finite element method in case of convection dominated problems.

The Mesh and Finite Element Spaces
Let 0 = t 0 < t 1 < · · · < t N = T be a partition of the time interval I = [0, T ] into time steps I n = (t n−1 , t n ] of length k n = t n − t n−1 for n = 1, 2, . . . , N . Let P 1 (I n ) be the space of We define the following sets of elements and the corresponding sets For an illustration of the sets in two dimensions see Fig. 2. We let V 0,h be the space of piecewise linear continuous functions defined on K 0,h . On each space-time slab S n = (N n B,h × I n ) × (N n B,h × I n ) we define the spaces where t ∈ I n and v B,j and v S,j , j = 0, 1 can be written as Here ξ B ij , ξ S ij ∈ R are coefficients, ϕ i (x) is the standard nodal basis function associated with mesh vertex i, N B and N S are the number of nodes in N n B,h and in N n S,h respectively.

Weak form
Given u h (x, t − n−1 ) and u 0 (see equation (2.8)) find u h = (u B , u S ) ∈ W n h and λ ∈ R, such that where j h (u, v) is a stabilizing term of the form τ B , τ S are positive parameters and, letting [x]| F denote the jump of x over the face F , with F S,h the set of internal faces (i.e. faces with two neighbors) in N n S,h and F B,h the set of faces that are internal in N n B,h and also belong to an element in N n S,h , see Fig. 2. Finally, note that we use a Lagrange multiplier to impose the condition (2.8).
Remark 4.1. The stabilization terms j B and j S that appear in the method are consistent and are needed to control the condition number of the system matrix so that the resulting algebraic system is well conditioned independently of the position of the interface relative to the mesh. The stabilization term j S also ensures a stable discretization of the surface PDE in case the problem is convection dominated [25]. However, one may need to apply the stabilization term j B on all faces that are internal in N n B,h to guarantee a stable discretization of the bulk PDE in case the problem is convection dominated. One may also apply the bulk stabilization only on faces in F B,h and add other stabilization methods like for example SUPG [30] when the problem is convection dominated.
Remark 4.2. The method is of optimal convergence order also without the terms containing λ and µ. However, these terms are added to guarantee good conservation of the total mass of surfactants.

Implementation
Since the bulk and surface surfactant forms are coupled through a nonlinear term, see (2.7), the method (4.6) leads to a nonlinear system of equations in each time step, which we solve using Newton's method. To formulate Newton's method we define the residual F and the Jacobian DF as follows With this notation the nonlinear problem resulting from (4.6) takes the form: find u h ∈ W n h and λ ∈ R such that F (u h , λ) = 0, and the corresponding Newton's method reads: 1. Choose initial guesses u h,0 and λ 0 2. while ||(w,λ)|| > tol For t ∈ I n we choose the initial guess u h,0 to be the solution at t − n−1 i.e. u h,0 (x, t) = u h (x, t − n−1 ). Recall from equation (4.4) that a function v ∈ W n h can be written as where t ∈ I n and v B,0 and v B,1 are functions in V 0,h | N n B,h (the space of restrictions to N n B,h of functions in V 0,h ) and v S,0 and v S,1 are functions in V 0,h | N n S,h (the space of restrictions to N n S,h of functions in V 0,h ). Thus, for example the first term in DF (u, λ)(w,λ) can be written as (4.17) Our approach is to first use a quadrature rule in time and for the corresponding quadrature points compute the integrals in space. Since the geometry changes in time the contributions to the Jacobian and the residual are computed in each quadrature point for the current geometry at that point in time. Using a quadrature formula in time with quadrature weights ω q and quadrature points t q , the first term in DF (u, λ)(w,λ) is approximated by 18) The other terms in F (u, λ) and DF (u, λ)(w,λ) are treated in the same way. We use Simpson's quadrature rule for the time integration which corresponds to choosing n q = 3, quadrature points t 1 = t n−1 , t 2 = t n−1 +tn 2 , and t 3 = t n , and quadrature weights ω 1 = ω 3 = kn 6 and ω 2 = 4kn 6 . Another approach, used for example in [31], is to explicitly construct three dimensional triangulations, in case of two space dimensions, or four dimensional triangulations, in case of three space dimensions, of the space-time subdomains.
Finally, we obtain a 2(N B + N S ) + 1 × 2(N B + N S ) + 1 linear system of equations We use a direct solver to solve this linear system of equations.

Numerical examples
In all the computations in this section we use a uniform underlying mesh K 0,h consisting of triangles of size h and a constant time step of the form k = Ch. The stabilization constants τ B and τ S in the stabilization term j h are 10 −2 . We measure the convergence by using consecutive refinements of the underlying mesh and study (u B,h − u B,2h ) Ω h,1 (0.5) and (u S,h − u S,2h ) Γ h (0.5) at time t = 0.5. In cases where, e.g., u S,2h is not defined at a point in Γ h (0.5), we go to the closest point in Γ 2h (0.5) to find the concentration u S,2h . Our results show that the proposed method (see Section 4) is second order accurate and the condition number of the algebraic system of equations is bounded independently of how the interface cuts the underlying mesh. Also, the total mass of surfactants is conserved. The examples we consider in this section are from [16,13,14].

Example 1
As a first example we consider a surface convection-diffusion problem with a moving interface modeling insoluble surfactants. Thus, surfactants only exist on the interface. We consider Example 6 of [16]. Initially the interface Γ is a circle centered at the origin with radius r 0 = 1 and the velocity field β = ( (y+2) 2

3
, 0). The initial interfacial surfactant concentration u S (0, x, y) = y/r 0 + 2. The non-dimensional numbers are Pe S = 1. The computational domain is chosen as Ω = [−2, 6.4] × [−2, 2] and we use 148 gridpoints in the x-direction and 71 gridpoints in the y-direction which yields a meshsize h ≈ 0.06. The time step k = h/8 as in [16]. The surfactant concentration on the moving interface at times t = 0, 1, 2 is shown in Fig 3. In Fig. 4 the relative error in the total surfactant mass versus time is shown and we see that the error is of the order of machine epsilon. In Fig. 5 we show the relative change of the area enclosed by the interface. We observe a change in the area by less than 0.005% at time t = 2 for h ≈ 0.06.

Example 2
We now consider a coupled bulk-surface problem. The same example as in Section 4.1 of [13]. Initially the interface Γ is a circle centered in (0, 1) with radius r 0 = 0.5 and the velocity field The bulk and interfacial surfactant concentrations on the moving interface at time t = 0.5 are shown in Fig. 6. The meshsize h = 2/50 and the time step k = 0.625h. We see in Fig. 7 that the relative error in the total surfactant mass is of the order of machine epsilon and hence the total surfactant mass is accurately conserved. In Fig. 8 and Fig. 9 we show the convergence of (u B,h − u B,2h ) Ω h,1 (0.5) and (u S,h − u S,2h ) Γ h (0.5) , respectively. We observe second order convergence both in the L 2 norm and the L 1 norm.

Example 3
We consider here the same coupled bulk-surface problem as in Section 5.3 of [14]. The initial interface is a circle with radius r 0 = 0.3 centered in (0.1, 0) and the velocity field is given by (1 + cos(πx)) sin(πy), 1 2 (1 + cos(πy)) sin(πx) The bulk and surface surfactant concentrations on the moving interface at times t = 0.5, 1, 1.5, 2 are shown in Fig. 10 and Fig. 11, respectively. We see in Fig. 12 that the total surfactant mass is accurately conserved. In Fig. 13 we show the convergence of (u B,h − u B,2h ) Ω h,1 (0.5) (represented by circles) and (u S,h − u S,2h ) Γ h (0.5) (represented by stars) in the L 2 norm. We observe second order convergence which is optimal for the elements we use. We see in Fig. 14 where the condition number is shown as a function of time that as the interface evolves the condition number is bounded, independently of how the interface cuts through the mesh.

Conclusion
We have presented a new finite element method for the coupled bulk-surface problem modeling the evolution of surfactants. The model is given by a time-dependent convection-diffusion equation on the surface (or interface) coupled with a time-dependent convection-diffusion equation in the bulk. The concentration of surfactants in the bulk and the concentration of surfactants on the surface are coupled through a nonlinear term. For this problem we have proposed a space-time CutFEM. The variational formulation is written over the space-time domain. The space-time domain is divided into space-time slabs. Interpolation functions that are continuous in space but discontinuous from one space-time slab to another has been used and the discrete equations have been solved one space-time slab at a time. We believe the space-time formulation is very convenient for multiphase flow problems since good stability is expected and also unstructured meshes in time slabs can easily be used which is useful when adaptive schemes are used.
We have added consistent stabilization terms in the weak formulation which guarantee that the system matrices have bounded condition number independently of the interface position relative to the background grid and that the method is stable in case the problem is convection dominated. The numerical results demonstrate that the proposed method has optimal convergence order, leads to linear systems with bounded condition number, and conserves the total mass of surfactants independently of how the interface cuts through the fixed background mesh.
In this paper, the velocity field has been given analytically. However, by using a computed velocity field the method can be coupled to a Stokes or Navier-Stokes flow solver. In future work we will couple the method to the CutFEM proposed for the Stokes equations in [21].