A decoupled scheme with second-order temporal accuracy for magnetohydrodynamic equations

In this paper, we propose and analyze a temporally second-order accurate, fully discrete finite element method for the magnetohydrodynamic (MHD) equations. A modified Crank–Nicolson method is used to discretize the model and appropriate semi-implicit treatments are applied to the fluid convection term and two coupling terms. These semi-implicit approximations result in a linear system with variable coefficients for which the unique solvability can be proved theoretically. In addition, we use a decoupling method in the Stokes solver, which computes the intermediate velocity field based on the gradient of the pressure from the previous time level, and enforces the incompressibility constraint via the Helmholtz decomposition of the intermediate velocity field. This decoupling method greatly simplifies the computation of the whole MHD system. The energy stability of the scheme is theoretically proved, in which the decoupled Stokes solver needs to be analyzed in details. Optimal error estimates are provided for the proposed numerical scheme, with second-order temporal convergence and 𝒪 (h^r+1) spatial convergence, where r is the degree of the polynomial functions. Numerical examples are provided to illustrate the theoretical results.

Authors

• 70 publications
• 1 publication
• 1 publication
• 9 publications
10/29/2020

Analysis of Chorin-Type Projection Methods for the Stochastic Stokes Equations with General Multiplicative Noises

This paper is concerned with numerical analysis of two fully discrete Ch...
12/23/2021

Robust error bounds for the Navier-Stokes equations using implicit-explicit second order BDF method with variable steps

This paper studies fully discrete finite element approximations to the N...
04/14/2021

A second order ensemble method based on a blended BDF timestepping scheme for time dependent Navier-Stokes equations

We present a second order ensemble method based on a blended three-step ...
07/23/2020

An integral-based spectral method for inextensible slender fibers in Stokes flow

Every animal cell is filled with a cytoskeleton, a dynamic gel made of i...
12/21/2020

A novel structure preserving semi-implicit finite volume method for viscous and resistive magnetohydrodynamics

In this work we introduce a novel semi-implicit structure-preserving fin...
03/02/2021

A Variational Integrator for the Discrete Element Method

