Error Analysis of a PFEM Based on the Euler Semi-Implicit Scheme for the Unsteady MHD Equations

In this article, we mainly consider a first order penalty finite element method (PFEM) for the 2D/3D unsteady incompressible magnetohydrodynamic (MHD) equations. The penalty method applies a penalty term to relax the constraint “∇·u=0”, which allows us to transform the saddle point problem into two smaller problems to solve. The Euler semi-implicit scheme is based on a first order backward difference formula for time discretization and semi-implicit treatments for nonlinear terms. It is worth mentioning that the error estimates of the fully discrete PFEM are rigorously derived, which depend on the penalty parameter ϵ, the time-step size τ, and the mesh size h. Finally, two numerical tests show that our scheme is effective.

In this paper, we mainly consider the 2D/3D unsteady incompressible MHD equations. This model is a coupled strongly nonlinear system, and it is a saddle point problem due to the incompressible constraint. Therefore, it is necessary to construct unconditionally stable and decoupled algorithms for our model. For the time discretization, it is well-known that simple discretizations, like fully explicit or implicit type schemes, can lead to considerable instabilities or suffer from costly time expense. Recently, Euler semi-implicit schemes for some evolution differential equations have been given in [4,13,14]. This method is unconditionally stable. For the saddle point problem, there are many methods to release the incompressibility constraint for incompressible flow such as the projection method, the pressure stabilization method, the artificial compressibility method, and the penalty method (see also [15][16][17][18][19][20]).
It is worth mentioning that the penalty method is the simplest and the most basic of these methods mentioned above. For the penalty method, it can be traced back to [21]. Then, the optimal error estimate of the unsteady NS system based on the penalty finite element method (PFEM) was given in [22]. A PFEM of a Euler implicit/explicit scheme for the unsteady NS system was proposed in [23]. An error estimate of the unsteady NS system based on the P 1 nonconforming PFEM was given in [24]. The authors study the PFEM of the steady MHD equations in [25]. A decoupling PFEM for the steady incompressible MHD equations was given in [18].
The aim of this paper is to develop a first order linear and decoupled scheme. We adopt an implicit scheme for the linear terms and semi-implicit treatments for nonlinear terms. Meanwhile, the penalty method is used for fluid equations. This method decouples the MHD equations into two small equations; one is the equations of the velocity and magnetic field (u, B), and the other is the equation of pressure p. Then, the error estimates for the developed scheme are present, which depend on the penalty parameter , the time-step size τ, and the mesh size h. Finally, we give two numerical tests to verify the theoretical results of our method.
The structure of the paper is as follows: in Section 2, we present the model and estimates for the solutions of penalty MHD equations. In Section 3, we give the Euler semi-implicit scheme in the case of time discretization and the convergence rate of the semi-discrete solutions. In Section 4, we give the PFEM based on the Euler semi-implicit scheme. In Section 5, we give the error estimates of the fully discrete solutions. In Section 6, the error estimates are obtained for our scheme. In Section 7, two numerical tests show that our scheme is effective. Finally, we give some conclusions.

Functional Setting of the Unsteady MHD Equation
In this paper, we consider the unsteady incompressible MHD equations as follows: Here, u, p, B are the velocity, the pressure, and the magnetic field, f is the external force term, g is the known applied current with ∇ · g = 0, and n denotes the outward normal on ∂D. For the physical parameters, ν −1 = R e (fluid Reynolds number), µ −1 = R m (magnetic Reynolds number) and S is the coupling coefficient.
Then, we give some notations and estimates for MHD equations. For 1 ≤ r ≤ ∞, L r (D) denotes the usual Lebesgue space on D with the norm · L r . The inner product of the space L 2 (D) is denoted by (·, ·) that is (u, v) = D uvdx, and the norm of the space L 2 (D) is denoted by · . For all non-negative integers k and r, W k,r (D) stands for the standard Sobolev space equipped with the standard Sobolev norm · k,r . The norm of the space W k,2 (D) is represented by · k . The functions and spaces of vectors are represented in boldface.

