Convergence Analysis of the Fully Decoupled Linear Scheme for Magnetohydrodynamic Equations

Abstract

In this paper, we propose a fully decoupled and linear scheme for the magnetohydrodynamic (MHD) equation with the backward differential formulation (BDF) and finite element method (FEM). To solve the system, we adopt a technique based on the “zero-energy-contribution” contribution, which separates the magnetic and fluid fields from the coupled system. Additionally, making use of the pressure projection methods, the pressure variable appears explicitly in the velocity field equation, and would be computed in the form of a Poisson equation. Therefore, the total system is divided into several smaller sub-systems that could be simulated at a significantly low cost. We prove the unconditional energy stability, unique solvability and optimal error estimates for the proposed scheme, and present numerical results to verify the accuracy, efficiency and stability of the scheme.

Share and Cite:

Xia, Z. and Xu, Q. (2022) Convergence Analysis of the Fully Decoupled Linear Scheme for Magnetohydrodynamic Equations. Journal of Applied Mathematics and Physics, 10, 3462-3474. doi: 10.4236/jamp.2022.1011228.

1. Introduction

The MHD system, consisting of the Navier-Stokes equation and Maxwell’s equation [1] [2], is widely applied into plasma and liquid metal processing, and its numerical simulation is also of great importance [3] [4]. The incompressible MHD system reads as

μ t H + σ 1 × ( × H ) μ × ( u × H ) = 0, (1.1)

t u + u u ν Δ u + p + μ H × ( × H ) = 0, (1.2)

u = 0, (1.3)

over Ω × ( 0, T ] , where Ω is a bounded and convex polyhedral domain in 3 (polygonal domain in 2 ). In the above system, unknowns H , u and p denote the magnetic filed, the velocity field, and the pressure, respectively. Coefficient σ denotes the magnetic Reynolds number, ν is the viscosity of the fluid and μ = M 2 ν σ 1 , where M is the Hartman number. All the constants here are positive. The initial data and boundary conditions are as follows

H | t = 0 = H 0 , u | t = 0 = u 0 , in Ω , (1.4)

H × n = 0, u = 0, on Ω × ( 0, T ] . (1.5)

We assume that the initial data satisfies

H 0 = u 0 = 0. (1.6)

Taking divergence of (1.1) leads to H = 0 for all t > 0 .

There are already many works dedicated to theoretical analysis of the MHD equations [5] [6] [7]. For the numerical schemes of MHD system, since the weak solution of the system lies in H 1 ( Ω ) space, much research work is devoted to the application of H 1 ( Ω ) conforming elements. Li et al. used H 1 ( Ω ) conforming elements in the constrained transport type framework and reconstructed a completely non-divergence magnetic field. Gunzburger et al. studied the approximation method of finite element and described three iterative methods for solving finite element model. By using H 1 ( Ω ) conforming element, the optimal error estimates for stationary MHD systems are analyzed [8]. He et al. extended this result to the time-dependent MHD model [9]. More research work on H 1 ( Ω ) conforming elements can be seen in [10] [11] [12] [13].

There are increasing interests devoted to the design of temporally discrete schemes for the MHD system recently. By combining the projection method of Navier-Stokes equations with some subtle implicit and explicit techniques on nonlinear coupling terms, Zhang et al. developed a linear and fully decouple first-order time advance scheme for solving MHD equations, and established unconditional energy stability and strict error estimates for the scheme [14] theoretically. Yang et al. proposed a second-order fully decoupled unconditional energy stabilization finite element algorithm for the incompressible MHD equation [15]. Li et al. proposed a stable mixed finite element method for solving three-dimensional incompressible MHD equations, and gave an iterative format for discrete problems with its adaptability. This method can guarantee long time stable simulation and allow large time step simulation [16]. Wang et al. designed a second-order accurate, fully discrete numerical scheme based on the decoupled Crank-Nicolson method and H 1 ( Ω ) -conforming finite element for the MHD system, and proved the unconditional energy stability, unique solvability and optimal error estimates [17].

