A conforming-nonconforming mixed immersed finite element method for unsteady Stokes equations with moving interfaces

1.   Introduction

Let $\Omega\subset \mathbb{R}^2$ be an open bounded domain separated by a time-dependent smooth interface $\Gamma(t)$. The evolving interface $\Gamma(t)$ divides the domain $\Omega$ into two open subdomains $\Omega^+(t)$ and $\Omega^-(t)$ such that $\Omega = \Omega^+(t)\cup\Omega^-(t)\cup \Gamma(t)$, see Figure 1. Consider the following initial-boundary-value problems of the Stokes equation

where $\mathbf{u}$ and $p$ denote the flow velocity and the pressure, respectively. Functions $\mathbf{f}$, $\mathbf{u}_0$, and $p_0$ are the given body force, the initial velocity, and the initial pressure, respectively. $\mathbb{I}$ denotes the identity tensor. The movement of the interface is assumed to be guided by a given velocity field $\mathbf{v}(\mathbf{x}, t)$ as follows

The viscosity function $\mu(\mathbf{x})$ is assumed to have a finite jump across the interface $\Gamma(t)$. For simplicity, we assume that $\mu(\mathbf{x})$ is a piecewise constant function

where $\mu^{\pm}>0$ and $\mathbf{x} = (x, y)$. At any time $t$, the velocity and the stress tensors satisfy the following homogeneous interface jump conditions

where the jump $\left[\!\left[{\mathbf{v}(\mathbf{x})}\right]\!\right]_\Gamma : = \mathbf{v}^+(\mathbf{x})|_\Gamma-\mathbf{v}^-(\mathbf{x})|_\Gamma$, and $\mathbf{n}$ denotes the unit normal vector to the interface $\Gamma$ pointing from $\Omega^-(t)$ to $\Omega^+(t)$.

Numerical approximations of the Stokes equation have been extensively studied for many years due to its wide applicability to model natural phenomena such as airflow, water flow, and ocean currents. The family of Taylor-Hood finite elements [35] uses conforming $P_k$-$P_{k-1}$ pairs to approximate the velocity and pressure requiring the polynomial degree $k\ge 2$. To use low-order polynomials and preserve the mass conservation property, nonconforming finite element methods [4] and discontinuous Galerkin methods [33] have been widely used. Crouzeix and Raviart introduced the lowest-order nonconforming $P_1$-$P_0$ finite element [6], which is well-known as the CR finite element. On quadrilateral meshes, Rannacher and Turek developed a nonconforming rotated-$Q_1$ element in [32]. These nonconforming finite elements make use of low-order polynomials and they are elementwise divergence-free [2,23]. A mixed conforming-nonconforming finite element space was introduced in [21]. In this discretization, one velocity component is approximated by the conforming $P_1$ element and the other one is approximated by the nonconforming $P_1$ element, while the pressure is piecewise constant. This mixed FEM construction has advantages over the CR element in handling the Neumann boundary condition [21], and it is computationally less expensive than the CR element. This mixed conforming-nonconforming FEM has also been extended to the 3D Stokes equation [38]. For more details regarding numerical methods for Stokes equations, we refer readers to [8,19].

Traditional numerical methods use interface-fitted meshes for solving interface problems. For fluid flow interface problems, the arbitrary Lagrangian-Eulerian (ALE)-based finite element is a popular numerical method [7,22,37]. Recently, there has been a growing interest in developing unfitted-mesh numerical methods for solving a variety of interface problems, see Figure 2. Comparing with conventional fitted-mesh methods, such as classical FE and DG methods, the unfitted-mesh methods do not require the alignment of the mesh with a prescribed nontrivial interface; hence it is more desirable for time-dependent problems with moving interfaces. In the past decades, several unfitted-mesh methods have been developed for solving Stokes interface problems, such as CutFEM [15], Nitsche's FEM [36], XFEM[9], fictitious domain FEM [31,34], to name only a few. The immersed finite element method (IFEM) [24,26,18,11,14,30] is a class of unfitted-mesh finite element methods for solving interface problems. The main idea of IFEM is to incorporate the interface jump conditions in the construction of IFE basis functions. Unlike other aforementioned unfitted-mesh methods, the IFE space is isomorphic to the standard FE space with no interface. Consequently, not only is the mesh independent of the interface in an IFEM, but also the number and the location of the degrees of freedom are interface-independent. For time-dependent interface problems with a moving interface, the linear system has the same size at each time level and the nonzero entries remain at the same locations [10,12,13,16,17]. Moreover, the method-of-lines technique can be utilized together with IFEM for solving moving interface problems [25].

