Error estimates

# Error estimates for structure-preserving discretization of the incompressible MHD system 111This material is based upon work supported in part by the US Department of Energy Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics program under Award Number DE-SC0014400 and by Beijing International Center for Mathematical Research of Peking University, China.

Yicong Ma Department of Mathematics,The Pennsylvania State University, University Park, PA 16802, USA Jinchao Xu Department of Mathematics,The Pennsylvania State University, University Park, PA 16802, USA  and  Guodong Zhang School of Mathematics and Statistics, Xi’an Jiaotong University, Xi’an 710049, China
###### Abstract.

In this paper, we carry out the error analysis for the structure-preserving discretization of the incompressible MHD system. This system, as a coupled system of Navier-Stokes equations and Maxwell’s equations, is nonlinear. We use its energy estimate and the underlying physical structure to facilitate the error analysis. Under certain CFL conditions, we prove the optimal order of convergence. To support the theoretical results, we also present numerical tests.

\reserveinserts

28

## 1. Introduction

An incompressible magnetohydrodynamic (MHD) system is a coupled partial differential equation system resulting from the incompressible Navier-Stokes equations and the (reduced) Maxwell’s equations. Assuming is a simply connected open-bounded domain with Lipschitz boundary, the model problem we consider is

 (1.1) ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩ut+(u⋅∇)u−R−1eΔu−sj×B+∇p=f,Bt+∇×E=0,j−R−1m∇×μ−1rB=0,σr(E+u×B)=j,∇⋅u=0.

The coefficients in this system are the Reynolds number , the magnetic Reynolds number , the coupling number , the relative electric conductivity , and the relative magnetic permeability . The initial conditions for this set of equations are

 u(x,0)=u0(x),B(x,0)=B0(x),∀x∈Ω,

and the boundary conditions are

 u=0,n×E=0,n⋅B=0,∀x∈∂Ω,t>0.

As discussed in the literature, the variables , and , once known, uniquely determine and . There are many different numerical methods to discretize MHD. We now briefly examine some existing literature on some of the numerical methods and their error analysis for two types of MHD systems: the stationary MHD system [12, 28] and the evolutionary MHD system [21, 26, 14].

For the stationary MHD system, Gunzburger, Meir and Peterson  propose a formulation with finite element discretization for the magnetic field, and analyze its well-posedness and convergent behavior. Schötzau , who also works on the stationary MHD system, proposes a new formulation with discretization for the magnetic field, and proves its well-posedness and the optimal order of convergence. There are also many other methods for stationary problems, for example, [10, 11, 29].

For the evolutionary MHD system, Prohl  studies the coupled and decoupled schemes based on conforming discretization of the magnetic field. He proves that the discrete solution converges to the weak solution under a strong Courant-Friedrichs-Lewy (CFL) condition; that is, ( stands for the time step size, and for the mesh size). And He  studies the MHD system on a regular domain with conforming discretization of the magnetic field. He proves an unconditional optimal order of convergence.

In this paper, we study the convergence property of a structure-preserving discretization presented in . This method is based on the mixed formulation , which comes from the idea of FEEC (finite element exterior calculus) [1, 2] and DEC (discrete exterior calculus) . and conforming finite element discretization are used for the electric field and the magnetic field respectively. The advantage of this approach is that the important Gauss’s law for magnetic field is preserved exactly on the discrete level. Moreover, the incompressible MHD system we focus on is a time-dependent nonlinear problem. Therefore, to conduct the error analysis, we work on an evolutionary nonlinear saddle point problem. Before approaching the detailed analysis, we briefly review the existing literatures for the error estimates of the (evolutionary) saddle point problems and nonlinear problems.

Abstract error estimates exist for standard (linear, non-evolutionary) saddle point systems . Optimal order of convergence is ensured by the well-posedness of the discretization system and the approximation property of the finite element space. For the evolutionary saddle point problem, Boffi and Gastaldi  build a general framework for the semi-discretization of the evolutionary (linear) saddle point problem and provide sufficient conditions for a good approximation in the natural functional spaces.

For nonlinear saddle point problems, no abstract error estimate framework can be found in the literature. But various techniques have been developed for specific problems. For example, Temann  discusses the theory and numerical methods for NS equations. Heywood and Rannacher [17, 18, 19, 20] discuss the stability and error estimates of both semi-discretization and full discretization schemes for the NS systems. He  study linearized implicit-explicit schemes for this model.

General error estimates exist for nonlinear parabolic and elliptic problems. Thomée et al. [27, 24, 22, 23, 9, 31] investigate the error estimates of nonlinear parabolic problems intensively. Xu  uses the priori estimate to derive the error estimates of a general nonlinear elliptic problem. Brezzi, Rappaz and Raviart [6, 7, 8] build an abstract theory for finite element approximation of nonlinear problems.

