Distributed Lagrange Multiplier/Fictitious Domain Finite Element Method for a Transient Stokes Interface Problem with Jump Coefficients

The distributed Lagrange multiplier/fictitious domain (DLM/FD)-mixed finite element method is developed and analyzed in this paper for a transient Stokes interface problem with jump coefficients. The semiand fully discrete DLM/FD-mixed finite element scheme are developed for the first time for this problem with a moving interface, where the arbitrary Lagrangian-Eulerian (ALE) technique is employed to deal with the moving and immersed subdomain. Stability and optimal convergence properties are obtained for both schemes. Numerical experiments are carried out for different scenarios of jump coefficients, and all theoretical results are validated.

The regularity study for solutions to (1)- (10) is still open to the community of theoretical partial differential equations, especially when the interface Γ t deforms along the time, i.e., the shape of Γ t and Ω i t (i = 1, 2) depend on the primary unknowns (u, p). In this paper, in order to show a certain amount of convergence rate during the numerical experiments process for validating the convergence theorem of the developed DLM/FD finite element method, in what follows we assume a reduced regularity result for the solution (u, p) to (1)-(10) which is similar with that of the stationary Stokes interface problem in space [Shibataa and Shimizu (2003); Hansbo, Larson and Zahedi (2014); Olshanskii and Reusken (2006)], where, 3/2 < σ ≤ 2. Without loss of generality, in this paper we only study the immersed interface case by assuming Ω 2 t ⊂ Ω, as shown in the left of Fig. 1. The regularity property (11) is assumed to hold under the circumstance that Γ t does not deform but only rotates and/or translates with a prescribed domain velocity w(x, t), as defined in Section 2.3. Thus the shape and position of Ω i (i = 1, 2) are prescribed and do not depend on the primary unknowns, and the regularity results (11) can still be hypothesized, accordingly. Numerical results shown in Section 5 also support the regularity property of solution (u, p) defined in (11). In practice, (1)-(10) generally model a type of immiscible two-phase fluid flow problem, where two phases of the fluid are separated by an distinct interface, and both fluid phases are defined by Stokes/Navier-Stokes equations in terms of fluid velocity and pressure as sketched in (1)-(4). In this scenario, β i and ρ i (i = 1, 2) may stand for the fluid viscosity and density of different phases. Hence, the essential characteristic of the immiscible twophase fluid flow model is preserved in the transient Stokes interface problem (1)-(10), that is, two different types of fluid equations bearing with different viscosity and density are defined on either side of the moving interface Γ t . Some long existing body-unfitted mesh methods for interface problems such as the immersed interface method (IIM) [Deng, Ito and Li (2003); LeVeque and Li (1994); Li and Ito (2001)] and immersed finite element method (IFEM) [Li (1998) ;Ji, Chen and Li (2014)] are still far from satisfactory for solving the Stokes interface problem in either stationary or transient case. As for the representative body-fitted mesh method, the arbitrary Lagrangian-Eulerian (ALE) method [Hirth, Amsden and Cook (1974) ;Hughes, Liu and Zimmermann (1981); Huerta and Liu (1988); Nitikitpaiboon and Bathe (1993); Souli and Benson (2010)] is the most popular one for solving moving interface problems such fluidstructure interactions (FSI), where, the mesh on the interface is accommodated to be shared by both fluid and structure, and thus to automatically satisfy the interface conditions as sketched in (5) and (6). However, for large rotations and/or translations of the structure or inhomogeneous movements of the grid nodes, fluid elements tend to become ill-shaped, which reflects on the accuracy of the solution. In this case, re-meshing, in which the whole domain or part of the domain is spatially rediscretised, is then a common strategy. However, it could be very troublesome, time consuming and less accurate, and, the worst thing brought by the re-meshing is that the mesh connectivity is no longer preserved for ALE method and thus many properties of ALE method are lost. To overcome the above problems and to deliver an efficient and accurate numerical method for the transient Stokes interface problem in which the immersed phase may be engaged in a large translational/rotational motion, in this paper we develop a body-unfitted mesh method based upon the framework of the distributed Lagrange multiplier/fictitious domain (DLM/FD) method [Glowinski, Pana, Hesla et al. (1999); Wachs (2007); Glowinski and Kuznetsov (2007); Boffi and Gastaldi (2017); Wang and Sun (2017)], where, one fluid phase is smoothly extended into the other phase that is defined in the immersed subdomain, then occupies the entire domain Ω, and the Lagrange multiplier (physically a pseudo body force) is introduced to enforce the interior (fictitious) fluids in the immersed subdomain to satisfy the constraint of the immersed phase motion. The constraints are incorporated into the field equations to form an augmented matrix equation which involves the Lagrange multipliers as unknowns. Thus, the re-meshing in the fluid domain is no longer needed for DLM/FD method, and the possible failure of ALE method is completely avoided when the large translation/rotation occurs to the immersed phase motion [Auricchio, Boffi, Gastaldi et al. (2015); Shi and Phan-Thien (2005); Yu (2005); Glowinski, Pana, Hesla et al. (2001)]. The DLM/FD finite element method has been analyzed for the elliptic interface problem [Boffi, Gastaldi and Ruggeri (2014); Auricchio, Boffi, Gastaldi et al. (2015)], the parabolic interface problem [Wang and Sun (2017)], the stationary Stokes interface problem [Lundberg, Sun and Wang (2019)], but has not yet applied to the transient Stokes interface problem. As shown in Boffi et al. [Boffi, Gastaldi and Ruggeri (2014); Auricchio, Boffi, Gastaldi et al. (2015)], the DLM/FD method essentially produces a saddle-point problem in regard to the unknown of elliptic equation and Lagrange multiplier, so the existing Babuska-Brezzi's theory [Babuška (1971); Brezzi and Fortin (1991); Brezzi (1974); Brezzi and Pitkaranta (1984)] can be employed to analyze the well-posedness, stability and convergence properties of the corresponding saddle-point problem induced from the DLM/FD finite element method. However, for the stationary Stokes interface problem, which is the steady state of the transient Stokes interface problem (1)-(10), we can see that its corresponding DLM/FD formulation forms a nested saddle-point problem including two subproblems of saddle-point type: the inside one from Stokes equations regarding Stokes unknowns (velocity and pressure), and the outside one from the DLM/FD method itself regarding Lagrange multiplier and Stokes unknowns, of which the well-posedness, stability as well as convergence analyses are more sophisticated than those of the elliptic and the parabolic interface problems. In the authors' recent work [Lundberg, Sun and Wang (2019)], a modified DLM/FD finite element method is developed for a stationary Stokes interface problem that consists of a nested saddle-point problem, and its well-posedness, stability and optimal convergence properties are analyzed still by means of the Babuska-Brezzi's theory but a more complicated approach. So in this paper, we will be able to develop the DLM/FD finite element method for the transient Stokes interface problem (1)-(10) and analyze its stability and convergence properties based on our previous work.
The structure of the paper is the following: in Section 2 we introduce the fictitious fluid (Stokes) equations then derive weak formulations of a transient Stokes interface problem with and without the employment of DLM/FD method. Then we define the semi-discrete DLM/FD finite element approximation and analyze its stability and optimal convergence theorem in Section 3. The full discretization is defined and its stability and convergence properties are analyzed in Section 4. Numerical experiments are carried out in Section 5, where the theoretical convergence results are validated.
such that 2.2 Weak formulations.
If we add the fictitious fluid Eqs. (13)- (14), which are defined in Ω 2 t , to the Stokes Eqs. (1)-(2), which are defined in Ω 1 t , and integrate by parts, then On the other hand, if we subtract the fictitious fluid Eqs. (13) and (14) from the Stokes Eqs.
(3) and (4), and integrate by parts in Ω 2 t , then If we add (21) to (23) and (22) to (24), then the terms of the fictitious Stokes equations and the normal derivative terms on Γ t are all cancelled out, resulting in the original weak formulation of (1)-(10) as follows.
In the following theorem, we prove the equivalence between Weak Forms I and II. (25) and (26), and define λ ∈ Λ t by whereũ We have (44) due to (48) , then (25) can be easily proved by taking v ∈ V in (45) with v| Ω 2 t = v 2 , and simply adding (45) and (47) together to cancel all Lagrange multiplier terms. Due to (33) and (46), (26) is obvious.

