A Cut Finite Element Method for two-phase flows with insoluble surfactants

We propose a new unfitted finite element method for simulation of two-phase flows in presence of insoluble surfactant. The key features of the method are 1) discrete conservation of surfactant mass; 2) the possibility of having meshes that do not conform to the evolving interface separating the immiscible fluids; 3) accurate approximation of quantities with weak or strong discontinuities across evolving geometries such as the velocity field and the pressure. The new discretization of the incompressible Navier--Stokes equations coupled to the convection-diffusion equation modeling the surfactant transport on evolving surfaces is based on a space-time cut finite element formulation with quadrature in time and a stabilization term in the weak formulation that provides function extension. The proposed strategy utilize the same computational mesh for the discretization of the surface Partial Differential Equation (PDE) and the bulk PDEs and can be combined with different techniques for representing and evolving the interface, here the level set method is used. Numerical simulations in both two and three space dimensions are presented including simulations showing the role of surfactant in the interaction between two drops.


Introduction
Surfactants are present in many forms and play an important role in many applications, for example in drug delivery for their ability of solubilizing drugs and encapsulating microbubbles with a shell or in petroleum industries for promoting oil flow [1]. Computer simulation is an important tool in such applications and is used to study the effects of surfactants on drop deformation, breakup or coalescence in multiphase flow systems.
There are several challenges that computational methods need to handle both with respect to the representation and evolution of the interface separating the immiscible fluids and the approximation of quantities such as the surfactant concentration, the fluid velocity, and the pressure. Here we assume that we know the PDEs that model two-phase flow systems in the presence of insoluble surfactant and that we also have a numerical representation technique for representing and evolving the interface. Given that, we propose a strategy based on finite element methods to accurately discretize the given PDEs both in two and three space dimensions so that the total surfactant mass is conserved and with accurate approximation of discontinuities across evolving interfaces without conforming the mesh to the time dependent domains.
Many computational methods have been developed for simulations of two-phase flow systems with surfactants, see for example [2][3][4], where discretizations based on finite differences are pro-posed, [5][6][7][8] for strategies based on boundary integral methods, [9] for a conservative discretization based on finite volume methods, and [10,11] for different strategies using finite element methods.
Our method falls into the class of so called unfitted finite element methods. We use a regular mesh that covers the computational domain but does not need to conform to outer or internal boundaries (in this case the interface). We refer to this mesh as the background mesh. On the background mesh we start with appropriate standard finite element spaces. We then define 1) active meshes corresponding to the subdomains occupied by the immiscible fluids and the interface; 2) active finite element spaces; 3) a weak formulation so that the proposed method is accurate and robust independent of the interface position relative the background mesh. Active spaces are defined as the restriction of the standard finite element spaces defined on the background mesh to the active mesh. From the active spaces we build pressure and velocity spaces with functions that are double valued on elements in a region around the interface. Hence functions in our function spaces can be discontinuous across interfaces separating the two fluids. Therefore, we don't need to regularize discontinuities by introducing regularized delta functions to approximate the surface tension force as for example in [2,4,12] or/and regularized Heaviside functions to approximate densities and viscosities across the interface as for example in [13]. In the weak formulation of the Navier-Stokes equations coupled to the surfactant transport equation, we handle space and time similarly, physical interface conditions are imposed weakly, and we add stabilization terms in the weak form that are vital for the proposed method and its straightforward implementation.
This work is a further development of the space-time Cut Finite Element Method (CutFEM) proposed in [14] to two-phase Navier-Stokes flows where also surfactant is present. We note that a Cut-FEM for the two-phase Navier-Stokes equations without surfactant has also been proposed in [15]. However, the discretization in [15] is based on implicit Euler, linear elements for both velocity and pressure, and a first order approximation of the mean curvature proposed in [16] and is therefore limited to first order accuracy. For the discretization of the incompressible Navier-Stokes equations we here start with a regular mesh and the inf-sup stable Taylor-Hood element P2-P1 pair for velocity and pressure. See [17] and [18] for stability results. We also use an integration by part result in the presence of surfactant and avoid explicit computation of the mean curvature of the interface. The proposed space-time strategy for the surfactant transport equation is similar to [19] but here we rewrite the weak formulation so that the surfactant mass is discretely conserved without enforcing the mass with a Lagrange multiplier. As in [19] and [14] the space-time integrals in the weak formulation are approximated using quadrature rules, utilizing that the stabilization term defines the extension we need to have a robust and stable method. The proposed method decouples the representation and evolution of the interface from the discretization of the PDEs and although we in this work use the standard level set method [20,21] we emphasize that other interface representation techniques can also be used. See for example [22] where a space-time cut finite element discretization is used together with an explicit representation of the interface using cubic splines.
The paper is organized as follows. We start with some basic notation in Section 2. In Section 3 we state the mathematical model and a weak formulation. In Section 4 we introduce the mesh and the finite element spaces we use in our finite element method. In Section 5 we propose a space-time cut finite element method for a surface PDE modeling the surfactant transport on an evolving interface but here we assume the velocity field is known and explicitly given. We show numerical examples and convergence studies. We then propose a space-time CutFEM for the coupled problem where the velocity field is given by the Navier-Stokes equations in Section 6 and we show simulations in both two and three space dimensions in Section 7. Finally, we conclude in Section 8.

