Solving the Incompressible Surface Stokes Equation by Standard Velocity-Correction Projection Methods

In this paper, an effective numerical algorithm for the Stokes equation of a curved surface is presented and analyzed. The velocity field was decoupled from the pressure by the standard velocity correction projection method, and the penalty term was introduced to make the velocity satisfy the tangential condition. The first-order backward Euler scheme and second-order BDF scheme are used to discretize the time separately, and the stability of the two schemes is analyzed. The mixed finite element pair (P2,P1) is applied to discretization of space. Finally, numerical examples are given to verify the accuracy and effectiveness of the proposed method.


Introduction
Fluid equations on surfaces are used in mathematical modeling of emulsions [1], foams [2], and biofilms [3] where curved fluid membranes behave transversely as efficient viscoelastic media, as studied by Rahimi et al. [4]. Typically, such models consist of surface (Navier)-Stokes equations. These equations are also studied as an interesting mathematical problem in their own right. At present, the discretization methods of fluid equation on a curved surface are mainly the finite element method of a scalar elliptic equation and a parabolic partial differential Equation [4,5].
Due to the importance of surface incompressible flow in practical applications, many scholars have carried out in-depth research on it. Many achievements have also been made in the study of discretization methods of vector equations. Nitschke, Reuther et al. proposed the finite element method of Navier-Stokes equations of curved surfaces based on curl operators [6,7]. Reuther and Fries discretized the incompressible Navier-Stokes equations using a curved finite element method [8,9]. Reuther considered the P 1 − P 1 finite element that does not satisfy inf-sup condition, and used the punishment technique to force the velocity field to be tangent to the surface. In contrast, Fries used the P 2 − P 1 element, combined with the Lagrange multiplier method, to force the velocity to meet the tangential condition. However, the numerical analysis of the Navier-Stokes finite element method for curved surfaces was not involved in these two papers. Later, scholars began to use the finite element method to analyze the vector Laplace equation on the surface. This was the first step in extending the analysis of scalar problems to surface Navier-Stokes equations. They cross-analyzed a curved finite element method combined with a penalty technique to impose tangential constraints, and gave the results of stability analysis and error analysis, and also explained the influence of geometric errors on solving the problem [10]. A nonfitting finite element method, the trace finite element method, was proposed by Hansbo. When combining it with the Lagrange multiplier method, the forced velocity field satisfies the tangential condition, and the stability and optimal error estimation of the method are also proved [11].
The structure of this paper is as follows: The definition of surface differential operator is introduced in Section 2. In Section 3, the theory of surface finite element approximation and the Stokes equation are introduced. In Section 4, the standard format of first-order and second-order velocity correction methods is introduced and the stability analysis of semidiscrete time is proved. In Section 5, numerical examples are given to verify the accuracy and validity of the proposed method. The sixth part is the conclusion and prospects.

Surface Differential Operator
In this paper, we consider that the surface Γ ⊂ R 3 is closed and sufficiently smooth, and can be defined by a level set function satisfying φ(x) = 0, so that The unit normal vector of the surface Γ at node x can be defined as On the surface Γ, the orthogonal projection P(x) ∈ R 3 is defined by the normal vector Define the neighborhood of the surface Γ as In addition, for the symbol distance function defined on the surfaced(x) = dist(x, Γ), The normal vector n is continuously extended along the normal direction in the neighborhood, and can be obtained by n(x) = n(a(x)) = ∇d(x).
From the above property, the projection P : U δ → Γ of the nearest point in the neighborhood U δ of the surface Γ is clearly defined. Thus, for a scalar function u(x), differentiable on the surface Γ, the tangential gradient operator ∇ Γ is defined as Then, the gradient operator of a vector u(x) = (u(x), v(x), w(x)) T is defined as follows: Definition of the strain tensor of surface [12]: Similarly, for any vector u and tensor A, the surface divergence operator is defined as follows: div Γ u := tr(∇ Γ u) = tr(P(∇u)P) = tr((∇u)P) = ∇ Γ × u, In this section, we review some surface differential operators of smooth functions and derive some properties of these operators. Therefore, tangent vector function space is introduced: Then, the curl operator can be defined as In the Gaussian curvature κ(x) on the surface of Γ, for all the ψ ∈ C 1 (Γ), u ∈ C 1 t (Γ) 3 , with the following formula: In the literature [13,14], the following equation is found: When u ∈ C 2 t (Γ) 3 satisfies div Γ u = 0, then we have where: κ is Gaussian curvature.