Remark 1.
Taking the L 2 inner product of the first equation in (1) with v, the second equation in (1) with q and summing up the two relations, we obtain the first equation of (4). Taking the L 2 inner product of the third equation in (1) with C, we obtain the second equation of (4). Equation (5) can be obtained similarly.

Remark 2.
For φ ∈ H 2 0 (D), taking C = ∇φ in the second equation of (4), we easily obtain (∇ · B) t = 0, which implies ∇ · B(t) = ∇ · B 0 = 0. In addition, the second equation of (4) has an equivalent form (cf. [1,7]) in which ∇ · B is used as a penalty term. For ∀t ∈ (0, T ), we choose C = ∇φ(t) ∈ W in the second equation of (4); here, φ(t) is generated by the boundary value problem We can obtain 1 2 Using the operators A 1 , A 2 , we can rewrite the penalized system (5) as Referring to [2,6,13,28], the following estimates hold: where c(D) > 0 is a constant, which have different values in different cases. In this paper, C > 0 denotes a constant depending on (ν, µ, S, D, T , u 0 , B 0 , f , g), which may have different values in different cases. We make the following assumptions for (1), which specify the regularity of the data and the smoothness of the domain D (cf. [1,4]). Assumption 1. The initial data u 0 ∈ X 0 ∩ H 2 (Ω) and B 0 ∈ W 0 ∩ H 2 (Ω), the external force f , and the applied current g satisfy the following bound: Assumption 1 ensures that there is a unique strong solution for (1) over some time interval [0, T ); we have (cf. [5]) such that u t , B t ∈ L 2 (0, T ; L 2 (D)), and Equation (4) holds for almost all t ∈ [0, T ). If the data u 0 , B 0 , f and g are sufficiently small, then the solution exists for any T > 0 and satisfies Assumption 2. The problem (4) has a weak solution (u(t), p(t), B(t)) satisfying u ∈ L 2 (0, T ; X 0 ), p ∈ L 2 (0, T ; M) and B ∈ L 2 (0, T ; W 0 ) such that where H −2 (D) is the dual space of H 2 (D) ∩ X, and · −2 is the corresponding norm.
For the proof of these results, we can refer to [1,4].

Theorem 2.
Under Assumptions 1-3, the solution (u (t), p (t), B (t)) of the problem (3) satisfies the estimates We can finish the proof by a similar technique used in the proof of Theorem 1.

Theorem 3.
Under Assumptions 1-3, we have the following estimate (cf. [26,29]) Proof. We first consider the linear form of MHD equations. Then, we subtract the penalized linear MHD equations from the linear MHD equations to obtain their error equations. The error estimate in linear form is obtained through its dual problem Next, we obtain the following error estimate by choosing an appropriate L 2 inner product for the error equations sup 0≤t≤T Finally, we transform the nonlinear MHD equations into an intermediate linear equations, and then obtain Theorem 3 by applying a suitable L 2 inner product to this system and using the previous result.
In addition, we need the following discrete Gronwall lemma (cf. [4,30,31]). Lemma 2. Let a n , b n , d n and C 0 be nonnegative numbers for integer n ≥ 0, such that Then