A lot of work has been done to improve computational efficiency by dealing with nonlinear and coupled terms. In general, as long as the system of the PDE system is energy-stable, one would like to maintain stability in the design of numerical schemes. Therefore, a common approach is to linearize the nonlinear terms using implicit-explicit technique, that is, using a combination of interpolation and extrapolation. It can lead to a linear system with variable coefficients, which thus avoids the introduction of additional methods such as Newton’s iterative methods in solving the nonlinear discrete system. As a result, the computational efficiency is indeed improved without affecting the unconditional energy stability and accuracy. Furthermore, in the MHD equations, ( u ) u , H × ( × H ) and × ( u × H ) would vanish by taking inner products with some appropriate terms in the energy stability analysis, which is so called the “zero-energy-contribution” property. Based on this, Yang [18] proposed a fully decoupled, second-order temporally accuracy scheme for a binary fluidsurfactant phase-field model, and proved the energy stability of the numerical scheme. Chen et al. proposed two completely discrete time propulsion schemes, which are linear, decoupled, unconditionally energy stable and second-order. The energy stability of discrete problems was proved [19]. The above method can generate more smaller subsystems than existing methods, so the application of it also allows us to employ finer grids, which means it is possible to obtain more accurate results. In particular, the incompressible MHD equation has been solved by this method in [20]. More work on the “zero-energy contribution” property can be seen in [21] [22] [23] [24] [25] and so on. It is noted that these works only demonstrate the accuracy through numerical examples without rigorous theoretical analysis, while we manage to fill in this gap this time.

In our work, we design a fully decoupled scheme with the FEM and a second-order BDF, and give a rigorous proof of unconditional energy stability, unique solvability and optimal rates of convergence. We adopt H 1 ( Ω ) conforming elements to solve MHD system (1.1)-(1.3). The pressure projection method combined with the second-order BDF is utilized in temporal discretizations, and the nonlinear terms are treated as the explicit way based on the “zero-energy-contribution” terms, which ensures the discrete system to be linear and fully decoupled. We prove the unconditional energy stability, and then the corresponding homogeneous equations only admit zero solutions, which leads to the unique solvability immediately. In practical computation, instead of solving the whole system together, we implement the discrete scheme step by step. Therefore, we could utilize the conjugate gradient method to calculate the magnetic and velocity fields, which is significantly efficient, and the pressure field is obtained by solving a Poisson-type equation. Compared with the original system, the sizes of each matrix in actual calculation are much smaller. All these features make this method possess higher efficiency and finer mesh than the traditional pressure projection method. In addition, we provide the rigorous proof analysis for the numerical scheme with errors O ( h r + 1 + τ 2 ) in l ( 0, T ; L 2 ( Ω ) ) -norm, where h and τ are the spatial and temporal sizes respectively, and r is the degree of the polynomial functions.

The organizational structure of this paper is as follows. In Section 2, we introduce variational formulation of the MHD system, and propose a fully decoupled discrete scheme based on the second-order BDF-FEM. In Section 3, we provide the theoretic results of the energy stability, unique solvability and optimal error estimates for the numerical scheme. Finally, in Section 4 some numerical examples are illustrated to verify the accuracy, efficiency and stability.

2. Variational Formulation and Numerical Scheme

2.1. Variational Formulation

We adopt the conventional Sobolev space W k , p ( Ω ) of functions defined on Ω for k 0 and 1 p , and denote L p ( Ω ) = W 0, p ( Ω ) and H k ( Ω ) = W k ,2 ( Ω ) . Then we take the notations W 0 1, p ( Ω ) representing the space of functions in W 1, p ( Ω ) with zero traces on the boundary Ω , and naturally H 0 1 ( Ω ) : = W 0 1,2 ( Ω ) . The corresponding vector spaces are given by

L p ( Ω ) = [ L p ( Ω ) ] d , W k , p ( Ω ) = [ W k , p ( Ω ) ] d , W 0 1, p ( Ω ) = [ W 0 1, p ( Ω ) ] d , H 0 1 ( Ω ) = W 0 1,2 ( Ω ) , H 1 ( Ω ) = { v H 1 ( Ω ) : v × n | Ω = 0 } ,

where d is the dimension of space. As usual ( , ) denotes the inner product of L 2 ( Ω ) .

Consequently, we obtain that the exact solutions of (1.1)-(1.3) satisfy

μ ( t H , w ) + σ 1 ( × H , × w ) μ ( u × H , × w ) = 0, (2.1)

( t u , v ) + b ( u , u , v ) + ν ( u , v ) ( p , v ) + μ ( H × ( × H ) , v ) = 0, (2.2)

( u , q ) = 0, (2.3)