There have been some IFE methods developed for steady-state Stokes interface problems. In [1,3] the $Q_1$-$Q_0$ immersed DG method was introduced. The velocity is approximated by the broken $Q_1$ functions while the pressure is approximated by the piecewise constant functions. The computational framework is based on the interior penalty DG method [33]. Based on the nonconforming finite element framework [6,32], a class of nonconforming IFE approximations was developed [20]. Recently, a $P_2$-$P_1$ Taylor-Hood IFE space was introduced in [5]. The partially penalized IFE scheme is used with ghost penalty for enhancing the stability of numerical scheme especially for the pressure approximation.

The goal of this paper is two-fold. First, we develop a lowest-order conforming-nonconforming mixed IFE space for the Stokes equation based on [21]. Comparing with the IDG method [1] and the Taylor-Hood IFE method [5], our new IFE method has no additional consistency and stability terms, so the numerical formulation is much simpler to implement. Comparing with the CR-$P_0$ IFE space [20], there are significantly less degrees of freedom due to the conformity of one velocity component. In fact, on the same triangular mesh, only two-thirds of degrees of freedom are required for velocity in this new mixed IFEM. Besides, the mixed conforming-nonconforming finite element is robust for handling both Dirichlet and Neumann boundary conditions, while the CR finite element space is only stable for Dirichlet boundary conditions [21].

The second goal is to apply this mixed IFE method for solving unsteady Stokes equations with a moving interface. We will use the new vector-valued IFE spaces for semi-discretization, and use the prototypical backward-Euler and Crank-Nicolson scheme for full-discretization. Our method does not require re-meshing at any time level. Since the degrees of freedom are also independent of the interface, there is no need to overhaul the global matrices at each time level. Instead, only local modification is carried out on elements where the interface configuration changed during two consecutive time steps.

The rest of the paper is organized as follows. In Section 2, we construct the new mixed IFE spaces for Stokes equations. In Section 3, we report some fundamental properties of the new IFE spaces. In section 4, we present the semi-discrete and the full-discrete IFE method for solving unsteady Stokes interface problems with a moving interface. Some numerical examples are reported in Section 5. A brief conclusion is given in Section 6.

2.   Mixed conforming-nonconforming immersed finite element spaces

In this section, we introduce the mixed conforming-nonconforming IFE spaces for Stokes equations. Let $\mathcal{T}_h = \{T\}$ be an interface-unfitted triangulation of a polygonal domain $\Omega$. Let $\mathcal{N}_h$ and $\mathcal{E}_h$ denote the collections of nodes and edges of the mesh $\mathcal{T}_h$, respectively. Elements in $\mathcal{T}_h$ are divided into two categories: an interface element if $T$ is cut through by the interface $\Gamma$, and a non-interface element otherwise. The collections of interface elements and non-interface elements are denoted by $\mathcal{T}^i_h$ and $\mathcal{T}^n_h$, respectively. Similarly, for each edge $e\in\mathcal{E}_h$, if $e$ intersects the interface, it is called an interface edge; otherwise it is a non-interface edge. The collections of interface edges and non-interface edges are denoted by $\mathcal{E}_h^i$ and $\mathcal{E}_h^n$, respectively. Additionally, we let $\mathring{\mathcal{E}}_h$ and $\mathcal{E}_h^b$ be the collections of internal edges and boundary edges, respectively. Let $\mathring{\mathcal{N}}_h$ and $\mathcal{N}_h^b$ be the collections of internal nodes and boundary nodes, respectively. We also assume that the triangulation $\mathcal{T}_h$ satisfies the following hypotheses [28]:

● (H1) The interface $\Gamma$ cannot intersect an edge of any element at more than two points unless the edge is part of $\Gamma$.

● (H2) If $\Gamma$ intersects the boundary of an element at two points, these intersection points must be on different edges of this element.

● (H3) The interface $\Gamma$ is a piecewise $C^2$-continuous function, and the mesh $\mathcal{T}_h$ is formed such that the subset of $\Gamma$ in every interface element is $C^2$-continuous.