Due to the nonlinearity and the loss of coercivity of the MHD model, the error estimate becomes difficult. To estimate the error of nonlinear problems, we usually need to prove that the norm (or a stronger norm) of the numerical solution is bounded. Generally, there are two ways to obtain this bound, one is using the mathematical induction method [13, 15], the other is introducing a semi-discrete problem [B.Li;W.Sun2013a, B.Li;W.Sun2013b]. Moreover, due to the loss of coercivity, we cannot use Cea’s Lemma to derive the error estimates directly.

In our analysis, we take advantage of the energy estimate of the structure-preserving discretization instead of estimating the norm of the numerical solution. We prove the unconditional error estimates for the velocity , the magnetic field , and the volume current density . And under certain constrains on the time-step size, we derive error estimates for the electric field and the pressure . Numerical tests support the theoretical results.

We organize this paper as follows. In §2, we introduce useful notation for our analysis. In  §3, we go over the discretization schemes and their energy estimates. We carry out detailed error estimates in §4, and present numerical experiments to demonstrate the optimal order of convergence in §5.

## 2. Magnetohydrodynamics model

In this section, we introduce some notation which follows mostly . We first define the usual inner product

 (u,v)=∫Ωu⋅vdx,

and the norm

 ∥u∥=(∫Ω|u|2dx)1/2.

For the sake of simplicity, we write both and as .

Given a linear operator , we define

 H(D,Ω)={v∈L2(Ω), Dv∈L2(Ω)},

and

 H0(D,Ω)={v∈H(D,Ω), tDv=0 on ∂Ω},

where is the trace operator defined by

We define

 L20(Ω):={v∈L2(Ω), ∫Ωv=0}.

When , we typically write instead of , and instead of . In the analysis, we also use the spaces and with norms

 ∥v∥m,p=⎛⎝∑0≤|α|≤m∫Ω|Dαv|p⎞⎠1/p,∥v∥0,∞=esssupx∈Ω|v|,∥v∥−1=supϕ∈H10(Ω)(v,ϕ)∥∇ϕ∥.

Another useful norm in the analysis is the discrete norm

 \VERTu\VERT2m,∗=km∑n=1∥un∥2∗,

where is the time step size (as is the in the following context). For example, stands for the discrete norm and for the discrete norm.

Next, we introduce some useful function spaces in the discretization.

 X=H10(Ω)3×H0(div,Ω)×H0(curl,Ω),Q=L20(Ω), V=H10(Ω)3, Vd=H0(div,Ω), Vc=H0(curl,Ω).

We use (, , or ) to denote the dual space of , and to denote the finite element space of . The divergence-free subspace of is defined as

 V0={v∈V, ∇⋅v=0}.

For the sake of simplicity, we assume that in the analysis. The Hilbert spaces and are equipped with norms and , which are defined as

 ∥ξ∥2X=∥v∥21+∥C∥2div+∥F∥2curl,∀ξ=(v,C,F)∈X, ∥q∥2Q=∥q∥2,∀q∈Q.

Here,

 ∥v∥21=∥v∥2+∥∇v∥2, ∀v∈H1(Ω), ∥C∥2div=∥C∥2+∥∇⋅C∥2, ∀C∈H(div,Ω), ∥F∥2curl=∥F∥2+∥∇×F∥2, ∀F∈H(curl,Ω).

We also use Sobolev space and , which are defined as

 Hr(div,Ω)={v∈Hr(Ω), ∇⋅v∈Hr(Ω)}, Hr(curl,Ω)={v∈Hr(Ω), ∇×v∈Hr(Ω)}.

The corresponding norms are denoted by and , which are defined as

 ∥C∥2r,div=∥C∥2r,2+∥∇⋅C∥2r,2, ∥F∥2r,curl=∥F∥2r,2+∥∇×F∥2r,2.

To facilitate the analysis, we also introduce a tri-linear form of , namely,

 (2.1) c(ϕ,u,v)=12(ϕ⋅∇u,v)−12(ϕ⋅∇v,u).

Based on the above notations, the variational formulation for system (1.1) is: find and such that for any and ,

 (2.2) ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩(ut,v)+c(u,u,v)+R−1e(∇u,∇v)−s(j×B,v)−(p,∇⋅v)=(f,v),α(Bt,C)+α(∇×E,C)=0,s(j,F)−α(B,∇×F)=0,(∇⋅u,q)=0,

where , and .

## 3. Finite element discretization

In this section, we briefly go over the finite element discretization of (1.1). For the temporal discretization, we use the backward Euler method. For the spacial discretization, we recall the formulation of both nonlinear and linearized discretization here. These discretization formulations are reasonable in the sense that they inherit the energy estimate from the continuous level. At the end of this section, we go over their energy estimates.

###### Algorithm 3.1.