Notation
Consider a bounded convex domain Ω in R d , d = 2, 3, with polygonal boundary ∂Ω. We assume that at any time instance t in the time interval I = [0, T ] (T ∈ (0, ∞)) the domain Ω contains two subdomains Ω i (t ) ⊂ Ω, i = 1, 2 that are separated by a sufficiently smooth closed d − 1 dimensional internal interface Γ(t ) = ∂Ω 1 (t ) ∩ ∂Ω 2 (t ). The time dependent interface is transported with a sufficiently smooth velocity field u : I × Ω → R d . We assume that for all t ∈ I , the subdomain Ω 2 (t ) is enclosed by Γ(t ) andΩ =Ω 1 (t ) ∪Ω 2 (t ) so ∂Ω ⊂ ∂Ω 1 . Let n be the unit normal to ∂Ω 1 which is outward directed with respect to Ω 1 . Thus, at the interface Γ, the normal n is pointing from Ω 1 into Ω 2 . See Figure 1 for an illustration in two space dimensions (d = 2).
For a scalar function v : I × Ω → R we use the notation v i = v| Ω i (t ) , i = 1, 2 and define the jump We denote by v t the jump of a function v in time at time instance t , i.e., We do similarly for a vector field (i.e. v : I × Ω → R d ) or a matrix-valued function but we denote those functions in bold. The usual R d gradient is denoted by ∇ and ∇ Γ denotes the surface gradient on Γ and is defined by ∇ Γ = P Γ ∇, where P Γ = I − n ⊗ n, I ∈ R d ×d is the identity matrix, and n is the unit normal to the surface Γ. The surface divergence of a vector field v(x), is defined by Note that for tangent vector fields v on Γ we have the identity ∇ Γ ·v = ∇·v, also note that v⊗∇ = (∇⊗v) T . The Laplace-Beltrami operator is denoted by ∆ Γ = ∇ Γ · ∇ Γ . Recall the Sobolev spaces and where U ⊂ R d and is in this paper Ω, Ω i (t ), or Γ(t ) for t ∈ I . We will use the notation ). We use the notation Given the time interval I we use the notation G I ,Γ to denote the following space-time surface Figure 1: Illustration of the domain Ω ∈ R 2 and the two subdomains Ω 1 (t ) and Ω 2 (t ) separated by the interface Γ(t ) at a fixed time instance t ∈ [0, T ].