2.1.   Conforming-nonconforming FE spaces

Let $T\in\mathcal{T}_h^n$ be a non-interface element with vertices $A_1$, $A_2$, $A_3$ oriented counterclockwise. We label the edges of $T$ by $e_1 = \overline{A_1A_2}$, $e_2 = \overline{A_2A_3}$, and $e_3 = \overline{A_3A_1}$. Let $\lambda_{j, T}\in \mathcal{P}_1$ be the Lagrange linear nodal basis functions such that

where $\delta_{ij}$ is the Kronecker function. Define $\psi_{j, T} = 1- 2\lambda_{k_j, T}$ with $k_1 = 3$, $k_2 = 1$, and $k_3 = 2$. It can be verified that $\psi_{j, T}$ satisfies the mean-value conditions, namely,

Thus $\psi_{j, T}$, $j = 1, 2, 3$ are nonconforming-$\mathcal{P}_1$ (CR) basis functions on $T$. The pressure is approximated by the piecewise constant function space denoted by $\mathcal{P}_0$. On each non-interface triangle $T\in\mathcal{T}_h^n$, the vector-valued CR-$P_1$-$P_0$ finite element space can be written as $\mathbf{S}_{h}^n(T) = \mathcal{P}_1\times\mathcal{P}_1\times \mathcal{P}_0$, or quivalently, $\mathbf{S}_{h}^n(T) = \text{span}\{\mathit{\boldsymbol{\psi}}_{i, T}: 1\le i\le 7\}$ where the vector-valued basis functions are given below

Similarly, we can also form the $P_1$-CR-$P_0$ finite element space using conforming-$\mathcal{P}_1$ bases for the first component, and the nonconforming-$\mathcal{P}_1$ in the second component, then the basis functions are

The $P_1$-CR-$P_0$ finite element space is $\tilde{\mathbf{S}}_{h}^n(T) = \text{span}\{\mathit{\boldsymbol{\tilde\psi}}_{i, T}: 1\le i\le 7\}$. Note that these two spaces ${\mathbf{S}}_{h}^n(T)$ and $\tilde{\mathbf{S}}_{h}^n(T)$ are identical, both equal $\mathcal{P}_1\times\mathcal{P}_1\times \mathcal{P}_0$. However, the degrees of freedom of ${\mathbf{S}}_{h}^n(T)$ and $\tilde{\mathbf{S}}_{h}^n(T)$ are different, as indicated in (9) and (10). For more details of the conforming-nonconforming finite elements, we refer readers to [21].

2.2.   Mixed conforming-nonconforming IFE spaces

In this subsection, we extend these conforming-nonconforming finite elements to the IFE spaces on each interface triangle $T\in\mathcal{T}_h^i$. Let $A_i = (x_i, y_i)$, $i = 1, 2, 3$ be the vertices of $T$. Without loss of generality, we consider the reference triangle whose vertices are given by

Note that an arbitrary triangle with vertices $A_i = (x_i, y_i)$, $i = 1, 2, 3$ can be mapped to this reference triangle by the following mapping

To simplify the notation, we still use $x$, $y$, rather than $\hat x$, $\hat y$ on the reference triangle. According to the hypotheses (H1)-(H3), there are two distinct intersection points on each interface triangle, denoted by $D = (x_d, y_d)$ and $E = (x_e, y_e)$, on two different edges. There are generally three types of interface triangles as depicted in Figure 3. The line segment $\overline{DE}$ is used to approximate the actual interface curve $\Gamma\cap T$, and it divides the element $T$ into two subelements, denoted by $T^+$ and $T^-$. For example, on a Type Ⅰ interface element, $D = A_1+d(A_2-A_1)$ and $E = A_1+e(A_3-A_1)$ where $0<d, e<1$.

We construct the vector-valued IFE shape functions in terms of the FE functions $\mathit{\boldsymbol{\psi}}_{i, T}$ in (11). To be more precise, we have

It can be observed that each vector-valued IFE shape function $\mathit{\boldsymbol{\phi}}_{j, T}$ has $14$ unknown coefficients $c_{ij}^s$, with $1\le i\le 7$ and $s = +, -$. These coefficients are determined by seven local degrees of freedom (prescribed nodal values, edge values, and the mean pressure value), six interface jump conditions, and a divergence free condition stated below.