The Euler Semi-Implicit Scheme and Its Error Estimates: Time Discretization
In this section, we consider a time discretization for the penalty MHD system (6). Let τ = T N be the time-step size and N > 0 is an integer. Then, t n = nτ, n = 1, 2, · · · , N denote the discrete time levels. The time-discrete approximations to (u (t n ), p (t n ), B (t n )) will be denoted by (u n , p n , B n ) for all 1 ≤ n ≤ N. Consider the Euler semi-implicit time-stepping algorithm: Given (u n−1 , B n−1 ) ∈ X × W, find (u n , B n ) ∈ X × W, p n ∈ M such that where (u 0 , B 0 ) = (u 0 , B 0 ), d t u n = 1 τ (u n − u n−1 ). We can rewrite (16) as Remark 5. Since ∇ · B 0 = 0, we take C = ∇φ with φ ∈ H 2 0 (D) and deduce from the second equation in (16) and the identity ∇ × ∇φ = 0 that ∇ · B n = 0 for all 0 ≤ n ≤ N.
Next, we give a priori bound of the scheme (17).
Proof. Taking the L 2 inner product of the first equation in (17) with 2u n τ, and the second equation in (17) with 2SB n τ, we obtain By using (7), (9), Lemma 1 and the Young inequality, we have Combining the above inequalities with (18), we obtain Summing the above inequality from 1 to m, we derive for all 1 ≤ m ≤ N. Using Assumption 1 and (19), we obtain Theorem 4.
Next, we establish the error estimates in time for the Euler semi-implicit scheme (17). To do this, subtracting (17) from (6) and setting e n u = u (t n ) − u n , e n B = B (t n ) − B n , we have d t e n u + νA 1 e n u +b(e n−1 u , u (t n )) +b(u n−1 , e n u ) + Se n−1 where We are now in a position to state and prove two error estimates for the Euler semiimplicit scheme (17).
Proof. Taking the inner product of (20) with 2d t e n u τ, (21) with 2d t e n B τ, we deduce that Using (7)- (12) and Lemma 1, we obtain Similarly, we have Combining the above inequalities with (28), and using Theorem 2, we can derive 1 ). Multiplying this inequality by σ(t n ) and taking the sum with respect to n from 1 to m, thanks to Theorems 2 and 5, we obtain Then, by applying Lemma 2 to this inequality and using Theorem 5, we have for all 1 ≤ m ≤ N. Finally, using (7)-(12) and the LBB condition (cf. [4,13]) we derive Multiplying this inequality by σ(t n ) and taking the sum with respect to n from 1 to m, due to Theorems 2 and 5 and (29), we obtain for all 1 ≤ m ≤ N. Using (29), (31), (9), and Lemma 1, we obtain Theorem 6.

PFEM for the MHD Equations
We further consider a spatial discretization for the penalty MHD system of time discretization in this section (cf. [4]). J h is a family of quasi-uniformly regular partitions of D into triangles or tetrahedron elements K with the diameter h K . Let the mesh size Assumption 4. The finite element space (X h , M h ) satisfies the discrete LBB condition (cf. [4,22]) where β 1 is a positive constant depending on D.
together with the inverse inequalities Next, to obtain an approximation of (u, p, B), we consider the following finite element pairs: {b} is a bubble function. Let b ∈ H 1 0 (K) take the value 1 at the barycentre of K and satisfy 0 ≤ b(x) ≤ 1, which is called a "bubble function" (cf. [22]). Furthermore, we denote the discrete subspace X 0h of X 0 as The finite element approximation for (16) based on X h × M h × W h is given as follows: where r 1h : L 2 (D) → X 0h and r 2h : L 2 (D) → W h are L 2 -orthogonal projectors. According to (34)-(36), these operators satisfy the following properties (cf. [3,4,22]): Here, we define the discrete Stokes operator A 1h = −r 1h ∆ h , which is defined by (see [4]) To obtain the error analysis of the scheme in the following section, we need the following discrete estimates which are obtained from [4,7]).
For the finite element space X h × M h × W h given above, problems (37) and (38) allow us to calculate velocity and pressure separately, i.e., (37) and (38) can be reduced as follows: In the following, we give the algebraic matrix form for the 2D case, and the algebraic matrix form for the 3D case is similar. In order to analyze the detailed form of the coefficient matrix for this scheme, we write the fluid velocity and magnetic field vectors and the corresponding test function vectors v n h = (v n 1h , v n 2h ), C n h = (C n 1h , C n 2h ), then, we expand (42) and obtain that Step Step 2. Find p n h from Here, ∂ x v 1h = ∂v 1h /∂x and ∂ y v 2h = ∂v 2h /∂y. Next, we assume the spaces X h and W h are combined with the basis functions where N and M denote the number of the basis functions in each of spaces. Then, Next, we show the relationship between the (u n−1 h , B n−1 h ) and (u n h , B n h ). Apparently, step 1 of the fully discrete PFEM generates an algebraic system as follows: Specifically, detailed calculation for PFEM gives Then, we obtain p n by step 2. Arguing in exactly the same way as in the proof of Theorem 4, using (9) and Lemma 1, we obtain a priori bound of schemes (37)-(39).