The Mathematical model
In this paper we consider two-phase flow systems with insoluble surfactants present. The subdomains Ω 1 and Ω 2 are occupied by incompressible and immiscible fluids with densities ρ i ∈ R >0 and viscosities µ i ∈ R >0 , i = 1, 2. Let u : I × Ω → R d denote the fluid velocity, p : I × Ω → R denote the pressure, and w : G I ,Γ → [0, w ∞ ) denote the surfactant concentration where w ∞ is the highest surfactant concentration on Γ. We assume the dynamics is governed by the incompressible Navier-Stokes equations coupled to a time-dependent convection-diffusion equation on the interface Γ(t ) that separates the two immiscible fluids: Here, ∂ t = ∂ ∂t , (u) = (∇u + (∇u) T )/2 is the strain rate tensor, σ is the surface tension which depends on the surfactant concentration, H is the mean curvature vector defined by H = ∆ Γ x Γ with x Γ : Γ x → x ∈ R d being the coordinate map or embedding of Γ into R d , D Γ > 0 is the diffusion coefficient and is assumed to be constant, f : I × Ω → R d is a given external force (e.g. the gravitational force), and g is a given function such that ∂Ω g · n ds = 0. We can also have other type of suitable boundary conditions on ∂Ω. The initial configuration, that is the interface Γ(0) and hence the subdomains Ω i (0), i = 1, 2, the initial data u 0 : Ω → R d , and w 0 : Γ(0) → [0, w ∞ ) are also given. In this work we do not consider high Reynolds number flows.
The evolution of the interface Γ(t ) is given by the following advection equation Here φ(t , x) : I × Ω → R is the signed distance function with positive sign in Ω 2 (t ), the subdomain enclosed by the interface, and the zero level set of this function represents the interface Γ(t ), i.e., The initial condition φ(0, x) = φ 0 (x), x ∈ Ω is given by the configuration of Γ(0). Finally, we use the following linear equation of state: where σ 0 ∈ R >0 and β ∈ R >0 are given parameters. Note that β = 0 corresponds to the surfactant-free (clean) interface with σ = σ 0 being constant. Other models then this linear model describing the relation between the surface tension σ : R >0 → R >0 and the surfactant concentration w can also be used, see [23]. The numerical method we will present is not restricted to this linear relation and we discuss this in Section 6.1.
To derive the weak form we use the following variant of Reynolds' transport theorem: For w, r ∈ W we have See Lemma 5.2 of [24]. Assume that (u, p, w) is a sufficiently smooth solution to (7)- (14) and so that Reynolds' transport theorem (23) holds. Then, it follows that Following the derivation in [14] we also have that for t ∈ I Recalling the definition of the mean curvature vector H, i.e., assuming Γ(t ) is sufficiently smooth, and using integration by parts we have see also Chapter 7.6.1 in [25]. Using (27) in (25) we arrive at our weak formulation.

The mesh and finite element spaces
Here, we define the mesh and the finite element spaces we will need later when we define the finite element method. Let {T h } be a regular family of simplicial meshes of Ω with h ∈ (0, h 0 ]. We denote by h K the diameter of element K ∈ T h and we let h be the piecewise constant function that is equal to h K on element K . The mesh T h , which we will refer to as the background mesh, conforms to the fixed polygonal boundary ∂Ω but does not need to conform to the interface and hence at any time instance t ∈ I the interface Γ(t ) may cut through the mesh arbitrarily, see Figure 2. However, the background meshes we use are not completely independent of the interface as we typically use graded meshes with more elements in interface regions. For each time instance t ∈ I we denote by T h,i (t ) the collection of elements that exhibit a nonempty intersection with Γ(t ) for i = 0 and with Ω i (t ), for i = 1, 2, i.e., Let 0 = t 0 < t 1 < · · · < t N = T be a partition of I into time intervals I n = (t n , t n+1 ] of length ∆t n = t n+1 − t n , n = 0, 1, . . . , N − 1. For each time interval I n we define active meshes T n h,i , i = 0, 1, 2 and corresponding domains N n h,i . The active mesh T n h,i contains those elements in T h that cover the domain N n h,i where Denote by F n h the set of faces in the active mesh T n h,0 . Faces in F n h that are shared by two elements in the active mesh T n h,i constitute a set denoted by F n h,i . For an illustration, in two space dimensions, of the sets introduced in this section see Figure 2. Let P k (D) be the space of polynomials in D ⊂ R d , d = 1, 2, 3, of degree less or equal to k ≥ 0. Denote by V m h the space of continuous piecewise polynomials of degree less than or equal to m ≥ 1 defined on the background mesh T h , i.e., On the space-time slab I n × N n h,0 we define the space In this paper, we will propose a finite element method for the discretization of the surfactant transport equation using linear polynomials in both space and time, m = k = 1. We use the notation W n h = W n,1 h,1 . Note that a function w h ∈ W n h can be represented in the following form and for each j where ξ i , j ∈ R are coefficients and φ i (x) is the standard nodal basis function associated with mesh vertex i in the finite element space V 1 h . For the discretization of the incompressible Navier-Stokes equations we start from the LBB stable Taylor-Hood element P2-P1 pair, (V h ,Q h ), on the background mesh T h , i.e., On each of the space-time slabs I n × N n h,i , i = 1, 2, we define for k ≥ 0 the spaces is Ω 1 (t ) and the domain to the right is Ω 2 (t ).
Finally, let For k = 1 we denote these spaces by V n h = V n,1 h and Q n h = Q n,1 h . Note that functions in V n h and Q n h are double valued on elements in T n h,0 since those elements exist in both active meshes T n h,1 and T n h,2 . For an illustration, in one space dimension, of a function p h ∈ Q n h at a time instance t ∈ I n see Figure 3. Before we formulate the finite element method for the two-phase flow model we propose and study a space-time cut finite element method for the surface PDE modeling surfactant transport.