● Three edge-value conditions:

● Three nodal-value conditions:

● One mean-pressure-value condition:

● Four continuity conditions of the velocity to incorporate (7):

● Two stress continuity conditions to incorporate (8):

● One continuity of the divergence condition to incorporate (2):

Here, in (18)-(21), the scalar function $\phi_{i, j}$ denotes the $i$-th component of $\mathit{\boldsymbol{\phi}}_{j, T}$. More precisely, we have $\mathit{\boldsymbol{\phi}}_{j, T} = (\phi_{1, j}, \phi_{2, j}, \phi_{p, j})$ such that $\mathit{\boldsymbol{\phi}}_{j, T}|_{T^s} = \mathit{\boldsymbol{\phi}}_{j, T}^s = (\phi_{1, j}^s, \phi_{2, j}^s, \phi_{p, j}^s)$ $\in \mathcal{P}_1\times \mathcal{P}_1\times \mathcal{P}_0$, with $s = +, -$. Combining the conditions (15)-(21) yields a linear system of fourteen unknowns. On Type Ⅰ interface element, we have

where the coefficient matrix $M_\text{I}$ is written as the first seven columns and the next seven seven columns due to width limit of the page:

and

with $\rho = \mu^+/\mu^-$ being the jump ratio. The unknown vector $\mathbf{c}_j$ and the right-hand-side vector $\mathbf{e}_j$ take the form

We can obtain the vector-valued IFE shape functions $\mathit{\boldsymbol{\phi}}_{j, T}$ by solving for $\mathbf{c}_j$ with each vector $\mathbf{e}_j$, $j = 1, 2, \cdots, 7$. Note that the matrices for Type Ⅱ and Type Ⅲ interface elements, denoted by $M_\text{II}$ and $M_\text{III}$, can be derived in a similar fashion; hence, we omit the details in this paper.

As an illustration, we plot the three components of the CR-$P_1$-$P_0$ IFE shape function $\mathit{\boldsymbol{\phi}}_{4, T}$ in Figure 4. As a comparison, we plot the standard CR-$P_1$-$P_0$ FE shape function $\mathit{\boldsymbol{\psi}}_{4, T}$. We note that both FE and IFE shape functions are such that their second velocity components have the value one at the node $A_1$. However, due to the coupled stress jump condition (8), the first velocity component and the pressure component of the IFE shape function $\mathit{\boldsymbol{\phi}}_{4, T}$ are not completely zero, as the FE shape function. This is a similar phenomenon that also occurs in other vector-valued IFE functions [1,20,27,29].

The local CR-$P_1$-$P_0$ IFE space is formed by $\mathbf{S}_{h}^i(T) = \text{span}\{\mathit{\boldsymbol{\phi}}_{j, T}: 1\le j \le 7\}$, and the global CR-$P_1$-$P_0$ IFE space is defined to be

C1: $\mathit{\boldsymbol{v}}|_T\in \mathbf{S}_h^n(T), \; \forall T\in \mathcal{T}_h^n$, and $\mathit{\boldsymbol{v}}|_T\in \mathbf{S}_h^i(T)\; \forall T\in \mathcal{T}_h^i$.

C2: $\int_e \left[\!\left[{v_1}\right]\!\right] ds = 0, \; \; \forall e\in \mathring{\mathcal{E}}_h$.

C3: $v_2$ is continuous at every internal point $(x, y)\in\mathring{\mathcal{N}}_h$.

We can construct the $P_1$-CR-$P_0$ IFE space in a similar manner. In this case the edge-value conditions (15) are imposed on the second velocity component, and the nodal-value conditions (16) will apply to the first velocity component. The remaining conditions (17)-(21) are the same. Let $\tilde{\mathbf{S}}_{h}^i(T)$ be the local $P_1$-CR-$P_0$ IFE space, then the corresponding global IFE space $\tilde{\mathbf{S}}_{h}(\mathcal{T}_h)$ is defined as follows

C4: $\mathit{\boldsymbol{v}}|_T\in \mathbf{\tilde S}_h^n(T), \; \forall T\in \mathcal{T}_h^n$, and $\mathit{\boldsymbol{v}}|_T\in \mathbf{\tilde S}_h^i(T)\; \forall T\in \mathcal{T}_h^i$.

