The magnetohydrodynamic equations have been widely applied into metallurgy and liquid-metal processing, and the numerical solutions are of great significance in practical scientific and engineering applications; see app1 and app2 . Such an MHD system could be formulated as
over , where is a bounded and convex polyhedral domain in (polygonal domain in ). In the above system, , and denote the velocity field, the magnetic filed, and the pressure, respectively; and are the given source terms ( denotes a scalar function in ); denotes the magnetic Reynolds number, the viscosity of the fluid , and , where is the Hartman number. The initial data and boundary conditions are given by
It is assumed that the initial data satisfies
By taking the divergence of (1.1), one can easily get , which together with the above divergence-free initial condition implies that
The existence and uniqueness of the weak solution for this problem has been theoretically proved in Max1991 ; PDE1 . More regularity analysis of the MHD system could be referred in PDE7 ; PDE5 ; PDE6 ; PDE4 , etc.
There have been many existing works on the numerical approximations for the incompressible MHD system. In bounded and convex domains, the solutions of the MHD model are generally in ( denotes the dimension of ) and therefore people often use -conforming finite element methods (FEMs) to solve the MHD equations numerically. For example, Gunzburger, Meir, and Peterson Max1991 used -conforming FEMs for solving the stationary incompressible MHD equations with an optimal error estimate being established. Later, He developed -conforming FEMs in FEMr3 for solving the time-dependent MHD equations and proved error estimates of the numerical scheme. More works on -conforming FEMs for the MHD equations can be found in FEMr1 ; Gerbeau2000A ; FEMr2 ; FEMr4 ; main .
While the spatial approximation has always been important, the temporal discretization also plays a significant role for solving the MHD system. There have been quite a few existing stability and convergence analyses for the first-order temporally accurate numerical schemes first ; FEMr3 ; FDM1 ; ratecon ; FDM2 . In most of these works, the stability and convergence analyses have been based on a Stokes solver at each time step, i.e., the computation of the pressure gradient has to be implemented with the incompressibility constraint being enforced, which in turn leads to a non-symmetric linear system, and the computation costs turn out to be extremely expensive. To overcome this difficulty, some “decoupled” techniques have been introduced. In shen2016A , Zhao, Yang, Shen, and Wang dealt with a binary hydrodynamic phase field model of mixtures of nematic liquid crystals and viscous fluids by designing a decoupled semi-discrete scheme, which is linear, first-order accurate in time, and unconditionally energy stable. In particular, a pressure-correction scheme Guermond2006 was used so that the pressure could be explicitly updated in the velocity equation by introducing an intermediate function and thus two sub-systems are generated. In shen2015 , Liu, Shen, and Yang proposed a first-order decoupled scheme for a phase-field model of two-phase incompressible flows with variable density based on a “pressure-stabilized” formulation, which treats the pressure term explicitly in the velocity field equation, and only requires a Poisson solver to update the pressure. These works have mainly focused on the design of energy-preserving schemes without presenting the convergence analysis. Meanwhile, the first-order temporal accuracy may not be sufficient in the practical computations of the MHD system, and therefore higher-order temporal numerical approximations have been highly desired. In FEM3 , a second-order conditionally stable backward differential formulation (BDF) algorithm was presented for a reduced MHD model at small magnetic Reynolds number, in which the coupling terms are explicitly updated, and other terms are implicitly computed. Another kind of BDF scheme with second-order accuracy in time was proposed in Heister2016 , and the method was proved to be unconditionally stable and convergent at an optimal order. In second2 , a second-order scheme with Newton treatment of the nonlinear terms was proposed, where the unconditional stability and optimal error estimates were obtained. Recently, a fully discrete Crank–Nicolson (CN) scheme was studied in FEMr5 , where the unconditional energy stability and convergence (without error estimates) were proved.
), and the following properties will be theoretically established: unique solvability, unconditional energy stability, and optimal rate convergence analysis. In particular, a modified Crank–Nicolson method with an implicit Adams–Moulton interpolation in the form of, instead of the standard Crank–Nicolson approximation, is applied to discretize the magnetic diffusion term. Such a technique leads to a stronger stability property of the numerical scheme, as will be demonstrated in the subsequent analysis. An explicit Adams–Bashforth extrapolation formula is used in the approximation of the nonlinear convection term. In addition, an intermediate velocity function is introduced in the numerical scheme, and its computation is based on the pressure gradient at the previous time step. After solving the intermediate velocity field, we decompose it into the divergence-free subspace by using the Helmholtz decomposition. This yields the velocity field at the same time level. Meanwhile, certain semi-implicit treatments have been adopted to the two coupled terms, and these treatments with the decoupled Stokes solver lead to a linear numerical system with variable coefficients. The energy stability is theoretically proved for such a numerical system, in which a rigorous estimate of the decoupled Stokes solver plays an important role. Then, the unique solvability of the proposed scheme follows immediately by the fact that the corresponding homogeneous equations only admit zero solutions. Furthermore, optimal error estimate of is established in -norm, with being the degree of the polynomial functions, where and denote the time step size and space mesh size, respectively. To our knowledge, this is the first work on preserving such four theoretical properties for a second-order numerical scheme of the MHD system: decoupled Stokes solver, unique solvability, unconditional energy stability and optimal rate convergence analysis. In more details, the semi-implicit treatment of the fluid convection part and the two coupled nonlinear terms have played an important role in the stability and convergence analysis.
This paper is organized as follows. In Section 2, a variational formulation and some preliminary results are reviewed. The fully discrete finite element scheme is introduced in Section 3, and its unconditional energy stability is established in details. Section 4 provides the rigorous proof of the unique solvability and optimal error estimates. Several numerical examples are presented in Section 5. Finally, some concluding remarks are provided in Section 6.
2 Variational formulation and stability analysis
For and , let be the conventional Sobolev space of functions defined on , with abbreviations and . Then, we denote by the space of functions in with zero traces on the boundary , and denote
. The corresponding vector-valued spaces are
where denotes the dimension of . As usual, the inner product of is denoted by .
for any test functions , where we have defined the trilinear form as
and denotes the Euclidean scalar product in . Notice that the trilinear form
is skew-symmetric with respect to its last two arguments, so that we further have
where is an arbitrary constant. Due to the zero boundary condition of in (1.5), we have . Since can be arbitrarily small, we obtain the following energy estimate
If the sources terms , we further get
which implies the total energy is decaying.
3 Numerical method and theoretical results
3.1 Numerical method
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 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 on for all and . To approximate the magnetic field , we introduce the finite element space defined by
Let denote a uniform partition of the time interval , with a step size , and . For any sequences and , we define
Here, the subscripts of denote the partial derivative to variable .
Suppose that the system (1.1)-(1.3) has a unique solution satisfying (3.5). Then there exist positive constants and such that when , , and , the fully discrete decoupled FEM system (3.1)-(3.4) admits a unique solution , , which satisfies that
where is a positive constant independent of and .
In (3.1)-(3.4), we have used a modified Crank–Nicolson method for temporal discretization, where the term is approximated by . This enables us to obtain error estimates for the term at certain time steps, instead of an average of those at two consecutive time levels; see (3.7). Such a modified Crank–Nicolson scheme has been extensively applied to various gradient flow models chen14 ; cheng16a ; diegel16 ; guo16 . An application of this approach to the incompressible MHD system is reported in this work, for the first time.
In (3.1), we have added a stabilization term to validate the coercivity of the magnetic equation, with which, optimal error estimates for the magnetic field in energy-norm can be proved.
Another feature of the proposed numerical scheme (3.1)-(3.4) is associated with its decoupled nature in the Stokes solver. Motivated by the second order incremental pressure projection method of van Kan type vankan1986 , we introduce an intermediate velocity to decouple the problem, and thus build two systems and both of them consist of two unknowns. More precisely, we first obtain and through (3.1)-(3.2), while treating the gradient of pressure explicitly. Then, we substitute into (3.3)-(3.4), so that and could be efficiently computed via solving a Darcy problem. In comparison with the classical coupled solver that the full system contains three unknowns , and , which have to be solved simultaneously, such a decoupled approach will greatly improve the efficiency of the numerical scheme.
There have been extensive analyses of decoupled numerical schemes for incompressible Navier–Stokes equations; see the pioneering works of A. Chorin Chorin1968 , R. Temam Temam1969 , and many other related studies BCG1989 ; E1995 ; E2002 ; Kim1985 ; STWW2003 ; Wang2000 , etc. Most existing works focused on the pure fluid model, while this article is the first one to report a decoupled temporally second-order accurate numerical scheme for the incompressible MHD system, with energy stability and optimal error estimates to be established later.
It is noted that the numerical solutions at two previous time levels are needed for the implementation of (3.1)-(3.4). The starting values at time steps and are assumed to be given and satisfy the estimates (3.6)-(3.7). An example of constructing the numerical scheme for starting values at is an application of the backward Euler method at the first time step, which satisfies the optimal convergence given in (3.6)-(3.7).
3.2 Stability analysis of numerical scheme
Through the definition of the discrete gradient operator , we can rewrite the equation (3.3) in the following equivalent form:
The abstract form (3.9) will be useful in the stability analysis of numerical scheme.
respectively, where we have used the fact that , and
In turn, a substitution of in (3.3) yields
where we have used the divergence-free condition for being .
Next, we choose in (3.3) and obtain
Furthermore, we get the following result from (3.9)
4 Optimal error estimates
4.1 Preliminary results
We introduce several types of projections. Let denote the projection which satisfies
For the sake of brevity, if (defined in (4.1)) is a vector function in , we still use to denote the projection over the finite element space . Furthermore, let denote the Stokes projection of satisfying
We also introduce the Maxwell projection operator , by
The following estimates are valid for the projection, Stokes projection, and Maxwell projection:
for , , , and