The surfactant transport equation
In this section we assume that the velocity field u is known and explicitly given and we consider the following convection-diffusion problem: Following the derivation in Section 3.1 we have for w being a sufficiently smooth solution to the differential equation (38)-(39), r ∈ W , and t ∈ I that Assuming w − (t n , x) is given, integrating in the time interval I n , and using that w In particular, note that with r = 1 in (41) we get

The Finite element method
We are now ready to present a discretization of the convection-diffusion equation on the interface. Assume the active mesh T n h,0 and the solution w − h (t n , x) from the previous space-time slab are known: where and The last term S is a stabilization term which is added to ensure that the method is stable and robust independent of how the interface is positioned relative the mesh. There are several choices for s w , see e.g. [26][27][28][29] and references therein. In this work we use where c F, j , and c Γ, j are positive constants. This stabilization was analyzed in connection with a Cut-FEM for surface PDEs on stationary surfaces in [27]. For the space-time CutFEM it is important to define the set F n h,0 , on which ghost penalty stabilization is applied, appropriately (see Figure 2). In particular note that in a space-time slab the set F n h,0 does not depend on time. The proposed stabilization term yields an extension and therefore allows us to propose a cut finite element scheme based on directly approximating the space-time integrals in the weak formulation using quadrature rules, see the next section. The proposed stabilization can be used also when the polynomial degree m > 1. For m = 1 an alternative is to extend the stabilization proposed in [26] to the space-time CutFEM, thus only having ghost penalty terms evaluated on the faces in F n h,0 , i.e., If the problem is convection dominated other stabilization terms such as streamline diffusion stabilization can be added, see e.g. [28,29].