C5: $v_1$ is continuous at every internal point $(x, y)\in\mathring{\mathcal{N}}_h$.

C6: $\int_e \left[\!\left[{v_2}\right]\!\right] ds = 0, \; \; \forall e\in \mathring{\mathcal{E}}_h$.

Remark 1. In many cases, the momentum equation (1) of the Stokes system is written as

where the stress is expressed using the strain tensor $\mathbf{\epsilon}(\mathbf{u}) = (\nabla \mathbf{u} + (\nabla\mathbf{u})^t)/2$. In this setup, the stress jump condition (8) should also be changed to

Since the viscosity coefficient $\mu(\mathbf{x})$ is a (piecewise) constant, the incompressibility condition (2) yields

Hence, these two equations are equivalent in this case. In construction of IFE shape functions, only (19)-(20) need to be replaced by the following two conditions

The local CR-$P_1$-$P_0$ IFE functions, denoted by $\mathit{\boldsymbol{\phi}}^\epsilon_{j, T}$, and local $P_1$-CR-$P_0$ IFE functions, denoted by $\mathit{\boldsymbol{\tilde\phi}}^\epsilon_{j, T}$ can be constructed accordingly. The corresponding global IFE spaces are denoted by $\mathbf{S}_{h}^\epsilon(\mathcal{T}_h)$ and $\mathbf{\tilde S}_{h}^\epsilon(\mathcal{T}_h)$.

3.   Properties of the mixed conforming-nonconforming IFE spaces

In this section, we present some basic properties of the mixed conforming-nonconforming IFE spaces.

Theorem 3.1 (Unisolvency). The CR-$P_1$-$P_0$ IFE shape functions $\mathit{\boldsymbol{\phi}}_{j, T}$, $1\le j\le 7$ can be uniquely determined by the prescribed edge values, the nodal values, and the mean pressure value, regardless of the interface locations and the jumps of viscosity coefficients $\mu^\pm>0$.

Proof. We show the unisolvency by considering the invertibility of the coefficient matrices $M_\text{I}$, $M_\text{II}$, and $M_\text{III}$. For the Type Ⅰ interface triangle, by direct calculation we have

For the Type Ⅱ interface element, we have

where

and with $s = 1-d$, we have

For the Type Ⅲ interface element, we have

where

and with $s = 1-d$ and $t = 1-e$,

The determinants of coefficient matrices are uniformly nonzero for all $0\le d\le 1$, $0\le e\le 1$, $\rho>0$ and for all three types of interface elements. This ensures the unisolvency of the IFE functions.

The following theorems provide basic properties of the new IFE functions. The proofs of these results can be verified by direct calculation, hence we omit the proof in this paper. For more details, we refer the readers to some earlier references [1,20].

Theorem 3.2 (Consistency). Let $T\in\mathcal{T}_h^i$ be an interface triangle.

● If $\mu^+ = \mu^-$, the IFE shape functions $\mathit{\boldsymbol{\phi}}_{j, T}$ become the FE shape functions $\mathit{\boldsymbol{\psi}}_{j, T}$, $1\le j\le7$.

● If the interface moves out of a triangle $T$, i.e.,

the IFE shape functions $\mathit{\boldsymbol{\phi}}_{j, T}$ become the FE shape functions $\mathit{\boldsymbol{\psi}}_{j, T}$, $1\le j\le7$.

Remark 2. The consistency (29) enables us to use IFE functions for solving Stokes moving interface problem efficiently. In fact, as the interface moves out of an element, the IFE functions smoothly convert to the FE functions. No extra condition is needed to enforce this transition.

Theorem 3.3 (Continuity of Velocity). Let $T\in\mathcal{T}_h^i$ be an interface element and $\mathit{\boldsymbol{\phi}}_{j, T}$ be the vector-valued shape functions. Then the velocity components $\phi_{i, j}\in C(T)$, for $i = 1, 2$, and $j = 1, 2, \cdots, 7$.

Theorem 3.4 (Partition of Unity). Let $T\in\mathcal{T}_h^i$ be an interface element. The vector-valued IFE shape functions $\mathit{\boldsymbol{\phi}}_{j, T}$, $j = 1, 2, \cdots, 7$, satisfy the partition of unity property, namely: for any $(x, y)\in T.$