Error Analysis for the Fully Discrete Euler Semi-Implicit Scheme
In this section, we establish the error estimates for (u n h , p n h , B n h ) of the fully discrete Euler semi-implicit scheme (37)-(39). To this end, subtracting (37)-(38) from (16), we have In order to derive estimates of the error, we need the following regularity results.

Lemma 4.
Under Assumptions 1-3, we have the following estimates: We refer to [4,32] for the proof of these results.

Theorem 9.
Under Assumptions 1-4, we have Proof. To obtain the error estimates of the PFEM, we give the Galerkin projector R 1h : for all (u, p) ∈ (X, M) with ∇ · u + ν p = 0. In addition, R 2h : W → W h is the H 1 -orthogonal projector defined by (cf. [4]) for all B ∈ W. Due to the properties of R 1h (u, p), Q h (u, p) and R 2h , we have and Letting ξ n u = R 1h (u n , p n ) − u n h , ξ n p = Q h (u n , p n ) − p n h , and ξ n B = R 2h B n − B n h , we derive from (43) and (44) that Taking v h = 2d t ξ n u τ and q h = 2ξ n p τ in (50) and C h = 2d t ξ n B τ in (51), we obtain By using (7)- (12), (48) and (49), we obtain Combining the above inequalities with (52), and using Theorems 4 and 8, we can derive . Multiplying this inequality by σ(t n ) and taking the sum with respect to n from 1 to m, thanks to Theorems 4 and 8 and Lemma 2, we obtain Moreover, by using (48) and (49) and Theorem 4, we obtain Hence, by combining (53) with (54), we obtain (46). Next, using (7)-(12), (33) and (43), we derive Multiplying this inequality by σ(t n ) and taking the sum with respect to n from 1 to m, due to Theorems 8 and 46, we have Using (55) and (48) and (43), we obtain (47). Thus, this proof is thus complete.
Taking v h = 2ξ n u τ and q h = 2ξ n p τ in (43) and C h = 2ξ n B τ in (44), and adding these equations together, we have By using (11), (12), and Lemma 3, we obtain Combining these inequalities with (61) and using (40) and (41), we have Multiplying this inequality by σ(t n ) and taking the sum with respect to n from 1 to m, thanks to Theorems 4 and 8 and Lemma 2, we obtain Then, by applying (40), (41), Theorem 4, and (62), we obtain Theorem 10.

Error Estimates
Combining Theorem 3 and the results in Sections 3-5, we obtain the following results on convergence of the fully discrete Euler semi-implicit scheme.
Theorem 11. Under Assumptions 1-4, we have the following error estimates

Numerical Example
In this part, we present two numerical tests to validate the accuracy and performance of our scheme. We use the P b 1 /P 1 element that satisfies the LBB condition for velocity and pressure (u, p), and the P 1 element for magnetic field B. The penalty parameter is selected as = τ in all the numerical tests.
We choose τ = h 2 and h = 1/n (n = 8, 16, 32, 64, 128 in R 2 or n = 4, 8, 12, 16, 20 R 3 ). The numerical errors and the space convergence rates of the PFEM based on the semi-implicit scheme at t n = 1s are presented in Tables 1 and 2. We observe the first order accuracy for H 1 errors of u, B, and the second order accuracy asymptotically for L 2 errors of u, B, which are consistent with our theoretical results. These convergence rates are consistent with the expected orders. Notice that L 2 errors of p has a faster convergence rate than the theoretical results.