Remark 1.
An internal form of the surface incompressible Stokes equation is independent of the selected coordinate system. Compared with the plane incompressible Stokes equation, not only is the corresponding operator replaced by the surface operator, but also, the Gaussian curvature has an additional contribution. This is because the surface strain tensor E s (u) = 1 2 (∇ Γ u +(∇ Γ u) T ) divergence of the cause uses the Codazzi-Mainardi equation and the incompressible condition

Finite Element Approximation of Surface
In this section, the standard Hilbert space and some of the symbols reused in the following sections are introduced.
There are a lot of studies about the numerical calculation of the two-dimensional Stokes equation in plane space [15][16][17], but there are few corresponding results on surfaces. The unsteady incompressible Stokes Equation [18] of a surface is shown below.
  where: Re denotes the surface Reynolds number, u ∈ V denotes the surface velocity, p ∈ Q denotes the surface pressure, f ∈ L 2 (Γ) 3 denotes source term, and f × n = 0. Pressure is constrained in the following form: In the Cartesian coordinate system, the curved Stokes equation is further modified by rotating the velocity field to reduce the complexity of the equation, and the original equation can be rewritten as   This formulation ensures the velocity to be tangential only weakly through the added penalty term and is equivalent to Equations (3) and (5) only if u × n = 0 [19].
For weak formats of Stokes problems, spaces with norms are used to define The corresponding space of the tangent vector field is defined as For u ∈ V, the velocity field u is decomposed into tangential and normal components using the following symbols [20]: Below, the time-dependent sequence's finite difference is denoted by the following symbol [21]: {φ}.
Finally, the following three identities are often used in theoretical analysis:

First-Order Velocity Correction Method
Firstly, the viscous velocity term is explicitly dealt with in the first step, and then it is modified in the second step. The corresponding scheme is as follows: u 0 = u 0 , and we select u 1 to better approximate u(∆t), when k ≥ 1; solve (u k+1 , p k+1 ; u k+1 ). ∆t is the time step and has the following first-order format.
Step 1. We solve for u k+1 ∈ V T , p k+1 ∈ Q from Step 2. We solve for u k+1 ∈ V T from

Implementation of the standard form
It is difficult to solve Equation (12) as a weak Poisson problem for pressure due to the existence of the second derivative of Equation (12). To avoid this difficulty, the algorithm can be rewritten into an equivalent form by algebraic substitution.
By subtracting (12) of k step format from (12) of k + 1 step format and by substituting (13) of k step format into the resulting equation, a more adequate form of the projection step (14) is obtained. The actual implementation steps are as follows.
Step 1. We solve for u k+1 ∈ V T , p k+1 ∈ Q from Equation (14) can be rewritten as Step 2. We solve for u k+1 ∈ V T from Note that in algorithm (14)- (16), projection velocity u k+1 be completely eliminated. Therefore, in the actual calculation, consider using approximate velocity u k+1 instead of u k+1 .

Second-Order Velocity Correction Method
Similarly, the time semi-discrete scheme of the second-order standard velocity correction method is given below. u 0 = u 0 and select u 1 to better approximate u(∆t); when k ≥ 1, solve (u k+1 , p k+1 ; u k+1 ); ∆t is the time step with the following second-order format.
Step 1. We solve for u k+1 ∈ V T , p k+1 ∈ Q from Step 2. We solve for u k+1 ∈ V T from 3 2∆t

Implementation of the standard form
It is difficult to solve Equation (17) as a weak Poisson problem for pressure due to the existence of the second derivative of the equation. To avoid this difficulty, the algorithm can be rewritten into an equivalent form by algebraic substitution.
Through (17) of k + 1, step format minus k, step format, and combining (18) of k step format, we get a new step projection of (19) type. The actual implementation steps are as follows.
Step 1. We solve for u k+1 Equation (19) can be rewritten as Step 2. We solve u k+1 ∈ V T from 1 2∆t Note once again that projection velocity u k+1 be completely eliminated from the algorithm (19)- (21). Therefore, in the actual calculation consider using approximate velocity u k+1 instead of u k+1 .