4.   Semi-discrete and full-discrete schemes

In this section, we first derive the weak form of the unsteady Stokes interface problem (1)-(8), and then develop the semi-discrete and full-discrete IFE schemes. We use $(\cdot, \cdot)_\omega$ to denote the $L^2$ inner product on a subset $\omega\subset\Omega$. We will omit the subscript $\omega$ if $\omega = \Omega$.

4.1.   Weak formulation

Taking the inner product with $\mathbf{v}\in [H_0^1(\Omega)]^2$ on the equation (1) and integrating by parts over $\Omega^-$ yields,

Here the second term is the inner product of two tensors $\mathbb{A} = [A_{ij}]$ and $\mathbb{B} = [B_{ij}]$, which is defined by $(\mathbb{A}, \mathbb{B}) : = \sum_{i, j} (A_{ij}, B_{ij})$. Note that $\mathbf{n}_\Gamma$ is pointing from $\Omega^-$ to $\Omega^+$ and $\mathbf{v}$ vanishes on the outer boundary $\partial\Omega$. We have

Similar argument applying to the subdomain $\Omega^+$ yields

Adding the above two equations together, and applying the interface jump condition (8), we have

Multiplying $q\in L^2(\Omega)$ to (2), and integrating by parts we have

Define the bilinear form and the linear form

Here, $L_0^2(\Omega) = \{q\in L^2(\Omega): \int_\Omega qd\mathbf{x} = 0\}$. The weak form of the unsteady Stokes interface problem (1)-(8) is given as follows.

Weak Form: Find $\mathbf{u}\in H^1(0, T;[H_0^1(\Omega)]^2)$ and $p\in L^2(0, T;L_0^2(\Omega))$ such that for each $t\in[0, T]$

and subject to the initial conditions $\mathbf{u}(\mathbf{x}, 0) = \mathbf{u}_0(\mathbf{x})$, $p(\mathbf{x}, 0) = p_0(\mathbf{x})$.

4.2.   Semi-discrete scheme

For semi-discretization in space, we use the CR-$P_1$-$P_0$ IFE space $\mathbf{S}_h(\mathcal{T}_h)$ to approximate to approximate $[H_0^1(\Omega)]^2 \times L^2(\Omega)$. We write the vector-valued IFE space $\mathbf{S}_h(\mathcal{T}_h) = U_{1h}\times U_{2h}\times W_h$. Then we propose the semi-discrete scheme as follows.

Semi-discrete IFE Scheme: Find $(\mathbf{u}_h, p_h): = (u_{1h}, u_{2h}, p_h)\in H^1(0, T;U_{1h})\times H^1(0, T;U_{2h}) \times L^2(0, T;W_h)$ such that

and subject to the initial conditions

where $\mathbf{u}_{0, h}$ and $p_{0, h}$ are some approximations (e.g. the interpolation) of $\mathbf{u}_{0}$ and $p_{0}$ in $U_{1h}\times U_{2h}$ and $W_h$. We rewrite the semi-discrete scheme in the following matrix form.

Matrix Form: Find the vector function $\mathbf{U}(t)$ such that

where $M(t)$ and $A(t)$ denote the IFE mass and stiffness matrices, and $\mathbf{F}(t)$ is the vector corresponding to the right-hand side of (36)-(37). The initial vector $\mathbf{U}^0$ takes the values of the coefficients of the interpolation $\mathbf{I}_{h} (\mathbf{u}_0, p_0)$. More details will be given in Section 5.

Remark 3. Since the interface $\Gamma(t)$ is a function of time $t$, the IFE spaces $\mathbf{S}_h(\mathcal{T}_h) = U_{1h}\times U_{2h}\times W_h$ depend on the interface location; hence they are time-dependent. Although the background mesh $\mathcal{T}_h$ is time-independent, the collections of interface elements $\mathcal{T}_h^{i(t)}$ and non-interface elements $\mathcal{T}_h^{n(t)}$ vary by time. That is why the mass matrix $M(t)$ and stiffness matrix $A(t)$ are both time-dependent.

4.3.   Full-discrete scheme

Let $0 = t_0<t_1<\cdots<t_{N-1} <t_{N} = T$ be a partition of the time interval $[0, T]$ with the uniform step size $\tau$, i.e., $\tau = T/N$, and $t_n = n\tau$. Evaluating (39) at $t = t_{n+\theta}: = t_n+\theta\Delta t$, we have