Find and such that for any and ,

 (3.1) ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩(¯∂unh,vh)+c(un−1h,unh,vh)+R−1e(∇unh,∇vh)−s(jnh×Bn−1h,vh)−(pnh,∇⋅vh)=(fn,v),α(¯∂Bnh,Ch)+α(∇×Enh,Ch)=0,s(jnh,Fh)−α(Bnh,∇×Fh)=0,(∇⋅unh,qh)=0,

where , , and .

The above formulation uses linearization as a discretization scheme. In fact, we can discretize the nonlinear system directly and solve the nonlinear equation by Picard or Newton iteration.

###### Algorithm 3.2.

Find and such that for any and ,

 (3.2) ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩(¯∂unh,vh)+c(unh,unh,vh)+R−1e(∇unh,∇vh)−s(jnh×Bnh,vh)−(pnh,∇⋅vh)=(fn,v),α(¯∂Bnh,Ch)+α(∇×Enh,Ch)=0,s(jnh,Fh)−α(Bnh,∇×Fh)=0,(∇⋅unh,qh)=0,

where .

As mentioned before, these above formulations admit desirable energy estimates. For the sake of completeness, we cite some of these estimates, which are established in .

###### Theorem 3.1 (Energy estimates).

For any , and that satisfies (3.2), the following energy estimates hold

 max0≤n≤m(∥unh∥2+α∥Bnh∥2)+R−1e\VERT∇uh\VERT2m,0+2s\VERTjh\VERT2m,0≤∥u0∥2+α∥B0∥2+Re\VERTf\VERT2m,−1,

where .

Here, and depend on the initial data. A similar energy estimate holds for (3.1).

## 4. Error estimates

Before starting the detailed analysis, first we recall Gronwall’s inequality , which is an important tool in our analysis.

###### Theorem 4.1 (Gronwall’s inequality ).

Let , , and , , , , for integers , be non-negative numbers such that

 an+kn∑i=0bi≤kn∑i=0γiai+kn∑i=0ci+B,∀ n≥0.

Suppose that (for all ), and set , then

 an+kn∑i=0bi≤exp(kn∑i=0σiγi)(kn∑i=0ci+B),∀ n≥0.

We choose to be the -th order polynomial space, the -th order Raviart-Thomas elements, the -th order Nédélec element, and the -th order polynomial space.

Define , , and

 a(ξ,ξ,η)=c(u,u,v)+R−1e(∇u,∇v)−s(j×B,v)+α(∇×(j−u×B),C) (4.1) +s(j,F)−α(B,∇×F),∀ξ, η∈X, (4.2) b(η,q)=−(∇⋅v,q),∀η∈X, q∈Q.