for any test functions ( w , v , q ) ( H 1 ( Ω ) , H 0 1 ( Ω ) , L 0 2 ( Ω ) ) . The trilinear operator b ( , , ) is defined as

b ( u , v , w ) : = ( u v , w ) + 1 2 ( ( u ) v , w ) = 1 2 [ ( u v , w ) ( u w , v ) ] , u , v , w H 0 1 ( Ω ) , (2.4)

and then it is natural that

b ( u , v , v ) = 0, u , v H 0 1 ( Ω ) . (2.5)

In order to treat the nonlinear term explicitly and ensure the energy stability, we introduce the function of “zero energy contribution”, which is defined as

Q t = μ ( u × H , × H ) + μ ( H × ( × H ) , u ) + b ( u , u , u ) 0, and Q 0 = 1. (2.6)

and then rewrite the Equations (2.1)-(2.3) as follows

μ ( t H , w ) + σ 1 ( × H , × w ) μ Q ( u × H , × w ) = 0, (2.7)

( t u , v ) + Q b ( u , u , v ) + ν ( u , v ) ( p , v ) + μ Q ( H × ( × H ) , v ) = 0, (2.8)

( u , q ) = 0, (2.9)

for any test functions ( w , v , q ) ( H 1 ( Ω ) , H 0 1 ( Ω ) , L 0 2 ( Ω ) ) . The variational Formulation (2.7)-(2.9) will be studied in this paper later.

2.2. Numerical Discretization

In this subsection, we propose a fully discrete decoupled finite element method for solving the system (1.1)-(1.3). Let h denote a quasi-uniform partition of Ω into tetrahedrons K j in 3 (or triangles in 2 ), j = 1,2, , M , with mesh size h = max 1 j M { diam K j } . To approximate u and p in the system (1.1)-(1.3), we introduce the Taylor-Hood finite element space X h × M h , defined by

X h = { l h H 0 1 ( Ω ) : l h | K j P r ( K j ) } , M h = { q h L 2 ( Ω ) : q h | K j P r 1 ( K j ) , Ω q h d x = 0 } ,

for any integer r 2 , where P r ( K j ) is the space of polynomials with degree r on K j for all K j h and P r ( K j ) : = [ P r ( K j ) ] d . To approximate the magnetic field H , we introduce the finite element space S h defined by

S h = { w h H 1 ( Ω ) : w h | K j P r ( K j ) } .

Let { t n = n τ } n = 0 N be a uniform partition of the time interval [ 0, T ] , and τ = T / N denotes the temporal step size. Furthermore, v n represents the value of v ( x , t n ) , and we adopt the abbreviation

v ˜ n + 1 : = 2 v n v n 1 .

Additionally, we adopt the following notations:

D t v h n + 1 = 3 2 v h n + 1 2 v h n + 1 2 v h n 1 τ ,

where v could be an arbitrary function except that u ^ h n + 1 , R h u n + 1 ^ and e ^ u n + 1 , for which we denote by

D t u ^ h n + 1 = 3 u ^ h n + 1 4 u h n + u h n 1 2 τ , D t R h u n + 1 ^ = 3 R h u n + 1 ^ 4 R h u n + R h u n 1 2 τ , D t e ^ u n + 1 = 3 e ^ u n + 1 4 e u n + e u n 1 2 τ .

Thus, we propose a fully decoupled discrete scheme with FEM and BDF for the incompressible MHD Equations (1.1)-(1.3) that, find Q h n + 1 and ( H h n + 1 , u h n + 1 , u ^ h n + 1 , p h n + 1 ) ( S h , X h , X h , M h ) such that

μ ( D t H h n + 1 , w h ) + σ 1 ( × H h n + 1 , × w h ) + σ 1 ( H h n + 1 , w h ) μ Q h n + 1 ( u ˜ h n + 1 × H ˜ h n + 1 , × w h ) = 0, (2.10)

( D t u ^ h n + 1 , v h ) + Q h n + 1 b ( u ˜ h n + 1 , u ˜ h n + 1 , v h ) + ν ( u ^ h n + 1 , v h ) ( p h n , v h ) + μ Q h n + 1 ( H ˜ h n + 1 × ( × H ˜ h n + 1 ) , v h ) = 0, (2.11)

( u h n + 1 u ^ h n + 1 τ , s h ) 2 3 ( p h n + 1 p h n , s h ) = 0, (2.12)

( u h n + 1 , q h ) = 0, (2.13)