Two-Sided Lid-Driven Square Cavity Flow
In this example, we test the 2D/3D two-sided driven cavity flow problem (cf. [33]). The problem we study is the incompressible viscous flow in a square cavity whose top and bottom walls move in the same (parallel) or opposite (antiparallel) direction. We take the initial values u 0 = B 0 = 0 and the source terms f = g = 0. In the 2D case, we set a computational domain as D = [0, 1] 2 . The two boundary conditions are shown below: We set h = 1 60 , τ = 1 3600 . First, we consider the upper and lower walls moving in the same direction at the same speed along the x-axis. Figure 1 shows the velocity streamlines for the fluid Reynolds number R e = 1, 100, 1000, the magnetic Reynolds number R m = 1, and the coupling coefficient S = 10. It can be observed that the velocity streamlines are symmetric lines parallel to these walls and pass through the center of the cavity. With the increase of the fluid Reynolds number R e , the centers of the two symmetric vortices move to the right, and the two symmetric vortices become four symmetric vortices. Figure 2 shows the velocity streamlines for the fluid Reynolds number R e = 100, the magnetic Reynolds number R m = 1 and the coupling coefficient S = 1, 100, 1000. With the increase of coupling coefficient S, the two symmetric large vortices can split into more and more small vortices.  Then, we consider the upper and lower walls moving in the opposite direction at the same speed along the x-axis. Figure 3 gives the velocity streamlines for the fluid Reynolds number R e = 1, 100, 1000, the magnetic Reynolds number R m = 1 and the coupling coefficient S = 1. We find that, with the increase of the fluid Reynolds number R e , the centers of the two symmetric vortices shift to the upper right corner and the lower left corner, respectively. Figure 4 presents the velocity streamlines for the fluid Reynolds number R e = 1, the magnetic Reynolds number R m = 1 and the coupling coefficient S = 100, 1000, 10000. It can be observed from the figure that, with the increase of coupling coefficient S, the two symmetric large vortices become four small vortices.   We set h = 1 12 , τ = 1 600 . First, we consider the top and bottom walls moving in the same direction at the same speed. Figure 5 shows the velocity streamlines at plane y = 0.5 for the fluid Reynolds number R e = 1, 100, 500, the magnetic Reynolds number R m = 1 and the coupling coefficient S = 1. We find that, with the increase of the fluid Reynolds number R e , the centers of the two symmetric vortices move to the right. Figure 6 shows the velocity streamlines at plane y = 0.5 for the fluid Reynolds number R e = 10, the magnetic Reynolds number R m = 1 and the the coupling coefficient S = 1, 100, 500. We find that the two symmetric large vortices can split into more and more small vortices.  Then, we consider the top and bottom walls moving in the opposite direction at the same speed. Figure 7 gives the velocity streamlines at plane y = 0.5 for the fluid Reynolds number R e = 1, 100, 500, the magnetic Reynolds number R m = 1, and the coupling coefficient S = 1. We can see the centers of the two symmetric vortices shift to the upper right corner and the lower left corner, respectively. Figure 8 presents the velocity streamlines at plane y = 0.5 for the fluid Reynolds number R e = 1, the magnetic Reynolds number R m = 1, and the coupling coefficient S = 1, 100, 500. With the increase of coupling coefficient S, the centers of the two symmetric vortices are slightly offset.

Conclusions
We present the fully discrete PFEM for the 2D/3D unsteady MHD equations in this paper. We introduce a penalty term to decouple the MHD equations into two small equations: one is the equations of velocity and magnetic field (u, B), and the other is the equation of pressure p. Furthermore, we derive the error estimates for our scheme. Finally, two 2D/3D numerical experiments are given to verify the theoretical results.