Therefore, we can write the MHD system (2.2) in the form of a saddle point problem. That is, find such that for any ,

 (4.3) {(Aξt,η)+a(ξ,ξ,η)+b(η,p)=⟨h,η⟩,b(ξ,q)=0,

where , and . Additionally, we can write (3.2) as: find such that for any ,

 (4.4) {(A¯∂ξnh,ηh)+a(ξnh,ξnh,ηh)+b(ηh,pnh)=⟨hn,ηh⟩,b(ξnh,qh)=0.

Before giving the detailed error estimates, we define the projections of first. Assume that is the Stokes projection of . That is, for any given , we define such that satisfies

 {R−1e(∇ρnu,∇vh)−(ρnp,∇⋅vh)=0,∀vh∈Vh,(∇⋅ρnu,qh)=0,∀qh∈Qh.

We choose as the projection of , and as the canonical interpolation of .

Define

 enw=wn−wnh,ρnw=wn−ˆwnh,θnw=ˆwnh−wnh,w=u, B, E, p, j.

By the definitions of and , we have

 enj =En+un×Bn−(Enh+unh×Bnh) =enE+un×Bn−unh×Bnh =enE+un×enB+enu×Bnh.

Similarly, we define

 ρnj=ρnE+un×ρnB+ρnu×Bnh,θnj=θnE+un×θnB+θnu×Bnh.

It follows that . Moreover, we define

 enξ=(enu,enB,enj),ρnξ=(ρnu,ρnB,ρnj),θnξ=(θnu,θnB,θnj).

For simplicity, we denote

 ∥ρnξ∥2L2×L2×curl=∥ρnu∥2+∥ρnB∥2+∥ρnE∥2curl, ∥ρnξ∥2H1×L2×curl=∥∇ρnu∥2+∥ρnB∥2+∥ρnE∥2curl, ∥A1(¯∂ξn−ξnt)∥2=∥¯∂un−unt∥2+∥¯∂Bn−Bnt∥2, ∥A1¯∂ρnξ∥2=∥¯∂ρnu∥2+∥¯∂ρnB∥2, ρ0=max1≤n≤N{∥A1(¯∂ξn−ξnt)∥2+∥A1¯∂ρnξ∥2+∥ρnξ∥2L2×L2×curl}, ρ1=max1≤n≤N{∥A1(¯∂ξn−ξnt)∥2+∥A1¯∂ρnξ∥2+∥ρnξ∥2H1×L2×curl},

where Noticing that

 c(un,un,vh)−c(unh,unh,vh)=c(enu,un,vh)+c(unh,enu,vh), (jn×Bn,vh)−(jnh×Bnh,vh)=(jn×enB,vh)+(enj×Bnh,vh),

we can rewrite the error equation as

 (4.5) ⎧⎪ ⎪⎨⎪ ⎪⎩(A¯∂θnξ,ηh)+ˆa(ξn,ξnh,θnξ,ηh)+b(ηh,θnp)=(A¯∂ξn−Aξnt,ηh)−(A¯∂ρnξ,ηh)−ˆa(ξn,ξnh,ρnξ,ηh)−b(ηh,ρnp),∀ηh∈Xh,b(θnξ,qh)=−b(ρnξ,qh),∀qh∈Qh,

where

 ˆa(ξn,ξnh,θnξ,ηh) =c(θnu,un,vh)+c(unh,θnu,vh)+R−1e(∇θnu,∇vh)−s(jn×θnB,vh) −s(θnj×Bnh,vh)+α(∇×(θnj−un×θnB−θnu×Bnh),Ch) +s(θnj,Fh)−α(θnB,∇×Fh).

### 4.1. Main results

We summarize main results of this paper for error estimates of (3.2) in the following theorem.

###### Theorem 4.2.

For any fixed time step such that , if is the solution to (4.3), and is the solution to (3.2), the following estimates hold:

1. There exists a constant , which only depends on , , and the computation domain, such that

 (4.6) ∥um−umh∥2+α∥Bm−Bmh∥2+R−1e\VERT∇θu\VERT2m,0+s\VERTj−jh\VERT2m,0≤Cρ0,

when the time step size is sufficiently small. And there also holds

 (4.7) k\VERTp−ph\VERT2m,0≤C(max1≤n≤N∥ρnp∥2+ρ0+h−1ρ20).
2. There exists a constant , which only depends on , , , and the computation domain, such that

 (4.8) \VERTE−Eh\VERT2m,0+k\VERT∇×E−∇×Eh\VERT2m,0≤Cρ0(1+h−1ρ0).
3. There exists a constant , which only depends on , , , and the computation domain, such that

 (4.9) \VERTp−ph\VERT2m,0≤C(ρ0+ρ20h−3+ρ1+ρ1ρ0h−1+max1≤n≤N∥ρnp∥2).

By similar arguments similar to Theorem 4.2, we can get the error estimates of the Picard linearization scheme (3.1) as follows.

###### Theorem 4.3.

For any fixed time step such that , we have the following error estimates of (3.1):

1. There exists a constant , only depending on the exact solution, such that

 ∥um−umh∥2+α∥Bm−Bmh∥2+R−1e\VERT∇θu\VERT2m,0+s\VERTj−jh\VERT2m,0≤Cρ0,

when the time step size is sufficiently small. And there also holds

 k\VERTp−ph\VERT2m,0≤C(max1≤n≤N∥ρnp∥2+ρ0+h−1ρ20).
2. There exists a constant only depending on exact solution such that

 \VERTE−Eh\VERT2m,0+k\VERT∇×E−∇×Eh\VERT2m,0≤Cρ0(1+h−1ρ0),

when the time step size is sufficiently small.

3. There exists a constant only depending on exact solution such that

 \VERTp−ph\VERT2m,0≤C(ρ0+ρ20h−3+ρ1+ρ1ρ0h−1+max1≤n≤N∥ρnp∥2).

when the time step size is sufficiently small.

As the proof of the above theorem is similar to that of theorem 4.2, we omitted its details in this paper.

### 4.2. Proof of Theorem 4.2

The basic idea of our proof is to verify that satisfies the Gårding condition, and the truncation error is bounded by some norm of . Then the conclusion follows by Gronwall’s inequality.

First of all, we prove that satisfies the Gårding condition.

###### Lemma 4.1.

The sum of bilinear form satisfies the Gårding condition. That is, for any , there exists , such that

 ˆa(ξn,ξnh,θnξ,ηh)≥β1∥∇θnu∥2+β1∥θnj∥2−β0∥θnu∥2−β0∥θnB∥2,

where and are positive constants that only depend on , , and the computation domain.

###### Proof.

By definition, we know that . Since , we choose . Noticing that , we get

 ˆa(ξn,ξnh,θnξ,ηh) = c(θnu,un,θnu)+c(unh,θnu,θnu)+R−1e(∇θnu,∇θnu)−s(jn×θnB,θnu) −s(θnj×Bnh,θnu)+s(θnj,θnj−un×