D t Q h n + 1 = μ ( u ˜ h n + 1 × H ˜ h n + 1 , × H h n + 1 ) + μ ( H ˜ h n + 1 × ( × H ˜ h n + 1 ) , u ^ h n + 1 ) + b ( u ˜ h n + 1 , u ˜ h n + 1 , u ^ h n + 1 ) , (2.14)

for any ( w h , v h , s h , q h ) ( S h , X h , X h , M h ) and n = 1,2, , N 1 .

Remark 2.1. We have added a stabilization term σ 1 ( H h n + 1 , w h ) to (2.10), which is consistent with (2.7) in view of H = 0 . It facilitates to prove the optimal error estimates for magnetic fields.

Remark 2.2. The pressure field appears explicitly in the velocity Equation (2.11), and then it could be computed by solving the simple linear Equation (2.12). This is so called the pressure projection method, i.e., thedecoupledtechnique, and it divides the whole numerical system into two smaller ones. In details, the incompressible restriction would be removed in the computation and therefore more numerical techniques could be potentially be utilized.

Remark 2.3. The scheme is a multi-step method, so the starting values at time steps t 0 and t 1 should be given. In this work we make use of the backward Euler method which will not affect the accuracy and energy stability.

2.3. The Detailed Implementation of the Scheme

Through Q h n + 1 , we denote H h n + 1 and u ^ h n + 1 by

H h n + 1 = H 1 h n + 1 + Q h n + 1 H 2 h n + 1 and u ^ h n + 1 = u ^ 1 h n + 1 + Q h n + 1 u ^ 2 h n + 1 , (2.15)

with which, (2.10) and (2.11) could be rewritten as

3 2 μ ( H 1 h n + 1 + Q h n + 1 H 2 h n + 1 τ , w h ) + σ 1 ( × ( H 1 h n + 1 + Q h n + 1 H 2 h n + 1 ) , × w h ) + σ 1 ( ( H 1 h n + 1 + Q h n + 1 H 2 h n + 1 ) , w h ) μ Q h n + 1 ( u ˜ h n + 1 × H ˜ h n + 1 , × w h ) = μ ( 4 H h n H h n 1 2 τ , w h ) , (2.16)

3 2 ( u ^ 1 h n + 1 + Q h n + 1 u ^ 2 h n + 1 τ , v h ) + Q h n + 1 b ( u ˜ h n + 1 , u ˜ h n + 1 , v h ) + ν ( ( u ^ 1 h n + 1 + Q h n + 1 u ^ 2 h n + 1 ) , v h ) ( P h n , v h ) + μ Q h n + 1 ( H ˜ h n + 1 × ( × H ˜ h n + 1 ) , v h ) = ( 4 u h n u h n 1 2 τ , v h ) . (2.17)

Step 1: We write the above formula as

3 2 μ ( H 1 h n + 1 τ , w h ) + σ 1 ( × H 1 h n + 1 , × w h ) + σ 1 ( H 1 h n + 1 , w h ) = μ ( 4 H h n + 1 H h n 1 2 τ , w h ) , (2.18)

3 2 μ ( H 2 h n + 1 τ , w h ) + σ 1 ( × H 2 h n + 1 , × w h ) + σ 1 ( H 2 h n + 1 , w h ) = μ ( u ˜ h n + 1 × H ˜ h n + 1 , × w h ) , (2.19)

and

3 2 ( u ^ 1 h n + 1 τ , v h ) + ν ( u ^ 1 h n + 1 , v h ) = ( p h n , v h ) + ( 4 u h n u h n 1 2 τ , v h ) , (2.20)

3 2 ( u ^ 2 h n + 1 τ , v h ) + ν ( u ^ 2 h n + 1 , v h ) = b ( u ˜ h n + 1 , u ˜ h n + 1 , v h ) μ ( H ˜ h n + 1 × ( × H ˜ h n + 1 ) , v h ) . (2.21)

Thus in step 1, by solving (2.18)-(2.21) we get the values of H 1 h n + 1 , H 2 h n + 1 , u ^ 1 h n + 1 and u ^ 2 h n + 1 .

Step 2: Substituting (2.15) into (2.14) leads to