Using the following finite-difference approximations in (41)

we can obtain the following full-discrete IFE scheme.

Full-discrete IFE Scheme: Given initial vector $\mathbf{U}^{0}$, find $\mathbf{U}^{n+1}$ for each $n = 0, 1, \cdots, N-1$ in

Note that when $\theta = 1$, the method becomes the Backward-Euler method:

When $\theta = \frac{1}{2}$, the method is the Crank-Nicolson method:

Remark 4. For the time-dependent Stokes interface problem with a stationary interface, i.e. $\Gamma$ is time-independent, the matrices $M$ and $A$ in the full-discrete scheme (45) will remain unchanged as time evolves. As a result, at each time level, only the vector $\mathbf{F}^n$ needs to be updated.

Remark 5. For the time-dependent Stokes interface problem with a moving interface, although the matrices $M^n$ and $A^n$ depend on the location of interface, which further depends on time, these matrices can be efficiently generated by locally modifying the matrices from the previous time step. A unique feature of IFEM is that not only is the computational mesh interface independent, but the number as well as the location of the unknowns also remain unchanged. In two consecutive time steps, only a small portion of elements change their interface configurations, as shown in Figure 5 marked in dark yellow color. Consequently, we only need to modify local stiffness and mass matrices on those elements. The majority of the global matrices remain unchanged. This feature is also important in the error analysis of IFE methods for moving interface problems, see [10].

5.   Numerical examples

In this section, we report some numerical experiments for the mixed conforming-nonconforming IFE methods for the Stokes interface problems. We test both the interpolation and the IFE solution with various configurations of the interface and coefficient jumps. All of our numerical experiments are performed on a family of Cartesian triangular meshes which are obtained by first partitioning the domain into $N_s\times N_s$ congruent rectangles, and then further dividing each rectangle into two triangles by its diagonal with the positive slope.

We investigate the approximation property of IFE space by the interpolation. Define the CR-$P_1$-$P_0$ IFE interpolation operator is defined to be $\mathbf{I}_h: H^1(\Omega)\times C(\Omega)\times L^2(\Omega) \to \mathbf{S}_h(\mathcal{T}_h)$ such that

where $\mathit{\boldsymbol{\phi}}_{j, T}$ and $\mathit{\boldsymbol{\psi}}_{j, T}$ are the local IFE/FE shape functions given in (14) and (11), respectively. For a fixed $t$, the coefficients $c_j$ take the values

and

where $A_j$ and $e_j$, $j = 1, 2, 3$ are the vertices and edges of the the triangle $T$, respectively. The $P_1$-CR-$P_0$ interpolation can be defined similarly. The errors of the IFE interpolations are measured in $L^2$ and semi-$H^1$ norms as follows

where $u_{1, I}$, $u_{2, I}$, $p_{I}$ are components of the vector-valued function $\mathbf{I}_{h} (\mathbf{u}, p)$. In the tables below, we report the convergence rate based on two consecutive meshes $\mathcal{T}_h$ and $\mathcal{T}_{h/2}$, as well as the overall convergence rate among all meshes using the linear regression.

Example 5.1 (Interpolation Accuracy). In this example, we test the approximation capability of the new vector valued IFE space using interpolation. Since the interpolation is a time-independent procedure, we use a steady-state solution given in [1,20] for this experiment. Let the domain be $\Omega = [-1, 1]^2$ and the interface be $\Gamma = \{(x, y): x^2+y^2 = 0.3\}$. The circular interface separates the domain $\Omega$ into two subdomains $\Omega^- = \{(x, y):x^2+y^2<0.3\}$ and $\Omega^+ = \{(x, y):x^2+y^2>0.3\}$. The exact solutions $u_1$, $u_2$ and ${p}$ are defined as follows:

We first test a moderate coefficient contrast with $\mu^- = 1$ and $\mu^+ = 10$. Tables 1-2 report the interpolation errors using CR-$P_1$-$P_0$ and $P_1$-CR-$P_0$ IFE functions, respectively. We can see from these tables that the accuracy of these two IFE spaces are similar. Both of these interpolation errors obey

where $i = 1, 2$.