A novel implicit integration scheme for the Discrete Element Method (DEM...
09/13/2020

Second-order invariant domain preserving approximation of the compressible Navier–Stokes equations

We present a fully discrete approximation technique for the compressible...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

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

 μ∂tH+σ−1∇×(∇×H)−μ∇×(u×H)=σ−1∇×J, (1.1) ∂tu+u⋅∇u−νΔu+∇p=\boldmathf−μH×(∇×H), (1.2) ∇⋅u=0, (1.3)

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

 H|t=0=H0,u|t=0=u0,in Ω, (1.4) H×\boldmathn=0,  u=0,on ∂Ω×(0,T]. (1.5)

It is assumed that the initial data satisfies

 ∇⋅H0=∇⋅u0=0. (1.6)

By taking the divergence of (1.1), one can easily get , which together with the above divergence-free initial condition implies that

 ∇⋅H=0. (1.7)

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.

In this article, we propose and analyze a temporally second-order accurate, fully discrete decoupled finite element method for the MHD system (1.1)-(1.3

), 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

 Lp(Ω)=[Lp(Ω)]d,Wk,p(Ω)=[Wk,p(Ω)]d, W1,p0(Ω)=[W1,p0(Ω)]d,   H10(Ω)=W1,20(Ω), \ringH1(Ω)={v∈H1(Ω):v×\boldmathn|∂Ω=0},

where denotes the dimension of . As usual, the inner product of is denoted by .

With the above notations, it could be seen that the exact solution of (1.1)-(1.3) satisfies

 (μ∂tH,\boldmathw)+(σ−1∇×H,∇×\boldmathw)−(μu×H,∇×\boldmathw)=(σ−1∇×J,% \boldmathw), (2.1) (∂tu,v)+(ν∇u,∇v)+b(u,u,v)−(p,∇⋅v)=(\boldmathf,v)−(μH×(∇×H),v), (2.2) (∇⋅u,q)=0, (2.3)

for any test functions , where we have defined the trilinear form as

 b(u,v,\boldmathw) :=(u⋅∇v,\boldmathw)+12((∇⋅u)v,\boldmathw) (2.4) =12[(u⋅∇v,\boldmathw)−(u⋅∇\boldmathw,v)],∀u,v,\boldmathw∈H10(Ω),

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

 b(u,v,v)=0,∀u,v,\boldmathw∈H10(Ω).

The energy stability of the continuous system (2.1)-(2.3) could be obtained in a straightforward manner. By taking in (2.1)-(2.3) and adding the resulting equations together, we get

 μ2ddt∥H∥2L2+σ−1∥∇×H∥2L2+12ddt∥u∥2L2+ν∥∇u∥2L2 =(J,σ−1∇×H)+(\boldmathf,u) ≤14σ∥J∥2L2+σ−1∥∇×H∥2L2+14ε∥\boldmathf∥2L2+ε∥u∥2L2,

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

 μ2ddt∥H∥2L2+12ddt∥u∥2L2≤14σ∥J∥2L2+14ε∥% \boldmathf∥2L2. (2.5)

If the sources terms , we further get

 μ2ddt∥H∥2L2+12ddt∥u∥2L2≤0, (2.6)

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

 Xh={\boldmathlh∈H10(Ω):\boldmathlh|Kj∈Pr(Kj)}, Mh={qh∈L2(Ω):qh|Kj∈Pr−1(Kj), ∫Ωqhdx=0},

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

 Sh={\boldmathwh∈\ringH1(Ω):\boldmathwh|Kj∈Pr(Kj)}.

Let denote a uniform partition of the time interval , with a step size , and . For any sequences and , we define

 \widecheckvn+12:=34vn+1+14vn−1,¯¯¯vn+12:=12ˆvn+1+12vn,˜vn+12:=32vn−12vn−1.

Then, a fully discrete second-order decoupled FEM for the incompressible MHD equations (1.1)-(1.3) is formulated as: find such that

 −μ(¯¯¯un+12h×˜Hn+12h,∇×\boldmathwh)=σ−1(∇×Jn+12,\boldmathwh), (3.1) (ˆun+1h−unhτ,vh)+ν(∇¯¯¯un+12h,∇vh)+b(˜un+12h,¯¯¯un+12h,vh)−(pnh,∇⋅vh) +μ(˜Hn+12h×(∇×\widecheckHn+12h),vh)=(\boldmathfn+12,vh), (3.2) (un+1h−ˆun+1hτ,\boldmathlh)−12(pn+1h−pnh,∇⋅\boldmathlh)=0, (3.3) (∇⋅un+1h,qh)=0, (3.4)

for any and . Here we have added a stabilization term to (3.1). This is consistent with (2.1) in view of (1.7).

In this paper, it is assumed that the system (1.1)-(1.3) admits a unique solution satisfying

 ∥Httt∥L∞(0,T;L2)+∥Htt∥L∞(0,T;H1)+∥Ht∥L∞(0,T;Hr+1)+∥uttt∥L∞(0,T;L2) +∥utt∥L∞(0,T;H1)+∥ut∥L∞(0,T;Hr+1)+∥ptt∥L∞(0,T;L2)+∥pt∥L∞(0,T;Hr)≤K. (3.5)

Here, the subscripts of denote the partial derivative to variable .

Next, we present our main results, i.e., optimal error estimates for scheme (3.1)-(3.4), in the following theorem.

Theorem 3.1

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

 max2≤n≤N(∥Hnh−Hn∥L2+∥unh−un∥L2)≤C0(τ2+hr+1), (3.6) (τN∑n=2(∥∇×(Hnh−Hn)∥2L2+∥∇(¯¯¯un−12h−¯¯¯un−12)∥2L2))12≤C0(τ2+hr), (3.7)

where is a positive constant independent of and .

Remark 1

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.

Remark 2

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.

Remark 3

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.

Remark 4

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).

In the following subsection, we analyze the energy stability of scheme (3.1)-(3.4). In this paper, we denote by a generic positive constant and by a generic small positive constant, which are independent of , , , and .

3.2 Stability analysis of numerical scheme

In this subsection, we present the energy stability analysis for the numerical system (3.1)-(3.4). Here, we introduce a discrete version of the gradient operator, , defined as

 (vh,∇hqh)=−(∇⋅vh,qh),∀vh∈Xh,qh∈Mh. (3.8)

Through the definition of the discrete gradient operator , we can rewrite the equation (3.3) in the following equivalent form:

 un+1h−ˆun+1hτ+12∇h(pn+1h−pnh)=0. (3.9)

The abstract form (3.9) will be useful in the stability analysis of numerical scheme.

Theorem 3.2

The numerical solution to the fully discrete linearized FEM (3.1)-(3.4) satisfies the following energy stability estimate

 μ2τ(∥Hn+1h∥2L2−∥Hnh∥2L2)+μ8τ(∥Hn+1h−Hnh∥2L2−∥Hnh−Hn−1h∥2L2) (3.10) +12τ(∥un+1h∥2L2−∥unh∥2L2)+τ8(∥∇hpn+1h∥2L2−∥∇hpnh∥2L2) ≤C(∥Jn+12∥2L2+∥% \boldmathfn+12∥2L2),

for , where is a positive constant independent of and .

Proof