The arbitrary Lagrange-Eulerian (ALE) formulation
We assume that there exists a bijective mapping X t ∈ H 1 (0, T ; (W 2,∞ (Ω 2 0 )) d ) such that for each t ∈ (0, T ], the mapping [Martín, Smaranda and Takahashi (2009) is invertible and X −1 Here y ∈ Ω 2 0 is so called the arbitrary Lagrange-Eulerian (ALE) coordinate, and x ∈ Ω 2 t is the spatial (or Eulerian) coordinate. We further introduce the domain velocity w, defined by In this paper we assume that Ω 2 t does not deform but only rotates/translates with a prescribed domain velocity, w.
is a given velocity function with the assumption that where c denotes a constant independent of any discretization parameters in the rest of the paper. so the ALE mapping function X t (y, t) ∈ H 1 (0, T ; (W 2,∞ (Ω 2 0 )) d ), which can be considered as a prescribed displacement of domain motion for Ω 2 t .
We use dv dt | y to denote the temporal derivative on the ALE frame which is defined as follows: for any function v : Ω 2 t → R d regular enough and defined on the Eulerian frame, we set [Gastaldi (2001); Martín, Smaranda and Takahashi (2009) Now we need to redefine the space V 2 t on the ALE frame as follows: t . Based on the above definitions, we can reformulate (51)-(54) as the following ALE-type weak formulation of the DLM/FD method.
(ψ,ũ| Ω 2 is the following where, V h × Q h is a stable pair of P 2 P 1 mixed (Taylor-Hood) finite element space for Stokes equations. In fact, it is proved in [Lundberg, Sun and Wang (2019)] that such chosen mixed finite element spaces for a fixed time t is stable for the developed DLM/FD method for the stationary Stokes interface problem -the steady state case of the transient Stokes interface problem (1)-(10). Thus, based upon the DLM/FD weak form III (51) H,t is still adopted as the stable mixed finite element spaces for the DLM/FD finite element method of (1)-(10), as defined below.
Then, we have the following convergence theorem for the semi-discretization (64)-(67). Theorem 3.1. Suppose (ũ, u 2 , p, φ) is the solution to (59)-(62), (u h , u 2,H , p h , φ H ) is the solution to (64)-(67). With P 2 -P 2 -P 1 -P 2 mixed finite element to respectively discretize u h , u 2,H , p h , φ H , we have the following error estimate Note that on the right hand side of (74), all terms but inf can be directly estimated based on a priori interpolation error estimates and the regularity assumptions (18). To find an error estimate for inf Integrating by parts gives Similarly, we can pick (v, v 2 ) ∈ V × V 2 t such that v| Ω 2 t = v 2 and v = 0 outside Ω 2 t including ∂Ω 2 t . Add (51) to (53) and integrate by parts, yields Now, let v ∈ V and take v 2 = v| Ω 2 t in Ω 2 t . Add (51) to (53), we have then apply (75) and (76), we attain Further, apply (76) and (79) to (53), and note that n 1 = −n 2 , we have With the regularity assumptions (18), we can obtain the following estimates by doing an analogous analysis with [Auricchio, Boffi, Gastaldi et al. (2015); Lundberg, Sun and Wang To get an estimate for inf (28) and (30), we just need to find an error estimate for inf To that end, we first let π H be the L 2 projection of V 2 t into V 2 H,t , that is, for any It is easy to see that λ 4 ∈ L 2 (Ω 2 t ) and λ 4 0 ≤ β β2 ρ 2 −ρ 2 ∂u2 ∂t 0 . On the other hand, because of the choice of finite element space (63), our finite elements V 2 H,t and Λ H,t are contained in L 2 (Ω 2 t ), we can interpret the duality pairing as scalar product in L 2 (Ω 2 t ) [Auricchio, Boffi, Gastaldi et al. (2015); Boffi, Gastaldi and Ruggeri (2014)]. Thus, we can define P H λ 4 ∈ Λ H,t be the L 2 -projection of λ 4 onto Λ H,t such that So by (87) where, we apply (87) and (83). By applying the Cauchy-Schwartz inequality and the a priori interpolation error estimate for π H , we obtain Then there exists a constant c > 0 such that where, the trace estimate is applied to get u 2,H 0,Γt ≤ c u 2,H V 2 t . Integrate both sides of (93) in time from 0 to t, add u 2,H 2 0,Ω 2 t to both sides, choose a sufficiently small ε, and apply Poincaré inequality in Ω and Grönwall's inequality to (93), reads Then, we have the following stability theorem for the semi-discrete scheme.

Numerical experiments
In this section, we study the numerical performance of the developed DLM/FD finite element method for an example of the transient Stokes interface problem (1)-(10) defined in Ω = [0, 1] × [0, 1], where the circular subdomain Ω 2 t makes a translational motion and the position of ∂Ω 2 t , which is the interface Γ t , satisfies where, w = (w 1 , w 2 ) T denotes the moving velocity of Γ t . We properly choose the functions of coefficients, source terms and jump flux of (1)-(10), i.e., β i , ρ i , f i , (i = 1, 2) and τ such that the true solution (u, p) to (1)-(10), where u = (u, v) T , is defined by where, β = β i (x), ∀x ∈ Ω i (i = 1, 2) is chosen as a piecewise constant depending on the location of x. Clearly, such chosen solution (u, p) satisfies the following regularity property: In what follows, we take a constant moving velocity w = (0.1, 0.2) T , and let T = 1. The meshes T h (Ω) and T H (Ω 2 t ) are constructed independently and thus mismatched with each other. Convergence results of the velocity vector at the time t = T in its H 1 -and L 2 norm, i.e., u − u h (H 1 (Ω)) 2 and u − u h (L 2 (Ω)) 2 which are displayed in their component forms, and of the pressure in its L 2 norm, p−p h L 2 (Ω 1 T ) , are illustrated in Tabs. 1 and 2 for large jump coefficient cases. We can observe that: (1) the developed DLM/FD-mixed finite element discretization is stable and converges in all cases, little influence from the choice of the time step size; (2) the convergence results are relatively more sensitive to β 2 /β 1 , comparing with the jump ρ 2 /ρ 1 , noting that the exact solution u depends on β, but independent of ρ; (3) due to the reduced regularity property of the solution, and the discontinuity of the normal derivative of u across Γ t , the convergence rates of velocity errors in H 1 -and L 2 -norm decrease to 0.55 ∼ 0.9 and 1.0 ∼ 1.3, respectively, and the convergence rates of pressure errors in L 2 -norm keeps around 1.0 ∼ 2.0, which validate our theoretical conclusions, and also match with the convergence rates of other types of interface problems when the DLM/FD method is applied [Boffi, Gastaldi and Ruggeri (2014); Auricchio, Boffi, Gastaldi et al. (2015); Wang and Sun (2017); Lundberg, Sun and Wang (2019); Sun (2019)]. Next, we investigate the influence of time step size on the convergence rate of the developed DLM/FD finite element method. In order to let O(∆t) be the main part of the error in comparison with the part O(h σ−1 + H σ−1 ), we particularly pick up the case of β 2 /β 1 = 2, ρ 2 /ρ 1 = 2, and take f i , (i = 1, 2) and τ such that the true solution (u, p) to (1)-(10), where u = (u, v) T , is defined by u = (y − 0.3 − w 2 t)((x − 0.3 − w 1 t) 2 + (y − 0.3 − w 2 t) 2 − 0.01)(2t 9 − t 5 )/β, v = −(x − 0.3 − w 1 t)((x − 0.3 − w 1 t) 2 + (y − 0.3 − w 2 t) 2 − 0.01)(2t 9 − t 5 )/β, p = sin(πx) sin(πy)(2t 9 − t 5 ).

Conclusion and future work
We develop the DLM/FD-mixed finite element method for a generic transient Stokes interface problem and carry out numerical analyses for both semi-and fully discrete scheme on the convergence and stability properties. By using the Taylor-Hood (P 2 P 1 ) mixed finite element space, we are able to obtain a nearly optimal convergence rate for both the velocity and the pressure in their respective norms, subjecting to the reduced regularity assumption for the solution to the transient Stokes interface problem. Numerical experiments validate the theoretical results, showing that the convergence rates of the velocity with respect to the mesh size is the 0.5th in H 1 -norm, and the first order in L 2 -norm, at least, which is true even for larger jump coefficient cases up to 1:10000, relatively insensitive to different choices of jump coefficients and time step sizes. And, the first order convergence rate with respect to the time step size is also validated.