3 Q h n + 1 4 Q h n + Q h n 1 2 τ = μ ( H ˜ h n + 1 × ( × H ˜ h n + 1 ) , u ^ 1 h n + 1 + Q h n + 1 u ^ 2 h n + 1 ) + b ( u ˜ h n + 1 , u ˜ h n + 1 , u ^ 1 h n + 1 + Q h n + 1 u ^ 2 h n + 1 ) μ ( u ˜ h n + 1 × H ˜ h n + 1 , × ( H 1 h n + 1 + Q h n + 1 H 2 h n + 1 ) ) : = T 1 + Q h n + 1 T 2 ,

which yields

Q h n + 1 = 2 Q h n 1 2 Q h n 1 + τ T 1 3 2 τ T 2 , (2.22)

where we denote

T 1 = μ ( H ˜ h n + 1 × ( × H ˜ h n + 1 ) , u ^ 1 h n + 1 ) + b ( u ˜ h n + 1 , u ˜ h n + 1 , u ^ 1 h n + 1 ) μ ( u ˜ h n + 1 × H ˜ h n + 1 , × H 1 h n + 1 ) , T 2 = μ ( H ˜ h n + 1 × ( × H ˜ h n + 1 ) , u ^ 2 h n + 1 ) + b ( u ˜ h n + 1 , u ˜ h n + 1 , u ^ 2 h n + 1 ) μ ( u ˜ h n + 1 × H ˜ h n + 1 , × H 2 h n + 1 ) .

By adopting w h = H 2 h n + 1 in (2.19) and v h = u ^ 2 h n + 1 in (2.21), we have

T 2 = 3 2 τ ( μ H 2 h n + 1 L 2 2 + u ^ 2 h n + 1 L 2 2 ) + σ 1 ( × H 2 h n + 1 L 2 2 + H 2 h n + 1 L 2 2 ) + ν u ^ 2 h n + 1 L 2 2 0,

which means 3 2 τ T 2 > 0 , so that (2.22) is always solvable for Q h n + 1 .

Step 3: For now, through (2.15) we can obtain H h n + 1 and u ^ h n + 1 . Then u h n + 1 and p h n + 1 follow immediately by solving (2.12) and (2.13).

Remark 2.4. In the practical implementation, in Step 1 we obtain an elliptic equation with constant coefficients, which corresponds to a symmetric positive stiffness matrix in the linear algebraic system, so that we could employ the conjugate gradient method to solve it extremely efficiently. It is also the reason why we adopt thedecoupledtechnique, as stated in Remark 2.2. Step 2 is just a direct algebraic calculation, and Step 3 is to solve a Poission-type equation, which is relatively fast to compute. Hence the computation of this scheme is lower than the traditional decoupled scheme, and meanwhile the adoption of smaller sub-systems allows us to employ the finer grids. We will illustrate the efficiency and accuracy later.

3. Theoretical Results

For compactness, the proof in this section is referred to [26].

3.1. Discrete Energy Stability

In this subsection, the discrete energy stability of the scheme (2.10)-(2.14) will be proven. We define the discrete gradient operator h : M h X h as

( v h , h q h ) = ( v h , q h ) , v h X h , q h M h (3.1)

and give the following energy stability theorem.

Theorem 3.1. The numerical solution ( H h n , u h n , p h n ) to the fully discrete scheme (2.10)-(2.14) is uniquely solvable and satisfies the following energy estimate

ε h n + 1 ε h n 0 , n = 2 , 3 , , N 1 ,

where the discrete energy function ε h n is defined as

ε h n = 1 4 ( μ H h n L 2 2 + μ 2 H h n H h n 1 L 2 2 + u h n L 2 2 + 2 u h n u h n 1 L 2 2 + ( Q h n + 1 ) 2 + ( 2 Q h n + 1 Q h n ) 2 ) + τ 2 3 h p h n L 2 2 .

3.2. Optimal Error Estimates

We make the following regularity assumption that the continuous system (1.1)-(1.3) admits unique solutions satisfying

H t t t L ( 0, T ; L 2 ) + H t t L ( 0, T ; H 1 ) + H t L ( 0, T ; H r + 1 ) + H L ( 0, T ; H r + 2 ) + u t t t L ( 0, T ; L 2 ) + u t t L ( 0, T ; H 1 ) + u t L ( 0, T ; H r + 1 ) + u L ( 0, T ; H r + 1 ) + p t t L ( 0, T ; L 2 ) + p t L ( 0, T ; H r + 1 ) M , (3.2)

where v t = t v , and the optimal error estimates are stated in the following theorem.