Stability Analysis
In this section, the stability analysis of the first and second-order time semi-discrete schemes of the standard velocity correction projection method is established. Without loss of generality, we assume f ≡ 0.

First-Order Scheme Stability Analysis
A key step in establishing the stability analysis results is to reformulate the standard velocity correction scheme using the Gauge-Uzawa formula. To be more precise, we introduce the auxiliary variable w k and define it as follows.

Stability Analysis of Second-Order Schemes
As with the first-order proof stability analysis method, an auxiliary variable w k is introduced and defined as follows: We can rewrite (18) equation as

Theorem 2. Equations (17) and (18) represent unconditional energy stability. For all
where: is the correction energy of the time step k.
Proof. Take the inner product of (17) with 4∆tu k+1 , to get where: We deal with this term by using a similar treatment as in [22]. It can be rewritten I as Using Equations (7) and (10), we can rewrite I 1 as Thanks to Equations (29) and (30), we can rewrite I 2 as By Equations (6), (7) and (11), we can obtain Equation (9) tells us that Next, taking the inner product of both sides of Equation (30) with itself, we get Equation (37) can also be written as Equation (38) Summing up (32) and (38), and taking into account (33), (35) and (36), we can obtain

Fully Discrete Scheme
The first-order backward Euler and second-order BDF are used in discrete-time, and the mixed finite element pair is used in discrete space. Γ h = {K} is Γ on a consistent regular triangular mesh, and the grid size is h = max K∈Γ h {diam{K}}. Define the following discrete subspace (V Th , Q h ) ⊂ (V T , Q):

First-Order Fully Discrete Velocity Correction Method
Step 1.
Step 2. For ∀v ih ∈ V T , we solve for u k+1

Second-Order Fully Discrete Velocity Correction Method
Step 1. For ∀v ih ∈ V Th , ∀q h ∈ Q h , we solve for (u k+1 h , p k+1 h ) from Step 2. For ∀v ih ∈ V T , we solve for u k+1

Numerical Experiments
In this section, some numerical examples are designed to verify the effectiveness and accuracy of the velocity correction method; test the convergence and stability; and simulate the physical phenomena of circulating flow in a single fluid system.

Convergence Test
In order to verify the accuracy and effectiveness of this method, two different surfaces φ 1 and φ 2 were selected for a convergence test.
The spherical implicit function: The ring implicit function: The exact solution of the given surface Stokes problem is as follows.
For the first-order numerical discrete scheme, let Re = 10, t = 0.1, α = 1000/h 2 , time step for ∆t = h 2 , h = 0.25, 0.15, 0.125, 0.06, 0.03. The order of convergence is calculated by the following formula [23]: (43) Figure 1 shows convergence on a sphere and a ring with different surfaces of the first-order scheme. The results show that pressure p is in the L 2 norm of convergence in order to achieve second-order convergence, and velocity u is in the L 2 norm under second-order convergence. This is because of the existence of spacial geometric error.

A Circular Flow on a Ring
In this experiment, the unsteady curvature surface φ 2 is used to simulate the physical phenomenon of the flow velocity change of the fluid on the ring [24,25].
Set the initial value u 0 as follows.
Take Re =10, h = 0.3, α = 3000, and ∆t = 0.1 ( Figure 2) to simulate the single fluid system. The initial velocity u 0 will gradually evolve into periodic flow over time. When the penalty coefficient α = 0, the velocity field is not tangent to the surface, and periodic flow cannot be generated. Simulation results show that the proposed method can effectively simulate the physical phenomena of fluid flow and produce the expected intra-ring circulating flow.

Stability Test
The system is studied when the external force f = 0. Consider a sphere with constant curvature φ 1 (x, y, z) = x 2 + y 2 + z 2 − 1.

Conclusions
In the curved Stokes fluid model, the velocity and pressure are decoupled by the standard velocity correction method, and the equivalent elliptic partial differential equation is obtained. The first-order and second-order numerical schemes of (P 2 , P 1 ) are constructed by using curved surface mixed finite element matching for space discretization. The stability analysis of first-order and second-order time semi-discretization is established. Finally, three examples are given to verify the accuracy and effectiveness of the proposed method.