By taking in (3.1) and in (3.2), we get

 μ2τ(∥Hn+1h∥2L2−∥Hnh∥2L2)+μ8τ(∥Hn+1h−Hnh∥2L2−∥Hnh−Hn−1h∥2L2) +μ8τ∥Hn+1h−2Hnh+Hn−1h∥2L2+σ−1∥∇×\widecheckHn+12h∥2L2+σ−1∥∇⋅\widecheckHn+12h∥2L2 −μ(¯¯¯un+12h×˜Hn+12h,∇×\widecheckHn+12h) =σ−1(∇×Jn+12,\widecheckHn+12h), (3.11)

and

 12τ(∥ˆun+1h∥2L2−∥unh∥2L2)+ν∥∇¯¯¯un+12h∥2L2−(pnh,∇⋅¯¯¯un+12h) +μ(˜Hn+12h×(∇×\widecheckHn+12h),¯¯¯un+12h) =(\boldmathfn+12,¯¯¯un+12h), (3.12)

respectively, where we have used the fact that , and

 (Hn+1h−Hnhτ,% \boldmathwh) =12τ(∥Hn+1h∥2L2−∥Hnh∥2L2)+18τ(∥Hn+1h−Hnh∥2L2−∥Hnh−Hn−1h∥2L2) +18τ∥Hn+1h−2Hnh+Hn−1h∥2L2.

In turn, a substitution of in (3.3) yields

 12τ(∥un+1h∥2L2−∥ˆun+1h∥2L2+∥un+1h−ˆun+1h∥2L2)=0, (3.13)

where we have used the divergence-free condition for being .

Next, we choose in (3.3) and obtain

 −(∇⋅ˆun+1h,pnh)=τ4(∥∇hpn+1h∥2L2−∥∇hpnh∥2L2−∥∇h(pn+1h−pnh)∥2L2). (3.14)

Furthermore, we get the following result from (3.9)

 14∥∇h(pn+1h−pnh)∥2L2=1τ2∥un+1h−ˆun+1h∥2L2. (3.15)

 μ2τ(∥Hn+1h∥2L2−∥Hnh∥2L2)+μ8τ(∥Hn+1h−Hnh∥2L2−∥Hnh−Hn−1h∥2L2) (3.16) +μ8τ∥Hn+1h−2Hnh+Hn−1h∥2L2+σ−1∥∇×\widecheckHn+12h∥2L2+σ−1∥∇⋅\widecheckHn+12h∥2L2 +12τ(∥un+1h∥2L2−∥unh∥2L2)+ν∥∇¯¯¯un+12h∥2L2+τ8(∥∇hpn+1h∥2L2−∥∇hpnh∥2L2) ≤σ−1(∇×Jn+12,\widecheckHn+12h)+(\boldmathfn+12,¯¯¯un+12h).

For the right-hand side of (3.16), we can easily see that

 σ−1(∇×Jn+12,\widecheckHn+12h) =σ−1(Jn+12,∇×\widecheckHn+12h) ≤14σ∥Jn+12∥2L2+σ−1∥∇×\widecheckHn+12h∥2L2,

and

 (\boldmathfn+12,¯¯¯un+12h) ≤14ε∥\boldmathfn+12∥2L2+ε∥¯¯¯un+12h∥2L2≤14ε∥\boldmathfn+12∥2L2+ε∥∇¯¯¯un+12h∥2L2,

where is an arbitrarily small constant. Substituting the above estimates into (3.16), we get the desired result (3.10) immediately. This completes the proof of Theorem 3.2.

4 Optimal error estimates

We present the proof of the existence and uniqueness of numerical solution and the optimal error estimates (3.6)-(3.7) in Section 4.

4.1 Preliminary results

We introduce several types of projections. Let denote the projection which satisfies

 (v−Phv,qh)=0,v∈L2(Ω),∀qh∈Mh. (4.1)

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

 ν(∇(u−Rhu),∇vh)−(p−Rhp,∇⋅vh)=0, ∀vh∈Xh, (4.2) (∇⋅(u−Rhu),qh)=0, ∀qh∈Mh. (4.3)

We also introduce the Maxwell projection operator , by

 (∇×(H−ΠhH),∇×wh)+(∇⋅(H−ΠhH),∇⋅wh)=0,H∈\ringH1(Ω),∀wh∈Sh. (4.4)

For the above projections, the following estimates are recalled GR1987 ; Thomee2006 .

Lemma 1

The following estimates are valid for the projection, Stokes projection, and Maxwell projection:

 ∥Phv∥Wm,s≤C∥v∥Wm,s, (4.5) ∥v−Phv∥L2≤Chℓ+1∥v∥Hℓ+1, (4.6)

for , , , and

 ∥Rhu∥W1,s+∥Rhp∥Ls≤C(∥u∥W1,s+∥p∥Ls), (4.7) ∥u−Rhu