Theorem 3.2 Suppose that the exact solution ( H , u , p ) to the original Equations (1.1)-(1.3) satisfies the regularity assumption (3.2), and the numerical scheme (2.10)-(2.13) admits the unique solution ( H h n , u h n , p h n ) ,2 n N . Then for some positive constants τ 0 and h 0 , when τ < τ 0 , h < h 0 , we have

max 2 n N ( H h n H n L 2 + u h n u n L 2 ) C 0 ( τ 2 + h r + 1 ) , (3.3)

where C 0 is a positive constant independent of τ and h.

4. Numerical Examples

In this section, we carried out numerical experiments using the software FreeFEM++.

4.1. Accuracy Test

For sake of brevity, we consider the incompressible MHD equations

μ t H + σ 1 × ( × H ) μ × ( u × H ) = J , (4.1)

t u + u u ν Δ u + p + μ H × ( × H ) = f , (4.2)

H = 0, u = 0, (4.3)

in a two-dimensional domain [ 0,2 π ] × [ 0,2 π ] , with the initial and boundary conditions (1.4)-(1.5), where the source terms J and f are chosen correspondingly to the exact solutions

u = t 8 ( sin 2 x sin ( 2 y ) sin ( 2 x ) sin 2 y ) , H = t 6 ( sin y cos x sin x cos y ) , p = t 5 sin 2 x sin 2 y . (4.4)

All the coefficients in (4.1)-(4.3) are chosen to be 1, and we take the final time T = 1 . Note that the above exact solutions u and H satisfy the divergence-free conditions.

We first solve the MHD system (4.1)-(4.3) by the scheme (2.10)-(2.13) with a quadratic finite element approximation for H and u , and a linear finite element approximation for p. To emphasize the convergence rate in time, a sufficiently small spatial mesh size h = 1 / 120 is chosen such that the spatial discretization error can be relatively negligible. The time step is τ = T / N with N = 40 , 80 , 160 . We show the numerical results at time T = 1 in Table 1, which indicate that the proposed scheme is convergent with second-order temporally accuracy.

Then we solve the problem (4.1)-(4.3) by the scheme (2.10)-(2.13) with a sufficiently small temporal step τ = 1 / 5000 , to observe the spatial convergence rates. Take spatial size as h = 1 / 10,1 / 20,1 / 40,1 / 80 . Again, a quadratic finite element approximation for H and u is adopted, combined with a linear finite element approximation for p. Numerical results at T = 1 are presented in Table 2.

It is clear to see the third-order accuracy, which is consistent with the theoretical analysis in Theorem 3.2.

4.2. Efficiency Test

As a comparison, we present another numerical scheme for the MHD based on the traditional “decoupled” technique, i.e., finding ( H h n + 1 , u h n + 1 , u ^ h n + 1 , p h n + 1 ) ( S h , X h , X h , M h ) such that it holds that

Table 1. Temporal accuracy.

Table 2. Spatial accuracy.

μ ( D t H h n + 1 , w h ) + σ 1 ( × H h n + 1 , × w h ) + σ 1 ( H h n + 1 , w h ) μ ( u ˜ h n + 1 × H ˜ h n + 1 , × w h ) = 0, ( D t u ^ h n + 1 , v h ) + b ( u ˜ h n + 1 , u ˜ h n + 1 , v h ) + ν ( u ^ h n + 1 , v h ) ( p h n , v h ) + μ ( H ˜ h n + 1 × ( × H ˜ h n + 1 ) , v h ) = 0, ( u h n + 1 u ^ h n + 1 τ , s h ) 2 3 ( p h n + 1 p h n , s h ) = 0, ( u h n + 1 , q h ) = 0,

for any ( w h , v h , s h , q h ) ( S h , X h , X h , M h ) and n = 1 , 2 , , N 1 . Using the same parameters and initial boundary conditions as those in Subsection 4.1, we can obtain the following results.

As seen in the above Table 3 and Table 4, the errors generated by the scheme provided in this subsection are almost same as those in Subsection 4.1, but it obviously takes longer time than before.

4.3. Energy Stability

Finally, we carry out the numerical experiment to verify the discrete energy stability. Choose the initial data as

u 1 = u 0 = ( sin 2 x sin ( 2 y ) sin ( 2 x ) sin 2 y ) , H 1 = H 0 = ( sin y cos x sin x cos y ) , p 0 = sin ( 2 x ) sin ( 2 y ) , (4.5)

