Convergence Analysis of the Fully Decoupled Linear Scheme for Magnetohydrodynamic Equations ()
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
(1.1)
(1.2)
(1.3)
over
, where
is a bounded and convex polyhedral domain in
(polygonal domain in
). In the above system, unknowns
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
, where M is the Hartman number. All the constants here are positive. The initial data and boundary conditions are as follows
(1.4)
(1.5)
We assume that the initial data satisfies
(1.6)
Taking divergence of (1.1) leads to
for all
.
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
space, much research work is devoted to the application of
conforming elements. Li et al. used
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
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
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
-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,
,
and
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
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
in
-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
of functions defined on
for
and
, and denote
and
. Then we take the notations
representing the space of functions in
with zero traces on the boundary
, and naturally
. The corresponding vector spaces are given by
where d is the dimension of space. As usual
denotes the inner product of
.
Consequently, we obtain that the exact solutions of (1.1)-(1.3) satisfy
(2.1)
(2.2)
(2.3)
for any test functions
. The trilinear operator
is defined as
(2.4)
and then it is natural that
(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
(2.6)
and then rewrite the Equations (2.1)-(2.3) as follows
(2.7)
(2.8)
(2.9)
for any test functions
. 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
denote a quasi-uniform partition of
into tetrahedrons
in
(or triangles in
),
, with mesh size
. To approximate
and p in the system (1.1)-(1.3), we introduce the Taylor-Hood finite element space
, defined by
for any integer
, where
is the space of polynomials with degree r on
for all
and
. To approximate the magnetic field
, we introduce the finite element space
defined by
Let
be a uniform partition of the time interval
, and
denotes the temporal step size. Furthermore,
represents the value of
, and we adopt the abbreviation
Additionally, we adopt the following notations:
where v could be an arbitrary function except that
,
and
, for which we denote by
Thus, we propose a fully decoupled discrete scheme with FEM and BDF for the incompressible MHD Equations (1.1)-(1.3) that, find
and
such that
(2.10)
(2.11)
(2.12)
(2.13)
(2.14)
for any
and
.
Remark 2.1. We have added a stabilization term
to (2.10), which is consistent with (2.7) in view of
. 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., the “decoupled” technique, 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
and
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
, we denote
and
by
(2.15)
with which, (2.10) and (2.11) could be rewritten as
(2.16)
(2.17)
Step 1: We write the above formula as
(2.18)
(2.19)
and
(2.20)
(2.21)
Thus in step 1, by solving (2.18)-(2.21) we get the values of
,
,
and
.
Step 2: Substituting (2.15) into (2.14) leads to
which yields
(2.22)
where we denote
By adopting
in (2.19) and
in (2.21), we have
which means
, so that (2.22) is always solvable for
.
Step 3: For now, through (2.15) we can obtain
and
. Then
and
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 the “decoupled” technique, 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
as
(3.1)
and give the following energy stability theorem.
Theorem 3.1. The numerical solution
to the fully discrete scheme (2.10)-(2.14) is uniquely solvable and satisfies the following energy estimate
where the discrete energy function
is defined as
3.2. Optimal Error Estimates
We make the following regularity assumption that the continuous system (1.1)-(1.3) admits unique solutions satisfying
(3.2)
where
, and the optimal error estimates are stated in the following theorem.
Theorem 3.2 Suppose that the exact solution
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
. Then for some positive constants
and
, when
,
, we have
(3.3)
where
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
(4.1)
(4.2)
(4.3)
in a two-dimensional domain
, with the initial and boundary conditions (1.4)-(1.5), where the source terms
and
are chosen correspondingly to the exact solutions
(4.4)
All the coefficients in (4.1)-(4.3) are chosen to be 1, and we take the final time
. Note that the above exact solutions
and
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
and
, and a linear finite element approximation for p. To emphasize the convergence rate in time, a sufficiently small spatial mesh size
is chosen such that the spatial discretization error can be relatively negligible. The time step is
with
. We show the numerical results at time
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
, to observe the spatial convergence rates. Take spatial size as
. Again, a quadratic finite element approximation for
and
is adopted, combined with a linear finite element approximation for p. Numerical results at
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
such that it holds that
for any
and
. 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
(4.5)
We fix
and
, and define the discrete energy function defined in Theorem 3.1. Meanwhile, we still adopt the quadratic elements for
and linear elements for p. Up to the final time
, we show the energy decaying curve as follows (Figure 1).
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.