Next, we test a larger coefficient jump ($\mu^- = 1$ and $\mu^+ = 200$) and a flipped coefficients case ($\mu^- = 10$ and $\mu^+ = 1$). We only report the $P_1$-CR-$P_0$ IFE interpolation, since the CR-$P_1$-$P_0$ IFE results are close. Errors for large jump and flipped jump cases are listed in Tables 3-4, respectively. The convergence rates are again consistent with (50).

Example 5.2 (Unsteady Stokes Equation with Fixed Interface). In this example, we consider a time-dependent Stokes equation with a fixed interface. The domain $\Omega$ and interface $\Gamma$ is the same as in Example 5.1. The time domain is set to be $[0, 1]$, and it is partitioned uniformly to $N_t$ subintervals. We use both backward-Euler and Crank-Nicolson schemes with the time step size $\tau = 2h$. The errors are measured at the final time $t = 1$. The initial data $\mathbf{u}_0$, $p_0$, the boundary condition, and the source term $\mathbf{f}$ are chosen so that the exact solutions of this problem are as follows

Table 5 and Table 6 report the backward-Euler and the Crank-Nicolson IFE solutions at the final time $t = 1$, respectively. The numerical results indicate that the errors of Crank-Nicolson are a little smaller than those of backward-Euler. They obey the expected convergence rates

where $i = 1, 2$, and $k = 1$ for backward-Euler, and $k = 2$ for Crank-Nicolson. As before, we report only the $P_1$-CR-$P_0$ IFE solutions, and the results for the CR-$P_1$-$P_0$ IFE solution are similar.

Example 5.3 (Unsteady Stokes Equation: Circular Moving Interface). In this example we test our mixed IFE method on a Stokes moving interface problem. The interface curve is a circle centered at origin with a varying radius. The function for the interface curve is given as

It can be seen that at time $t = 0$, the interface is the same as Example 5.2 with a radius of $r = 0.54772$. As the time $t$ increases, the radius will first increase, then decrease, and finally return to the original one. The maximum and minimum radius $r_{max} = 0.67082$ and $r_{min} = 0.3873$ occur at $t = 0.25$ and $t = 0.75$, as shown in Figure 6. The exact solution $\mathbf{u}$ is written in terms of the level-set interface function $\Gamma = 0$:

In this experiment, we set the time step size $\tau = h$. We first test the moderate jump case for this moving interface problem. Table 7 reports the errors at the final time level of the backward-Euler IFE solutions. The error decay is observed to converge in an optimal order, as stated in (52). Figure 6 shows the IFE solution $u_1$ and $u_2$ at time $t = 0.25$, $t = 0.75$, and $t = 1$, respectively, on the $64\times 64$ mesh. For a larger jump case, the errors are reported in Table 8.

The condition numbers of the IFE systems are reported in Tables 9 and 10. We monitor the condition numbers at $t = 0.25$, $t = 0.75$, and $t = 1$ which correspond to the interface circle listed in Figure 6. We test different contrast ratios by fixing the coefficient $\mu^- = 1$ and varying the other coefficient $\mu^+ = 0.01$, $0.1$, $1$, $10$, and $100$. Note that when $\mu^+ = 1$, there is no jump in coefficient, hence the IFE scheme becomes the standard FE scheme. Even in this no-jump case, we observe that the condition number is of order $\mathcal{O}(h^{-4})$. We also observe that the condition number increases as the jump ratio enlarges. No significant differences have been noticed for backward-Euler and Crank-Nicolson in terms of the conditions numbers.

6.   Conclusion

In this paper, we developed a mixed conforming-nonconforming immersed finite element method for unsteady Stokes interface problems. The proposed vector-valued IFE spaces use conforming $P_1$ approximation for one velocity component and nonconforming $P_1$ approximation for the other. The pressure is approximated by piecewise constant. Unisolvency of the vector-valued IFE functions is proved. Interpolation errors are observed to be optimal which indicates these new IFE spaces have sufficient approximation capabilities. This new IFE function space can be used to solve steady-state, and unsteady Stokes interface problems. In addition, we have extended the application to a moving interface problem, and our numerical results show the optimal convergence in some benchmark tests. The proposed IFE method can be extended to some more general fluid flow interface problems such as Navier-Stokes equations, which will be an interesting future topic to explore.