We fix τ = 0.1 and h = 1 / 200 , and define the discrete energy function defined in Theorem 3.1. Meanwhile, we still adopt the quadratic elements for ( H , u ) and linear elements for p. Up to the final time T = 10 , we show the energy decaying curve as follows (Figure 1).

Table 3. Temporal accuracy.

Table 4. Spatial accuracy.

Figure 1. Discrete energy evolution of the incompressible MHD system.

5. Conclusion

In this work, we have designed a fully decoupled second-order BDF scheme together with FEM for the incompressible MHD system (1.1)-(1.3), and then have carried out the rigorous proof of unconditional energy stability, unique solvability and optimal error estimates. The fully decoupled method adopted in this work is an efficient approach to deal with the incompressible constraint and nonlinear terms, and thus the technique could be applied to the other incompressible flow system, for example, the multi-phase MHD system. Although the method based on “zero-energy-contribution” property applied to the incompressible MHD equation has been presented in [20] already, we have managed to obtain the optimal convergence rates, which until now could not be found in the existing works.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1] Gerbeau, J., Le Bris, C. and Lelievre, T. (2006) Mathematical Methods for the Magnetohydrodynamics of Liquid Metals. Numerical Mathematics and Scientic Computation. Oxford University Press, Oxford. https://doi.org/10.1093/acprof:oso/9780198566656.001.0001
[2] London, S.D. (2018) Solitary Waves in Shallow Water Hydrodynamics and Magnetohydrodynamics in Rotating Spherical Coordinates. Geophys. Astrophys. Fluid Dyn., 112, 62-73. https://doi.org/10.1080/03091929.2017.1383407
[3] Asai, S. (2012) Magnetohydrodynamics in Materials Processing. Electromagnetic Processing of Materials. https://doi.org/10.1007/978-94-007-2645-1
[4] Unger, Y., Mond, M. and Branover, H. (2015) Liquid Metal Ows: Magnetohydrodynamics and Application. Progress in Astronautics and Aeronautics.
[5] Chen, F., Li, Y. and Zhao, Y. (2017) Global Well-Posedness for the Incompressible MHD Equations with Variable Viscosity and Conductivity. J. Math. Anal. Appl., 447, 1051-1071. https://doi.org/10.1016/j.jmaa.2016.10.059
[6] He, C. and Wang, Y. (2007) On the Regularity Criteria for Weak Solutions to the Magnetohydrodynamic Equations. J. Dierential Equations, 238, 1-17. https://doi.org/10.1016/j.jde.2007.03.023
[7] Ji, R. and Tian, L. (2021) Stability of the 3D Incompressible MHD Equations with Horizontal Dissipation in Periodic Domain. AIMS Math., 6, 11837-11849. https://doi.org/10.3934/math.2021687
[8] Gunzburger, M.D., Meir, A.J. and Peterson, J.S. (1991) On the Existence, Uniqueness, and Nite Element Approximation of Solutions of the Equations of Stationary, Incompressible Magnetohydrodynamics. Math. Comp., 56, 523-563. https://doi.org/10.1090/S0025-5718-1991-1066834-0
[9] He, Y. (2015) Unconditional Convergence of the Euler Semi-Implicit Scheme for the Three-Dimensional Incompressible MHD Equations. IMA J. Numer. Anal., 35, 767-801. https://doi.org/10.1093/imanum/dru015
[10] Adak, D. and Natarajan, S. (2021) On the H1 Conforming Virtual Element Method for Time Dependent Stokes Equation. Math. Comput. Sci., 15, 135-154. https://doi.org/10.1007/s11786-020-00473-1
[11] He, Y. and Zou, J. (2018) A Priori Estimates and Optimal Nite Element Approximation of the MHD Ow in Smooth Domains. ESAIM Math. Model. Numer. Anal., 52, 181-206. https://doi.org/10.1051/m2an/2018006
[12] Ravindran, S.S. (2018) Analysis of a Decoupled Time-Stepping Algorithm for Reduced MHD System Modeling Magneto-Convection. Numer. Methods Partial Dierential Equations, 34, 1953-1974. https://doi.org/10.1002/num.22270
[13] Li, L. and Zheng, W. (2017) A Robust Solver for the Nite Element Approximation of Stationary Incompressible MHD Equations in 3D. J. Comput. Phys., 351, 254-270. https://doi.org/10.1016/j.jcp.2017.09.025
[14] Zhang, G.-D., He, X. and Yang, X. (2020) Fully Decoupled, Linear and Unconditionally Energy Stable Time Discretization Scheme for Solving the Magneto-Hydrodynamic Equations. J. Comput. Appl. Math., 369, 112636. https://doi.org/10.1016/j.cam.2019.112636
[15] Yang, J. and Mao, S. (2021) Second Order Fully Decoupled and Unconditionally Energy-Stable Nite Element Algorithm for the Incompressible MHD Equations. Appl. Math. Lett., 121, Article No. 107467. https://doi.org/10.1016/j.aml.2021.107467
[16] Li, L., Zhang, D. and Zheng, W. (2021) A Constrained Transport Divergence-Free Nite Element Method for Incompressible MHD Equations. J. Comput. Phys., 428, Article No. 109980. https://doi.org/10.1016/j.jcp.2020.109980
[17] Wang, C., Wang, J., Xia, Z. and Xu, L. (2022) Optimal Error Estimates of a Crank-Nicolson Nite Element Projection Method for Magnetohydrodynamic Equations. ESAIM Math. Model. Numer. Anal., 56, 767-789. https://doi.org/10.1051/m2an/2022020
[18] Yang, X. (2021) On a Novel Fully Decoupled, Second-Order Accurate Energy Stable Numerical Scheme for a Binary Uid Surfactant Phase Eld Model. SIAM J. Sci. Comput., 43, B479-B507. https://doi.org/10.1137/20M1336734
[19] Chen, R. and Zhang, H. (2020) Second-Order Energy Stable Schemes for the New Model of the Cahn-Hilliard-MHD Equations. Adv. Comput. Math., 46, Article No. 79. https://doi.org/10.1007/s10444-020-09822-x
[20] Zhang, G.-D., He, X. and Yang, X. (2022) A Fully Decoupled Linearized Nite Element Method with Second-Order Temporal Accuracy and Unconditional Energy Stability for Incompressible MHD Equations. J. Comput. Phys., 448, Article No. 110752. https://doi.org/10.1016/j.jcp.2021.110752
[21] Li, T., Liu, P., Zhang, J. and Yang, X. (2022) Ecient Fully Decoupled and Second-Order Time-Accurate Scheme for the Navier-Stokes Coupled Cahn-Hilliard Ohta-Kawaski Phase Eld Model of Diblock Copolymer Melt. J. Comput. Appl. Math., 403, Article No. 113843. https://doi.org/10.1016/j.cam.2021.113843
[22] Yang, X. (2021) Ecient and Energy Stable Scheme for the Hydrodynamically Coupled Three Components Cahn-Hilliard Phase Eld Model Using the Stabilized-Invariant Energy Quadratization (S-IEQ) Approach. J. Comput. Phys., 438, Article No. 110342. https://doi.org/10.1016/j.jcp.2021.110342
[23] Yang, X. (2021) A Novel Decoupled Second-Order Time Marching Scheme for the Two-Phase Incompressible Navier-Stokes/Darcy Coupled Nonlocal Allen-Cahn Model. Comput. Methods Appl. Mech. Engrg., 377, Article No. 113597. https://doi.org/10.1016/j.cma.2020.113597
[24] Yang, X. (2021) A Novel Fully Decoupled Scheme with Second-Order Time Accuracy and Unconditional Energy Stability for the Navier-Stokes Equations Coupled with Mass-Conserved Allen-Cahn Phase Eld Model of Two-Phase Incompressible Ow. Internat. J. Numer. Methods Engrg., 122, 1283-1306. https://doi.org/10.1002/nme.6578
[25] Yang, X. (2021) A Novel Fully-Decoupled, Second-Order Time-Accurate, Unconditionally Energy Stable Scheme for a Ow-Coupled Volume-Conserved Phase Eld Elastic Bending Energy Model. J. Comput. Phys., 432, Article No. 110015. https://doi.org/10.1016/j.jcp.2020.110015
[26] Xu, Q. (2022) Finite Element Methods with High-Order Temporally Discretizations for the Maxwell’s Equation and Its Coupled Model. University of Electronic Science and Technology of China, Master’s Thesis.

Copyright © 2024 by authors and Scientific Research Publishing Inc.

Creative Commons License

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.