Space-time CutFEM with quadrature in time and implementation
We propose a weak form where integrals in time are replaced with a n q -point quadrature rule with weights α n q and points t n q , q = 1, · · · , n q . Thus, we propose the following weak formulation: Given where Recall the representation of functions in W n h , see (32), and note that for example for the second term in A n h we have Furthermore, since the set F n h,0 , see definition of s w in (46), is on each space-time slab independent of time and the trial and the test functions, (w h and r h ) which are time dependent are polynomials, the quadrature rule can be chosen so that the integration of the face stabilization term is exact, i.e., In the numerical examples we use Simpson's quadrature rule. Hence the quadrature points are t n 1 = t n , t n 3 = t n+1 , t n 2 = t n +t n+1 2 , and the quadrature weights are α n 1 = α n 3 = ∆t n 6 , α n 2 = 4∆t n 6 . Note that taking r h = 1 in the proposed weak formulation (47) yields the following conservation property To implement the proposed method one needs to know the position of the interface at the quadrature points, Γ(t n q ), and also the set F n h,0 . Often we do not have the exact interface position. In this paper we use a piecewise planar approximation of the interface. Let T n = {t n , {t n q }, t n+1 }. For each t ∈ T n , we either directly make an approximation of the level set function in the space V 1 h if φ is known explicitly at all time or we discretize the advection equation (15) using the Crank-Nicolson scheme in time and continuous FEM in space with linear Lagrange elements and streamline diffusion stabilization, following [14]. By looking for sign change in each element K ∈ T h , i.e., checking the sign of this approximated level set function, φ h (t , x), x ∈ Ω, we find those elements in the mesh that are cut by the interface and thus belong to the set T h,0 (t ) for each t ∈ T n (see (28)). We compute a piecewise planar approximation Γ h (t ) at each time instance t ∈ T n by computing the intersection of the zero level set of φ h (t , x) with each cut element K ∈ T h,0 (t ).
We also need the set of faces F n h,0 where the stabilization is active. This set consist of all interior faces (faces with two neighbors) in the active mesh T n h,0 . Hence we need to know if an element K ∈ T h is in the active mesh T n h,0 or not. We say that an element K ∈ T h is in the active mesh if it is cut by the interface at some time instance t ∈ T n or if there are two time instances t k ∈ T n and t l ∈ T n where K ∩Ω 1 (t k ) = but K ∩Ω 1 (t l ) = (i.e. K ∩Ω 2 (t l ) = ). Using our approximation of the signed distance function, φ h (t , x), we say that an element K ∈ T h belongs to T n h,0 if it is in T h,0 (t ) at some time instance t ∈ T n or the sign of φ h (t k , K ) is different than φ h (t l , K ) for two time instances t k , t l ∈ T n .

Remark 5.1. Note that if Reynolds' transport theorem (23) holds, formulation (40) is equivalent to
and from (41) we have

Hence, one may instead of the bilinear form A n h in (48) use
This form was used in [30] [31][32][33].

Those formulations are based on reformulating space-time integrals in the weak formulation to surface integrals over the space-time manifold in R d +1 (when Γ is a surface in R d ). In our formulation we never approximate surface integrals in R d +1 . The implementation of the proposed space-time cut finite element method is straightforward starting from an implementation of CutFEM for stationary domains but relies on the stabilization terms and also on a method for
representing and evolving the interface.

Numerical examples
We consider one example in two space dimensions and two examples in three space dimensions and validate that the proposed unfitted finite element method is accurate and conserves the surfactant mass at all time instances t 0 , t 1 , · · · , t N . In the examples in this section, the velocity field u(t , x) is explicitly given, we always use a uniform time step size ∆t = T /N , a uniform background mesh of the computational domain, and solve the resulting linear systems with a direct solver. Linear elements in both space and time are used. The proposed method uses the weak formulation (47) with A n h (u, w h , r h ) as in (48) and is refered to as the conservative formulation. When A n h (u, w h , r h ) based on (55) is used we refer to the method as the non-conservative method. The stabilization parameters C F,1 and C Γ,1 are in all examples equal to 10 −2 . We study the accuracy of the proposed cut finite element method by studying the following L 2 -error at the final time We also report the following conservation error where we use the same quadrature rule as when we approximate the space-time integrals in the weak formulation, see Section 5.1.1. Hence the quadrature points t n q and weights α n q are from Simpson's quadrature rule.

Example 1 -2D
We consider the convection diffusion equation (38)-(39) with parameters and solution chosen as in [34]. Initially the interface, Γ(0), is a circle centered at origin with radius one. The interface is evolving with the velocity field : The corresponding level-set function is The diffusion coefficient D Γ = 1, the right hand side f and initial solution w 0 in equation (38) Figure 5. Results for both the proposed scheme and the non-conservative formulation are shown. For both formulations, using m = k = 1, we see the expected second order convergence in the L 2 (Γ h (t N ))-norm, similar to the results reported for the unfitted finite element method in [34] (see Table 6.3 of [34]). We also study the conservation of mass. The conservation error e c , defined in (57), of the proposed scheme is shown in Figure 6a and results using the non conservative method are shown in Figure 6b. For the proposed formulation the magnitude of the conservation error is of the order of machine epsilon. In the non-conservative formulation the error e c depends on the drop deformation, the mesh size, the time step size, and the quadrature rule used.

Example 2 -accuracy test in 3D
We follow Example 1 in [33] and study the accuracy of the proposed space and time discretization.  a sphere centered at origin with radius equal to 1.5 and the concentration on Γ(0) is w 0 (x) = 1 + x y z. The velocity field is The corresponding level-set function is The diffusion coefficient D Γ = 1 and the right hand side of equation ( In Figure 7a we show the error in the L 2 (Γ h (t N ))-norm (t N = 1), see (56), versus mesh size h for different fixed time step sizes ∆t while in Figure 7b we show the error versus the time step size ∆t for different fixed mesh sizes h. We see that our scheme has a second order convergence rate. We also see that in this example the error is often dominated by the error in the space discretization. In our scheme in particular the geometric error, the error coming from the approximation of the interface, is dominating in this example. The results in Figure 7a-7b can be compared with Figure 7.1 in [33] where a different second order space-time unfitted finite element method is used, see Remark 5.2. However, note that for the results we report here we have approximated the level set function and the interface on the same background mesh as our spaces are defined on while the results reported in [33] for h and ∆t use an approximation of the interface that is on one regular refinement of the space-time mesh (hence with h/2 and ∆t /2). We observe that therefore when the error in the space discretization dominates we have slightly larger errors than reported in [33] but when the error in the time discretization dominates (seen when larger time step sizes are used) the errors from the proposed scheme are smaller. We study the conservation error in the next example.

Example 3 -conservation test in 3D
In this last example, also from [33], we study the conservation property of the proposed scheme in an example in three space dimensions. The initial interface Γ(0) is given by the zero level set of the following function φ(0, x) = (x − z 2 ) 2 + y 2 + z 2 − 1 and the initial concentration on Γ(0) is w 0 (x) = 1 + x y z. The velocity field with which the interface is also transported is β(t , x) = (0.1x cos(t ), 0.2y sin(t ), 0.2z cos(t )) and the right hand side f is zero.  Figure 8 for h = 0.125. We show the discrete total mass M h (t ) = Γ h (t ) w h (t , x) ds as function of time and the conservation error e c defined in (57) (with f = 0) in Figure 9. Compared to the space-time unfitted finite element method in [33] (see Figure 7.5. of [33]) the proposed method conserves the discrete mass.

The Finite element method for the two-phase flow system
We are now ready to propose a discretization of the PDEs in Section 3 describing a two-phase flow system with insoluble surfactants. Starting from the initial condition u − h (t 0 , x) = u 0 (x) in Ω 1 (t 0 )∪Ω 2 (t 0 ) and w − h (t 0 , x) = w 0 (x) on Γ(t 0 ) we solve the following variational formulation, where space and time are treated similarly and interface and boundary conditions are imposed weakly, one space-time slab at a time: for I n , n = 0, 1, · · · , N − 1 given u − h (t n , x) and w − h (t n , x) find u h ∈ V n h , p h ∈ Q n h , and w h ∈ W n h such that for all v h ∈ V n h , q h ∈ Q n h , and r h ∈ W n h . Here the form A(u h , w h , r h ) is defined as in (44), where the forms a (t , u h , v h ) and b(t , v h , p h ) are defined as in equation (19)- (20) with the penalty parameters λ Γ , λ ∂Ω , and the weights k 1 , k 2 in the averaging operators (22) chosen as in [14]. Further, with l (t , v h , q h , r h ) defined as in (21) and we also have where the terms s p (p h , q h ), s u (u h , v h ), and s w (t , w h , r h ) are appropriate stabilization terms. We choose s w (t , w h , r h ) as in (46) and Here C p and C u,m are positive constants and v F denotes the jump of a function v at the face F and is defined as x ∈ F , and n F is a fixed unit normal to F . The stabilization terms s p and s u are similar as in [35] but the sets F n h,i are different in the space-time method. Note that for each space-time slab the sets F n h,i , i = 0, 1, 2 do not depend on t , see Figure 2 but the trial and test functions do depend on time. These stabilization terms define an extension and therefore including them in the weak form allows us to formulate our space-time method based on quadrature in time and have a stable and robust discretization independent of how the interface cuts through the background mesh.
Note that from the weak formulation we find functions p h = (p h,1 , p h,2 ) ∈ Q n h and u h = (u h,1 , u h,2 ) ∈ V n h which are double valued in the interface region, see Section 4. We use (u h , p h ) and define an approximate velocity fieldũ h and pressurep h by:

Implementation
The proposed discretization (58) leads to a non linear system which we solve using Newton's method: -Choose initial starting guess (u h,0 , p h,0 , w h, Here Since we use the linear equation of state (17), we have σ (w h ) = −σ 0 β and If for example the Langmuir model is used which is a nonlinear model that describe the relation between the surface tension and the surfactant concentration, i.e., where σ 0 ∈ R >0 , β ∈ R >0 , and w ∞ ∈ R >0 are given parameters, we instead have σ (w h ) = −β w ∞ −w h . In the first step of the above algorithm we choose the initial guess (u h,0 (t ), p h,0 (t ), w h,0 (t )) to be equal to the solution from the previous space-time slab at time t n , i.e., ). In the assembly of the linear system in step i. we approximate each space-time integral in F and DF by first using a quadrature rule in time and then at each quadrature point in time the integral in space is computed. We approximate all the integrals in time using Simpson's quadrature rule. See also Section 5.1.1.
Here, the exact level set function is not known explicitly except at time t = 0. We find an approximation φ h (t , x) at each time instance t ∈ T n = {t n , {t n q }, t n+1 } by discretizing the level set equation (15) using the Crank-Nicolson scheme in time and continuous FEM with quadratic Lagrange elements in space and streamline diffusion stabilization following [14]. As in Section 5.1.1 we find a piecewise planar Γ h (t ) as the zero level set of φ h (t , x) ∈ V 1 h , for all t ∈ T n . We then also have approximations Ω h, 1 and Ω h,2 of the two subdomains Ω 1 and Ω 2 . The contribution from integration on Ω h,i (t n q ) ∩ K is divided into contributions on one or several triangles in two dimensions and tetrahedra in three dimensions depending on how the interface cuts element K .

Numerical examples
We consider several examples in two space dimensions including drop-drop interactions and one example in three space dimensions. In the last simulation, in three space dimensions, we use the proposed scheme but instead of linear element in time we use piecewise constant polynomials in time, i.e. k = 0 (see Section 4). This corresponds to using the implicit Euler method in time (see Remark 3.2 in [14]). The stabilization parameters are all chosen to be 10 −2 . The resulting linear systems are always solved with a direct solver.

Rising drop
Consider a two-dimensional drop rising in a liquid column due to buoyancy with physical parameters as in Table 1. The force of gravity is f = ρ(0, −0.98). This is a benchmark test case from [36] but we also add surfactant as in [11]. The computational domain in space is Ω = [0, 1] × [0, 2] and the drop is Figure 10: Initial configuration initially a circle centered at (0.5, 0.5) with radius equal to 0.25. The no-slip boundary condition, u = 0, is imposed on the horizontal walls and the free slip condition, u · n = 0 , τ · 2ε(u)n = 0, is imposed on the vertical walls (here n and τ are the unit normal-and tangent vectors to ∂Ω). For an illustration of the initial setup see Figure 10. 1000  100  10  1 24.5 0.5 1 0.1 We use a background mesh with more elements in the region where the drop is evolving (the mesh can be seen in Figure 14). The coursest mesh size is h outer = 1/40 while the finest mesh size is h inner = 1/80. We choose a time step size ∆t = h outer /4. The evolution of the drop and the surfactant concentration is shown in Figure 11. One can see that surfactant accumulate on the bottom and increase the deformation of the interface compared to the case of a clean interface, β = 0.
We compute the following benchmark quantities: • Center of mass Here, the second component y c is of interest.
• Circularity where P b is the perimeter of the drop and P a is the perimeter of the circle which has an area equal to that of the drop. • Rise velocity Here, we are interested in the velocity component u 2 c which is in the direction opposite to the gravitational vector f.
In Figure 12 these quantities are shown versus time t both when surfactant is present and for a clean interface which we studied carefully in [14]. When surfactant is present, the deformation increases and the rise velocity decreases. This impact also the center of mass. In Table 2 we report some characteristic values. Our results agree well with the results reported for the parametric finite element method in [11] (See Table 2 and Figure 7 and 8 in [11]).  In Figure 13 we show the conservation error which is of the same order as machine epsilon. Thus, the proposed scheme conserves the total surfactant mass. Finally, in Figure 14   mated discontinuous pressure at t = 2. We see that the discontinuity across the evolving interface is accurately approximated.

Drop in shear flow
We now consider an initially circular drop, centered at the origin with radius one, in a shear flow experiment from [12]. We study the influence of surfactant on the deformation of the drop. The computational domain is Ω = [−5, 5] × [−2, 2] and the Dirichlet boundary condition u(t , x) = (0.5y, 0) is imposed on ∂Ω. The physical parameters are chosen as in Table 3 and f = 0.
The simulations in this section have all been done on the mesh shown in Figure 15. The coursest elements have size h outer = 0.13 and the smallest element size is h inner = 0.05. The time step has been chosen equal to 0.015. In Figure 16 the shape of the drop and the surfactant concentration at different time instances, t=0, 4, 8, 12, is shown. One can see that when β increases, the drop becomes more elongated and narrower. Our results with the proposed CutFEM agree well with Figure 1 of [12] where  the simulations are made with an immersed boundary method with h = 0.02 and ∆t = h/8 = 0.0025 and with Figure 11 of [11]. In Figure 17 we show different quantities of interest for test case 3 when β = 0.5. First we see in Figure 17a that the total surfactant mass is conserved during the simulation. We show the error in the area of the drop in Figure 17b. The level set method and the discretization we use does not conserve the area of the drop and for the mesh and time step size that has been used in the simulation the error in the conservation of the area is of order 10 −4 . Finally, Figure 17c shows the perimeter of the drop. These results can be compared with Figure 5 in [12].
The background mesh is shown in Figure 18, where close to the boundaries, ∂Ω, the mesh size is 1/25 and around the origin the finest elements are of size 1/100. In Figure 19 and 20 we show the drop shape in the two test cases. In the first case, with clean interfaces, the two drops coalesce during   contact and form a single drop but in the second case, in the presence of surfactant (β = 0.6), the two drops become more elongated and pass each other without merging. Note that here topological changes such as merging is easily handeled by the level set method.

Rising drop in 3D
In this section, we consider the first 3D example in [11]. The interface is initialy a sphere with radius 0. 25 Table 1.
We compute the same quantities as in 7.1 but with their natural extensions to three space dimensions. We define the circularity as where S a is the surface area of the sphere which has the same volume as the drop and S b is the surface area of the drop. The computations are done on a uniform mesh of size h = 1/30 with a time step ∆t = 0.02. One can see the domain and the final shape of the drops in Figure 21. As in two space dimensions the presence of surfactant increase the deformation of the drop and decrease the rise velocity. Table 5 shows some characteristic quantities obtained for the simulation with and without surfactant. The results agree with the results obtained in [11].

Conclusions
In this work our focus was on developing an accurate discretization for the surfactant transport equation coupled to the incompressible Navier-Stokes equations. The proposed discretization of the surface PDE and the bulk PDEs utilize the same computational mesh which does not need to conform to the evolving interfaces. We have presented a new space-time weak formulation for the surfactant transport equation which results in discrete conservation of surfactant mass without enforcing the mass with a Lagrange multiplier. In this paper we considered insoluble surfactants but in future work we will also consider soluble surfactants. In the presented scheme the computation of the mean curvature vector is avoided and this simplifies the computations. The presented space-time scheme is based on using quadrature rule in time to approximate the space-time integrals and has a straightforward implementation. We used linear elements in both space and time and a piecewise linear approximation of the interface from an approximate level set function and observed second order convergence. However, if better approximations of the interface is used the proposed space-time CutFEM method can be used with high order elements and result in higher order discretization than second order.
Our scheme relies on an accurate representation and evolution of the interface. We used the level set method which simplified the representation and evolution of the interface in three space dimensions and in the simulations when coalescence occured. However, we note that with our approximate level set function we do not conserve the total area of the drop and in future work we would like to improve on the representation of the interface.