Analysis of coupling iterations based on the finite element method for stationary magnetohydrodynamics on a general domain

Computers and Mathematics with Applications 68 (2014) 770–788

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

Analysis of coupling iterations based on the finite element method for stationary magnetohydrodynamics on a general domain✩ Guo-Dong Zhang, Yinnian He ∗ , Di Yang School of Mathematics and Statistics, State Key Laboratory of Multiphase Flow in Power Engineering, Xi’an Jiaotong University, Xi’an 710049, People’s Republic of China

article

info

Article history: Received 7 August 2013 Received in revised form 21 May 2014 Accepted 26 July 2014 Available online 22 August 2014 Keywords: Magnetohydrodynamics Coupling iteration Finite element methods Singular solution

abstract This paper considers finite element approximation and three coupled type iterations of stationary magnetohydrodynamics (MHD) equations on a general Lipschitz domain. An additional Lagrange multiplier r is introduced related to divergence free of magnetic field b and the b is analyzed in H (curl; Ω ) space, which is proposed in Schötzau (2004). In finite element discretization, the hydrodynamic unknowns are approximated by stable finite element pairs, and the magnetic unknown is discretized by curl-conforming Nédélec element spaces. The well-posedness of this formula and the optimal error estimate are provided. Based on this, for numerical implementation of this scheme, we propose and discuss three coupled type iterative methods which are stable and convergent under different conditions. Specifically, Iteration I is stable and convergent under strong condition. Iteration II is stable and convergent under weaker condition. Iteration III is unconditionally stable and convergent under the weakest condition. Finally, some numerical tests confirm our theoretical analysis. © 2014 Elsevier Ltd. All rights reserved.

1. Introduction Magnetohydrodynamics (MHD) is the theory of macroscopic interaction of electrically conducting fluid and electromagnetic fields. Many interesting MHD-flows involve a viscous, incompressible, electrically conducting fluid that interacts with an electromagnetic field. It is important in many engineering and physics fields, such as sustained plasma confinement for controlled thermonuclear fusion, liquid metal cooling of nuclear reactors, electromagnetic casting of metals and some applications in geophysics and astronomy [1,2]. The governing partial differential equations are obtained by coupling the incompressible Navier–Stokes equations with Maxwell equations. In the magnetic problem, the coupling with the hydrodynamic problem comes from the convective term. In the Navier–Stokes equations, the coupling with the magnetic problem comes from Lorentz force. In recent years, there are many papers studying the MHD equations from mathematical theory viewpoint, such as considering well-posedness, regularity of weak solutions and long-time behavior of solutions, see [3–14]. There are some studies devoting to numerical approximation for the MHD problems, see [15–27] and references therein. When the domain Ω ⊂ Rd is a convex polygon (polyhedron) or has C 1,1 boundary, the magnetic solution of MHD equations has H 1 regularity. The

✩ This work is supported by the NSFs of China (Nos. 11271298 and 11362021).

∗

Corresponding author. E-mail address: [email protected] (Y. He).

http://dx.doi.org/10.1016/j.camwa.2014.07.025 0898-1221/© 2014 Elsevier Ltd. All rights reserved.

G.-D. Zhang et al. / Computers and Mathematics with Applications 68 (2014) 770–788

771

functional setting of MHD equations is based on standard Sobolev spaces H1 (Ω ) and L2 (Ω ). In this situation, the magnetic variables of MHD can be discretized by H 1 -conforming finite element spaces, see [17,19,20,22,25]. However, when the domain Ω is a non-convex polygon (polyhedron), the magnetic variables may have regularity below H 1 , and the continuous nodal finite element spaces may fail to approximate the singular magnetic solution. There are some studies to deal with the non-regular case. For example, in [15], a new stabilized finite element formulation was proposed, which can always converge to the physical solution, even for singular ones. In [23], weighted regularization approach was applied to cope with the magnetic field singularity. In [21,27], through introducing an additional Lagrange multiplier r related to ∇ · b = 0, the magnetic unknown is analyzed in H (curl; Ω ) space and discretized by curl-conforming Nédélec elements which are only tangential continuous across inter-elements’ faces. This scheme makes this approach feasible in the situations of singular magnetic fields. For more details, please see [28–30]. In [27], the well-posedness of the continuous and discrete problem is obtained under small data assumptions. But, the author failed to prove Proposition 3.2 and Corollary 3.1 in this paper, which are indispensable to give an error estimate. In [21], the authors considered the velocity in divergence-conforming Brezzi–Douglas–Marini (BDM) finite element space, the magnet in the curl-conforming Nédélec elements. However, they did not present an optimal error estimate. Motivated by the previous work, in this article, we use the formula introduced in [27] to deal with the magnetic singularity on the general Lipschitz domain. We first give the accurate condition of well-posedness for continuous and discrete problem by a simple method. Then, we provide the optimal error estimate for this finite element scheme. Based on this, for numerical implementation of this scheme, we introduce three coupling iterative methods and study their stabilities and convergence. 5 , Iteration II is stable and convergent. Specifically, when 0 < σ ≤ 41 , Iteration I is stable and convergent. When 0 < σ ≤ 11 Iteration III is unconditionally stable and convergent in the case of 0 < σ < 1, where the σ is defined in Section 4. Finally, some numerical tests verify our theoretical results. The rest of paper is organized as follows. In Section 2, we give some preliminaries of MHD equations. In Section 3, the well-posedness of the continuous problem is discussed. In Section 4, the well-posedness of the discrete problem and the optimal error estimate are provided. In Section 5, we introduce three coupling iterations and study their stabilities and convergence under different conditions. In Section 6, we present some numerical experiments to validate the theoretical results. In Section 7, we end with some concluding remarks. 2. Function spaces and preliminary results Let Ω be a bounded Lipschitz domain in Rd , d = 2, 3. We consider stationary incompressible MHD problems

−Re−1 1u + u · ∇ u + ∇ p − S curl b × b = f, in Ω Rm−1 S curl(curl b) − S curl(u × b) − ∇ r = g, in Ω ∇ · u = 0, in Ω ∇ · b = 0, in Ω u = 0, b × n = 0, r = 0, on ∂ Ω where n denotes the outward normal unit vector on ∂ Ω , Re is hydrodynamic Reynolds number, Rm is the magnetic Reynolds number and S > 0 is the coupling coefficient. The function f ∈ H −1 (Ω )d , g ∈ L2 (Ω )d are given source terms, and ∇ · g = 0. The scalar function r, named virtual magnetic pressure, is the Lagrange multiplier associated to the constraint ∇ · b = 0. In fact, by taking divergence of the magnetic equation, we obtain −1r = 0, which combines the boundary condition r = 0 on ∂ Ω yields r = 0 in Ω . In this article, vector-valued functions and function spaces will be denoted in bold-face. Lr (Ω ) denotes the usual Lebesgue space on Ω with the norm ∥ · ∥Lr , 1 ≤ r ≤ ∞. The norm in L2 (Ω ) is denoted by ∥ · ∥0 . For a nonnegative real number s, we write H s (Ω ) for W s,2 (Ω ), and the corresponding norm is ∥φ∥s . We denote the vector and scalar function spaces that will be used as: H1 (Ω ) = H 1 (Ω )d , H10 (Ω ) = w ∈ H1 (Ω ) : w|∂ Ω = 0 ,



L20 (Ω ) =





q ∈ L2 (Ω ) :



 Ω

qdx = 0 ,

H (curl; Ω ) = b ∈ L2 (Ω )d : curl b ∈ L2 (Ω )d ,





H0 (curl; Ω ) = b ∈ H (curl; Ω ) : b × n |∂ Ω = 0 ,





H (div; Ω ) = b ∈ L2 (Ω )d : div b ∈ L2 (Ω ) ,





H0 (div; Ω ) = b ∈ H (div; Ω ) : b · n |∂ Ω = 0 .





⟨, ⟩ represents the duality pairing between H−1 (Ω ) and H10 (Ω ), (,) the inner product in L2 (Ω )d . For mathematical setting of the MHD equations, we introduce the following trilinear and bilinear forms, for w, u, v, ∈ H10 (Ω ), b, c ∈ H0 (curl, Ω ),

772

G.-D. Zhang et al. / Computers and Mathematics with Applications 68 (2014) 770–788

d ∈ X, q ∈ L2 (Ω ), s ∈ H 1 (Ω ):



as (u, v) = Re−1

∇ u∇ vdx,  curl b · curl cdx, am (b, c) = S Rm−1  Ω  1 1 c0 (w, u, v) = w∇ u · vdx − w∇ v · udx, 2 Ω 2 Ω  c1 (d, v, b) = S curl b × d · vdx,  Ω  bs (q, v) = − q∇ · vdx, bm ( s, c ) = − ∇ s · cdx, Ω

Ω

Ω

⟨L, (v, c)⟩ = ⟨f, v⟩ + (g, c), A(u, b; v, c) = as (u, v) + am (b, c), C (w, d; u, b; v, c) = c0 (w, u, v) − c1 (d, v, b) + c1 (d, u, c), B(q, s; v, c) = bs (q, v) + bm (s, c). It is easy to see the trilinear term C (·, ·; ·, ·; ·, ·) is skew-symmetric about the later two variables, i.e., C (w, d; u, b; u, b) = 0,

∀w, u, ∈ H1 (Ω ), b ∈ H (curl, Ω ), d ∈ X.

(2.1)

(Ω ) × H0 (curl; Ω ) and (Ω ) × (Ω ) are defined as  1  1 ∥(p, r )∥ = ∥p∥20 + S −1 ∥∇ r ∥20 2 , ∥(u, b)∥ = ∥∇ u∥20 + S ∥b∥2curl 2 ,  1 respectively, here, the curl-norm ∥b∥curl = ∥b∥20 + ∥curl b∥20 2 . We denote   Υ = (u, b) ∈ H10 (Ω ) × H0 (curl; Ω ) : B(q, s; u, b) = 0, ∀(q, s) ∈ L20 (Ω ) × H10 (Ω ) = Z × X,

The energy norms on product spaces

H10

L20

H01

where, Z = {u ∈ H10 (Ω ) : ∇ · u = 0}, X = {b ∈ H0 (curl; Ω ) : ∇ · b = 0}. The above relation can be deduced using integration by parts, also refer [27]. The following important inequalities are known, see Theorem 3.50 and Corollary 3.51 in [28], Part I Theorem 1.3 in [31]. 1 2 ∥b∥2curl ≤ λ− 0 ∥curl b∥0 ,

∥u∥L3 ≤ c ∥∇ u∥0 ,

∀ b ∈ X,

(2.2)

∥u∥L6 ≤ λ1 ∥∇ u∥0 ,

∥b∥L3 ≤ λ2 ∥b∥curl ,

∀u∈

H10

(Ω ),

∀ b ∈ X,

(2.3) (2.4)

here, the positive constants c , λi only depend on Ω , i = 0, 1, 2. With these notations, we have some estimates: Lemma 2.1. The following estimates for bilinear and trilinear terms hold:

|A(u, b; v, c)| ≤ max{Re−1 , Rm−1 }∥(u, b)∥∥(v, c)∥,

∀u, v ∈ H10 (Ω ), b, c ∈ H0 (curl; Ω ),

A(u, b; u, b) ≥ min{Re

H10

|B(q, s; v, c)| ≤

√

−1

, Rm λ0 }∥(u, b)∥ , −1

d∥(v, c)∥∥(q, s)∥,

c0 (w, u, v)

sup w,u,v∈H10 (Ω )

2

∥∇ w∥0 ∥∇ u∥0 ∥∇ v∥0

sup d∈X,∀v∈H10 (Ω ),b∈H0 (curl;Ω )

∀u ∈

(Ω ), b ∈ X,

∀v ∈ H10 (Ω ), c ∈ H (curl; Ω ), q ∈ L2 (Ω ), r ∈ H01 ,

= N1 ,

c1 (d, v, b) = SN2 , ∥d∥curl ∥∇ v∥0 ∥b∥curl

sup d∈X,b,c∈H0 (curl;Ω ),w,u,v∈H10 (Ω )

C (w, d; u, b; v, c) = Nˆ . ∥(w, d)∥∥(u, b)∥∥(v, c)∥

(2.5) (2.6) (2.7) (2.8)

(2.9)

(2.10)

Proof. The estimates of (2.5)–(2.7) can be easily derived by Hölder inequalities. For (2.8), using the Hölder inequality and (2.3), we have

|c0 (w, u, v)| ≤

1 2

1

∥w∥L3 ∥∇ u∥0 ∥v∥L6 + ∥w∥L3 ∥∇ v∥0 ∥u∥L6 ≤ c λ1 ∥∇ w∥0 ∥∇ u∥0 ∥∇ v∥0 , 2

G.-D. Zhang et al. / Computers and Mathematics with Applications 68 (2014) 770–788

773

which implies N1 ≤ c λ1 . For (2.9), using the Hölder inequality, (2.3) and (2.4), we get

|c1 (d, v, b)| ≤ S ∥b∥curl ∥d∥L3 ∥v∥L6 ≤ S λ1 λ2 ∥b∥curl ∥d∥curl ∥∇ v∥0 , √ which yields N2 ≤ λ1 λ2 . (2.10) is followed by (2.8) and (2.9), and Nˆ ≤ 2 max{N1 , N2 }.



3. The well-posedness of MHD equations In this section, we present the well-posedness of MHD equations. Having the above notations, the compact variational form of MHD equations reads as follows: find (u, b) ∈ H10 (Ω ) × H0 (curl; Ω ), (p, r ) ∈ L20 (Ω ) × H01 (Ω ) such that for all (v, c) ∈ H10 (Ω ) × H0 (curl; Ω ), (q, s) ∈ L20 (Ω ) × H01 (Ω ),



A(u, b; v, c) + C (u, b; u, b; v, c) + B(p, r ; v, c) = ⟨L; v, c⟩, B(q, s; u, b) = 0.

(3.1)

Below, we denote

∥f∥−1 = sup v∈H10 (Ω )

⟨f, v⟩ , ∥∇ v∥0

 1 ∥L∥∗ = ∥f∥2−1 + S −1 ∥g∥20 2 .

Lemma 3.1. There exists a β0 > 0 only depending on Ω such that B(q, s; v, c)

sup

∥(v, c)∥

(v,c)∈H10 (Ω )×H0 (curl;Ω )

≥ β0 ∥(q, s)∥,

∀(q, s) ∈ L20 (Ω ) × H01 (Ω ).

(3.2)

Proof. By the fact that ∇ H01 (Ω ) ⊂ H0 (curl; Ω ) and using Part I. Corollary 2.4 in [31], the result can be proved.



In the following, we give the well-posedness of problem (3.1) by a simple method. Our tool is the Banach fixed point theorem. Theorem 3.1. Suppose

Nˆ ∥L∥∗ (min{Re−1 ,Rm−1 λ0 })2

< 1, then, problem (3.1) admits unique solution (u, b) ∈ Υ , (p, r ) ∈ L20 (Ω )× H01 (Ω ).

Moreover, the stability holds

∥(u, b)∥ ≤

∥L∥∗ min{Re−1 , Rm−1 λ0 }

.

(3.3)

Proof. Given (u, b) ∈ Υ , the problem that finding (w, d) ∈ H10 (Ω ) × H0 (curl; Ω ), (p, r ) ∈ L20 (Ω ) × H01 (Ω ) such that for all (v, c) ∈ H10 (Ω ) × H0 (curl; Ω ), (q, s) ∈ L20 (Ω ) × H01 (Ω ) A(w, d; v, c) + C (u, b; w, d; v, c) + B(p, r ; v, c) − B(q, s; w, d) = ⟨L; v, c⟩, has unique solution (w, d) ∈ Υ , (p, r ) ∈ using (2.1), (2.6), we have the estimate:

L20

(Ω ) ×

H01

(3.4)

(Ω ), see Part I. Corollary 4.1 in [31]. Besides, setting v = w, c = d, and

min{Re−1 , Rm−1 λ0 }∥(w, d)∥ ≤ ∥L∥∗ .

(3.5)

Thus, we can define a map G : Υ → Υ, Choosing R =

(u, b) → (w, d).

∥L∥∗ , min{Re−1 ,Rm−1 λ0 }

setting BR = {(v, c) ∈ Υ : ∥(v, c)∥ ≤ R}, evidently, the map G satisfies G : BR → BR .

Let (ui , bi ) ∈ BR and (wi , di ) = G(ui , bi ), i = 1, 2, and (wi , di ), (ui , bi ) satisfy the equations: for all (v, c) ∈ H10 (Ω ) × H0 (curl; Ω ), (q, s) ∈ L20 (Ω ) × H01 (Ω ), A(wi , di ; v, c) + C (ui , bi ; wi , di ; v, c) + B(pi , ri ; v, c) − B(q, s; wi , di ) = ⟨L; v, c⟩,

i = 1, 2.

(3.6)

Subtracting (3.6) as i = 2 from (3.6) as i = 1, and taking (v, c) = (w1 − w2 , d1 − d2 ), (q, s) = (p1 − p2 , r1 − r2 ), using (2.1), we obtain A(w1 − w2 , d1 − d2 ; w1 − w2 , d1 − d2 ) + C (u1 − u2 , b1 − b2 ; w2 , d2 ; w1 − w2 , d1 − d2 ) = 0. Then, by (2.6), (2.10) and (3.5), we have

∥G(u1 , b1 ) − G(u2 , b2 )∥ = ∥(w1 , d1 ) − (w2 , d2 )∥ ≤

Nˆ ∥L∥∗

(min{Re−1 , Rm−1 λ0 })2

∥(u1 , b1 ) − (u2 , b2 )∥.

Thus, G is a contraction mapping on BR → BR . Through the Banach fixed point theorem, G has a fixed point in BR , which is the solution of (3.1).

774

G.-D. Zhang et al. / Computers and Mathematics with Applications 68 (2014) 770–788

Next, we verify the uniqueness of the solutions. Suppose (ui , pi , bi , ri ), i = 1, 2 are two arbitrary solutions of (3.1), i.e., for all (v, c) ∈ H10 (Ω ) × H0 (curl; Ω ), (q, s) ∈ L20 (Ω ) × H01 (Ω ), A(ui , bi ; v, c) + C (ui , bi ; ui , bi ; v, c) + B(pi , ri ; v, c) − B(q, s; ui , bi ) = ⟨L; v, c⟩,

i = 1, 2.

(3.7)

Subtracting (3.7) as i = 2 from (3.7) as i = 1 yields for all (v, c) ∈ H10 (Ω ) × H0 (curl; Ω ), (q, s) ∈ L20 (Ω ) × H01 (Ω ), A(u1 − u2 , b1 − b2 ; v, c) + C (u1 , b1 ; u1 , b1 ; v, c) − C (u2 , b2 ; u2 , b2 ; v, c)

+ B(p1 − p2 , r1 − r2 ; v, c) − B(q, s; u1 − u2 , b1 − b2 ) = 0.

(3.8)

Then, taking (v, c) = (u1 − u2 , b1 − b2 ), (q, s) = (p1 − p2 , r1 − r2 ) in the above equation, and using (2.1), (2.6), (2.10) as well as (3.5), we can obtain min{Re−1 , Rm−1 λ0 }∥(u1 − u2 , b1 − b2 )∥2

≤ A(u1 − u2 , b1 − b2 ; u1 − u2 , b1 − b2 ) = −C (u1 − u2 , b1 − b2 ; u2 , b2 ; u1 − u2 , b1 − b2 ) Nˆ ∥L∥∗ ˆ ∥(u1 − u2 , b1 − b2 )∥2 = ≤ NR ∥(u1 − u2 , b1 − b2 )∥2 , min{Re−1 , Rm−1 λ0 } which combines the condition H10

(Ω ) × H0 (curl; Ω ),

Nˆ ∥L∥∗ (min{Re−1 ,Rm−1 λ0 })2

< 1 yields (u1 , b1 ) = (u2 , b2 ). Next, from (3.8), we have for all (v, c) ∈

B(p1 − p2 , r1 − r2 ; v, c) = 0, which combines Lemma 3.1 gives (p1 , r1 ) = (p2 , r2 ). The proof is completed.



4. Finite element discretization Let {Th } be a family of triangulations (tetrahedrons) of Ω into affine-equivalent finite elements K with Ω = K ∈Th K , which is assumed to be regular in the usual sense. Let hk = diam K , we will assume that the inverse assumption on the mesh  h ≤ C , ∀K ∈ h Th , h = maxK ∈Th hk is satisfied such that we can use inverse inequalities. Pk (K ) denotes the space of h



k

polynomials of total degree at most k ≥ 0 on finite element K and P˜k (K ) denotes the space of homogeneous polynomials of degree k on K . The space Dk (K ) denotes the polynomials p in P˜k (K )d that satisfy p(x) · x = 0 on K . For k ≥ 1, we introduce the space

Nk (K ) = Pk−1 (K )d ⊕ Dk (K ). Note that Nk (K ) ⊂ Pk (K )d . For more details, please see [28–30]. We choose conforming finite element space Vh ⊂ H10 (Ω ), Qh ⊂ L20 (Ω ) to discrete velocity u, pressure p. In addition, Vh , Qh satisfy the uniform LBB condition: for each qh ∈ Qh , there exists a vh ∈ Vh such that bs (qh , vh ) = ∥qh ∥20 ,

∥∇ v∥0 ≤ β1 ∥qh ∥0 ,

(4.1)

here, β1 > 0 only depends on Ω . There are several stable finite element pairs which satisfy this relation, please refer Part II and Part IV in [31]. To approximate the magnetic field b and virtual magnetic pressure r, we use Nédélec’s first family of spaces, see [28–30], given by Ch = {ch ∈ H0 (curl; Ω ) : ch |K ∈ Nk (K ), K ∈ Th }, combined with standard H 1 -conforming space Sh = {sh ∈ H01 (Ω ) : sh |K ∈ Pk (K ), K ∈ Th }. The finite element approximate form of problem (3.1) is finding (uh , bh ) ∈ Vh × Ch , (ph , rh ) ∈ Qh × Sh such that for all (vh , ch ) ∈ Vh × Ch , (qh , sh ) ∈ Qh × Sh  A(uh , bh ; vh , ch ) + C (uh , bh ; uh , bh ; vh , ch ) + B(ph , rh ; vh , ch ) = ⟨L; vh , ch ⟩, (4.2) B(qh , sh ; uh , bh ) = 0. Similar to Lemma 3.1, through (4.1) and ∇ Sh ⊂ Ch , we can verify there exists a β ∗ > 0, only depending on Ω , such that sup

(vh ,ch )∈Vh ×Ch

B(qh , sh ; vh , ch )

∥(vh , ch )∥

≥ β ∗ ∥(qh , sh )∥,

∀(qh , sh ) ∈ Qh × Sh .

Define the discrete kernel space

Υh = {(uh , bh ) ∈ Vh × Ch : B(qh , sh ; uh , bh ) = 0, ∀(qh , sh ) ∈ Qh × Sh } = Zh × Xh ,

(4.3)

G.-D. Zhang et al. / Computers and Mathematics with Applications 68 (2014) 770–788

775

where Zh = {uh ∈ Vh : bs (qh , uh ) = 0, ∀qh ∈ Qh },

Xh = {bh ∈ Ch : bm (sh , bh ) = 0, ∀sh ∈ Sh }.

Noting that Xh ̸⊆ X, the Hodge mapping H : Xh → X is the bridge between Xh and X, which satisfies

∥d − Hd∥0 ≤ λ3 hτ ∥curl d∥0 , ∗

curl(Hd) = curl d,

∀d ∈ X h ,

here τ > , λ3 is a positive constant only depending on Ω , see Lemma 5.4 in [32]. There has the corresponding discrete Poincaré inequality on Xh , see Lemma 7.20 in [28], ∗

1 2

λ4 ∥c∥2curl ≤ ∥curl c∥20 ,

c ∈ Xh ,

where the λ4 is a positive constant only depending on Ω . With this inequality, the bilinear form A(·, ·) is Υh - elliptic, A(uh , bh ; uh , bh ) ≥ min{Re−1 , Rm−1 λ4 }∥(uh , bh )∥2 ,

(uh , bh ) ∈ Υh .

(4.4)

Furthermore, we have some estimates on finite element spaces. Lemma 4.1. There exist positive constants N3 , N4 , N˜ =

√

2 max{N1 , N3 , N4 }, only depending on Ω , such that

|c1 (dh , vh , b)| ≤ SN3 ∥dh ∥curl ∥∇ vh ∥0 ∥b∥curl ,

dh ∈ Xh , vh ∈ Vh , b ∈ H0 (curl; Ω ),

(4.5)

|c1 (dh , v, bh )| ≤ SN4 ∥dh ∥curl ∥∇ v∥0 ∥bh ∥curl ,

dh ∈ X h , v ∈

(4.6)

H10

(Ω ), b ∈ Ch ,

|C (wh , dh ; uh , bh ; vh , ch )| ≤ N˜ ∥(wh , dh )∥∥(uh , bh )∥∥(vh , ch )∥, ∀(wh , dh ) ∈ Υh , uh , vh ∈ Vh , bh , ch ∈ Ch , C (wh , dh ; uh , bh ; uh , bh ) = 0, ∀(wh , dh ) ∈ Υh , uh ∈ Vh , bh ∈ Ch .

(4.7) (4.8)

Proof. For dh ∈ Xh , vh ∈ Vh , b ∈ H0 (curl; Ω ), we have

|c1 (dh , vh , b)| ≤ |c1 (dh − Hdh , vh , b)| + |c1 (Hdh , vh , b)| ≤ S ∥b∥curl ∥dh − Hdh ∥0 ∥vh ∥L∞ + S ∥b∥curl ∥Hdh ∥L3 ∥vh ∥L6 ≤ λ3 Cinv hτ

∗− d

6

S ∥b∥curl ∥dh ∥curl ∥vh ∥L6 + λ1 λ2 S ∥b∥curl ∥Hdh ∥curl ∥∇ vh ∥0 −1

≤ Cinv λ1 λ3 S ∥b∥curl ∥dh ∥curl ∥∇ vh ∥0 + λ0 2 λ1 λ2 S ∥b∥curl ∥∇ vh ∥0 ∥dh ∥curl   − 12 = Cinv λ1 λ3 + λ0 λ1 λ2 S ∥dh ∥curl ∥∇ vh ∥0 ∥b∥curl , −1

we see N3 is well-defined and N3 ≤ Cinv λ1 λ3 + λ0 2 λ1 λ2 . For (4.6), we can verify similarly,

|c1 (dh , v, bh )| ≤ |c1 (dh − Hdh , v, bh )| + |c1 (Hdh , v, bh )| ≤ S ∥curl bh ∥L3 ∥dh − Hdh ∥0 ∥v∥L6 + S ∥curl bh ∥0 ∥Hdh ∥L3 ∥v∥L6 −1

d

′ ≤ Cinv λ1 λ3 hτ − 6 S ∥bh ∥curl ∥∇ v∥0 ∥dh ∥curl + λ0 2 λ1 λ2 S ∥bh ∥curl ∥∇ v∥0 ∥dh ∥curl   −1 ′ ≤ Cinv λ1 λ3 + λ0 2 λ1 λ2 S ∥dh ∥curl ∥∇ v∥0 ∥bh ∥curl . ∗

By (4.5), (4.6) and simple calculation, we can obtain (4.7). (4.8) can be easily derived from the definition of C . Denoting σ =

N˜ ∥L∥∗ , (min{Re−1 ,Rm−1 λ4 })2



we give the well-posedness of problem (4.2).

Theorem 4.1. Suppose σ < 1, then, the problem (4.2) is well-posed. Moreover, the following energy estimate holds,

∥(uh , bh )∥ ≤

∥L∥∗ . min{Re , Rm−1 λ4 }

(4.9)

−1

The proof is similar to Theorem 3.1, here, we omit it. From now on, C denotes a general positive constant, which may depend on Re, Rm, S , f, g and Ω , but not depend on mesh size h, which may stand for different values at different occurrences. Assume the solution of continuous problem (3.1) possesses the following regularity, u ∈ H1+γ (Ω ),

p ∈ H γ (Ω ),

b ∈ Hτ (Ω ),

curl b ∈ Hτ (Ω ),

r ∈ H 1+τ (Ω ),

γ,τ >

1 2

.

Having the regularity, an interpolant operator rh onto the Ch is well-defined and satisfies

∥b − rh b∥curl ≤ Chmin{τ ,k} (∥b∥τ + ∥curl b∥τ ) ,

(4.10)

776

G.-D. Zhang et al. / Computers and Mathematics with Applications 68 (2014) 770–788

see Theorem 5.41 in [28]. We choose the stable velocity and pressure finite element spaces Vh = {vh ∈ H10 (Ω ) : uh |K ∈ Pl (K )d , K ∈ Th }, Qh = {qh ∈ L20 (Ω ) : qh |K ∈ Pl−1 (K ), K ∈ Th },

l > 1.

An example is the ‘‘Taylor–Hood element’’, see Part II Section 4.2 in [31]. The following standard approximation properties hold, inf ∥u − vh ∥1 ≤ Chmin{γ ,l} ∥u∥1+γ ,

(4.11)

inf ∥p − qh ∥0 ≤ Chmin{γ ,l} ∥p∥γ ,

(4.12)

inf ∥r − sh ∥1 ≤ Chmin{τ ,k} ∥r ∥1+τ .

(4.13)

vh ∈Vh

qh ∈Qh

sh ∈Sh

In the below, we state the main result of this section: optimal error estimate. N˜ ∥L∥∗ (min{Re−1 ,Rm−1 λ0 })2

Theorem 4.2. Suppose the conditions of Theorems 3.1, 4.1 and

< 1 hold, let (u, b, p, r ) ∈ Υ × L20 (Ω ) ×

H01 (Ω ) and (uh , bh , ph , rh ) ∈ Υh ×Qh ×Sh be the solutions of the continuous problem (3.1) and discrete problem (4.2), respectively. Then, we have

∥(u − uh , b − bh )∥ ≤ C

inf

(vh ,ch )∈Vh ×Ch  min{γ ,l,τ ,k}

∥u∥1+γ

≤ Ch ∥(p − ph , r − rh )∥ ≤ C

∥(u − vh , b − ch )∥ + C

inf

(vh ,ch )∈Vh ×Ch  min{γ ,l,τ ,k}

≤ Ch

∥(p − qh , r − sh )∥  + ∥b∥τ + ∥curl b∥τ + ∥p∥γ + ∥r ∥1+τ ,

∥(u − vh , b − ch )∥ + C ∥u∥1+γ

inf

(qh ,sh )∈Qh ×Sh

∥(p − qh , r − sh )∥  + ∥b∥τ + ∥curl b∥τ + ∥p∥γ + ∥r ∥1+τ .

(4.14)

inf

(qh ,sh )∈Qh ×Sh

(4.15)

Proof. The error equation is A(u − uh , b − bh ; vh , ch ) + C (u − uh , b − bh ; u, b; vh , ch ) + C (uh , bh ; u − uh , b − bh ; vh , ch )

+ B(p − ph , r − rh ; vh , ch ) = 0,

∀(vh , ch ) ∈ Vh × Ch .

(4.16)

For (wh , dh ) ∈ Υh , (qh , sh ) ∈ Qh × Sh , we have A(uh − wh , bh − dh ; vh , ch ) + C (uh − wh , bh − dh ; u, b; vh , ch ) + C (uh , bh ; uh − wh , bh − dh ; vh , ch )

= A(u − wh , b − dh ; vh , ch ) + C (u − wh , b − dh ; u, b; vh , ch ) + C (uh , bh ; u − wh , b − dh ; vh , ch ) + B(p − qh , r − sh ; vh , ch ) − B(ph − qh , rh − sh ; vh , ch ), ∀(v, c) ∈ Vh × Ch . Taking (vh , ch ) = (uh − wh , bh − dh ) ∈ Υh in the above equation, and using (2.1), we have A(uh − wh , bh − dh ; uh − wh , bh − dh ) + C (uh − wh , bh − dh ; u, b; uh − wh , bh − dh )

= A(u − wh , b − dh ; uh − wh , bh − dh ) + C (u − wh , b − dh ; u, b; uh − wh , bh − dh ) + C (uh , bh ; u − wh , b − dh ; uh − wh , bh − dh ) + B(p − qh , r − sh ; uh − wh , bh − dh ).

(4.17)

Using (4.4), (2.8), (4.5), (4.6) and (3.3), we derive the left hand side of (4.17) has the following estimate

 L.H.S ≥

min{Re

−1

, Rm λ4 } −

N˜ ∥L∥∗

−1



min{Re−1 , Rm−1 λ0 }

where, the coefficient min{Re−1 , Rm−1 λ4 } −

N˜ ∥L∥∗ min{Re−1 ,Rm−1 λ0 }

∥(uh − wh , bh − dh )∥2 ,

(4.18)

> 0 since the assumptions. Before estimating the right hand

side of (4.17), we first bound the second term of the right hand side. Using (2.8), (2.3), the Hölder inequality, H τ (Ω ) ↩→ L3 (Ω ) and H 1+γ (Ω ) ↩→ L∞ (Ω ), we deduce C (u − wh , b − dh ; u, b; uh − wh , bh − dh )

= c0 (u − wh , u, uh − wh ) − c1 (b − dh , uh − wh , b) + c1 (b − dh , u, bh − dh ) ≤ N1 ∥∇(u − wh )∥0 ∥∇ u∥0 ∥∇(uh − wh )∥0 + S ∥curl b∥L3 ∥b − dh ∥0 ∥uh − wh ∥L6 + S ∥bh − dh ∥curl ∥b − dh ∥0 ∥u∥L∞ ≤ N1 ∥∇(u − wh )∥0 ∥∇ u∥0 ∥∇(uh − wh )∥0 + C ∥curl b∥τ ∥b − dh ∥curl ∥∇(uh − wh )∥0 + C ∥bh − dh ∥curl ∥b − dh ∥curl ∥u∥1+γ ≤ C ∥(uh − wh , bh − dh )∥∥(u − wh , b − dh )∥.

(4.19)

G.-D. Zhang et al. / Computers and Mathematics with Applications 68 (2014) 770–788

777

Then, using (2.5), (2.8), (4.5), (4.6), (2.7) and (4.19), we get the right hand side of (4.17) has the following estimate, R.H.S ≤ max{Re−1 , Rm−1 }∥(u − wh , b − dh )∥∥(uh − wh , bh − dh )∥

+ N1 ∥∇ uh ∥0 ∥∇(u − wh )∥0 ∥∇(uh − wh )∥0 + SN3 ∥bh ∥curl ∥∇(uh − wh )∥0 ∥b − dh ∥curl √ + SN4 ∥bh ∥curl ∥∇(u − wh )∥0 ∥bh − dh ∥curl + d∥(p − qh , r − sh )∥∥(uh − wh , bh − dh )∥ + C ∥(uh − wh , bh − dh )∥∥(u − wh , b − dh )∥ √ ≤ C ∥(uh − wh , bh − dh )∥∥(u − wh , b − dh )∥ + d∥(p − qh , r − sh )∥∥(uh − wh , bh − dh )∥.

(4.20)

Combining (4.18) with (4.20) arrives at

∥(uh − wh , bh − dh )∥ ≤ C ∥(u − wh , b − dh )∥ + C ∥(p − qh , r − sh )∥. Then, by the triangle inequality and taking infimum over Υh , Qh , Sh , we get

∥(u − uh , b − bh )∥ ≤ C

inf

(wh ,dh )∈Υh

∥(u − wh , b − dh )∥ + C

inf

(qh ,sh )∈Qh ×Sh

∥(p − qh , r − sh )∥.

(4.21)

Using the LBB condition (4.3) and Lemma 4.1 in Part I of [31], we obtain

√ 

 inf

(wh ,dh )∈Υh

∥(u − wh , b − dh )∥ ≤ 1 +

d

β∗

inf

(vh ,ch )∈Vh ×Ch

∥(u − vh , b − ch )∥.

(4.22)

Combining (4.21) with (4.22) yields

∥(u − uh , b − bh )∥ ≤ C

inf

(vh ,ch )∈Vh ×Ch

∥(u − vh , b − ch )∥ + C

inf

(qh ,sh )∈Qh ×Sh

∥(p − qh , r − sh )∥.

Then, using (4.10)–(4.13), we derive (4.14). We deduce the error estimate of p, r in the below. For any (qh , sh ) ∈ Qh × Sh , (vh , ch ) ∈ Vh × Ch , from (4.16), we have that B(ph − qh , rh − sh ; vh , ch ) = A(u − uh , b − bh ; vh , ch ) + C (u − uh , b − bh ; u, b; vh , ch )

+ C (uh , bh ; u − uh , b − bh ; vh , ch ) + B(p − qh , r − sh ; vh , ch ).

(4.23)

Using (2.8), (2.3), the Hölder inequality, H τ (Ω ) ↩→ L3 (Ω ) and H 1+γ (Ω ) ↩→ L∞ (Ω ), the second term of the right hand side of (4.23) can be bounded C (u − uh , b − bh ; u, b; vh , ch ) ≤ N1 ∥∇(u − uh )∥0 ∥∇ u∥0 ∥∇ vh ∥0 + S ∥curl b∥L3 ∥b − bh ∥0 ∥vh ∥L6

+ S ∥ch ∥curl ∥b − bh ∥0 ∥u∥L∞ ≤ N1 ∥∇(u − uh )∥0 ∥∇ u∥0 ∥∇ vh ∥0 + C ∥curl b∥τ ∥b − bh ∥curl ∥∇ vh ∥0 + C ∥ch ∥curl ∥b − bh ∥curl ∥u∥1+γ ≤ C ∥(u − uh , b − bh )∥∥(vh , ch )∥.

(4.24)

Using (2.5), (2.8), (4.5), (4.6), (2.7) and (4.24), we obtain B(ph − qh , rh − sh ; vh , ch ) ≤ C ∥(u − uh , b − bh )∥∥(vh , ch )∥ +

√

d∥(p − qh , r − sh )∥∥(vh , ch )∥.

Then, by the LBB condition (4.3), the triangle inequality and taking infimum over Qh × Sh , we get

∥(p − ph , r − rh )∥ ≤ C ≤C

inf

∥(p − qh , r − sh )∥ + C ∥(u − uh , b − bh )∥

inf

∥(p − qh , r − sh )∥ + C

(qh ,sh )∈Qh ×Sh (qh ,sh )∈Qh ×Sh

Next, using (4.10)–(4.13), we derive (4.15).

inf

(vh ,ch )∈Vh ×Ch

∥(u − vh , b − ch )∥.



Remark 4.1. In this theorem, b ∈ Hτ (Ω ), curl b ∈ Hτ (Ω ), τ > 12 , which may not belong to H1 (Ω ) and has some singularity. Our results show the formula (4.2) has optimal convergent rates in this general situation. 5. Coupling iterations In this section, we introduce three coupled type iteration algorithms for numerical implementation of the formula (4.2). For each step of the iterations, one need to solve (un , pn , Bn , r n ) simultaneously, a coupled linear system, thus, we call them coupling iterations.

778

G.-D. Zhang et al. / Computers and Mathematics with Applications 68 (2014) 770–788

Iteration I Give (un−1 , bn−1 ) ∈ Υh , solve (un , bn , pn , r n ) ∈ Vh × Ch × Qh × Sh from as (un , vh ) − c1 (bn−1 , vh , bn ) + bs (pn , vh ) − bs (qh , un ) = ⟨f, vh ⟩ − c0 (un−1 , un−1 , vh ), am (b , ch ) + c1 (b

n −1

n

, u , ch ) + bm (r , ch ) − bm (sh , b ) = (g, ch ), n

n

n

∀(vh , qh ) ∈ Vh × Qh ,

∀(ch , sh ) ∈ Ch × Sh .

(5.1) (5.2)

Iteration II Give (un−1 , bn−1 ) ∈ Υh , solve (un , bn , pn , r n ) ∈ Vh × Ch × Qh × Sh from as (un , vh ) + c0 (un−1 , un , vh ) + c0 (un , un−1 , vh ) − c1 (bn−1 , vh , bn ) − c1 (bn , vh , bn−1 ) + bs (pn , vh )

− bs (qh , un ) = ⟨f, vh ⟩ + c0 (un−1 , un−1 , vh ) − c1 (bn−1 , vh , bn−1 ), ∀(vh , qh ) ∈ Vh × Qh , am (bn , ch ) + c1 (bn−1 , un , ch ) + c1 (bn , un−1 , ch ) + bm (r n , ch ) − bm (sh , bn ) = (g, ch ) + c1 (bn−1 , un−1 , ch ), ∀(ch , qh ) ∈ Ch × Qh . Iteration III Give (un−1 , bn−1 ) ∈ Υh , solve (un , bn , pn , r n ) ∈ Vh × Ch × Qh × Sh from as (un , vh ) + c0 (un−1 , un , vh ) − c1 (bn−1 , vh , bn ) + bs (pn , vh ) − bs (qh , un ) = ⟨f, vh ⟩, am (bn , ch ) + c1 (bn−1 , un , ch ) + bm (r n , ch ) − bm (sh , bn ) = (g, ch ),

∀(vh , qh ) ∈ Vh × Qh , ∀(ch , qh ) ∈ Ch × Qh .

The initial value (u0 , b0 ) ∈ Υh is obtained from as (u0 , vh ) + bs (p0 , vh ) − bs (qh , u0 ) = ⟨f, vh ⟩,

∀(vh , qh ) ∈ Vh × Qh ,

(5.3)

∀(ch , sh ) ∈ Ch × Sh .

(5.4)

am (b , ch ) + bm (r , ch ) − bm (sh , b ) = (g, ch ), 0

0

0

It is easy to see (u , b ) is well defined and 0

∥(u0 , b0 )∥ ≤

0

∥L∥∗ , min{Re , Rm−1 λ4 } −1

∥(uh − u0 , bh − b0 )∥ ≤ σ

∥L∥∗ . min{Re , Rm−1 λ4 }

(5.5)

−1

In the rest of this section, we analyze the stabilities and convergence of the three iterations. Theorem 5.1. Assume 0 < σ ≤

∥(un , bn )∥ ≤

1 , 4

then Iteration I is stable and convergent, namely, for n ≥ 0,

2∥L∥∗ min{Re−1 , Rm−1 λ4 }

β ∗ ∥(pn , r n )∥ ≤ 2∥L∥∗ +

,

(5.6)

2 max{Re−1 , Rm−1 }∥L∥∗ min{Re−1 , Rm−1 λ4 }

∥(uh − un , bh − bn )∥ ≤ (3σ )n

,

∥L∥∗ 4 min{Re−1 , Rm−1 λ4 }

β ∗ ∥(ph − pn , rh − r n )∥ ≤ (3σ )n



(5.7)

,

max{Re−1 , Rm−1 }∥L∥∗ 4 min{Re−1 , Rm−1 λ4 }

(5.8)

+

5∥L∥∗ 16



.

Proof. We first prove (5.6) by the induction method. Observing that ∥(u0 , b0 )∥ ≤ bn−1 )∥ ≤

2∥L∥∗ min{Re−1 ,Rm−1 λ4 }

and prove ∥(un , bn )∥ ≤

2∥L∥∗ . min{Re−1 ,Rm−1 λ4 }

(5.9) 2∥L∥∗ , min{Re−1 ,Rm−1 λ4 } pn in (5.1), ch

Taking vh = un , qh =

we suppose ∥(un−1 ,

= bn , sh = r n in (5.2)

and adding the two equations, using (2.8), (4.4), we have

min{Re−1 , Rm−1 λ4 }∥(un , bn )∥2 ≤ A(un , bn ; un , bn ) = ⟨L; un , bn ⟩ − c0 (un−1 , un−1 , un )

≤ ∥L∥∗ ∥(un , bn )∥ + N1 ∥∇ un−1 ∥20 ∥∇ un ∥0 . Then using the induction assumption and σ ≤

1 , 4

we derive

N˜ ∥L∥∗ + ∥(un−1 , bn−1 )∥2 − 1 min{Re , Rm−1 λ4 } min{Re , Rm−1 λ4 } N˜ 4∥L∥2∗ ∥L∥∗ ≤ + min{Re−1 , Rm−1 λ4 } min{Re−1 , Rm−1 λ4 } min{Re−1 , Rm−1 λ4 }2 2∥L∥∗ ≤ . min{Re−1 , Rm−1 λ4 } Adding (5.1) and (5.2) and taking qh = sh = 0 gives

∥(un , bn )∥ ≤

−1

B(pn , r n ; vh , ch ) = ⟨L; vh , ch ⟩ − A(un , bn ; vh , ch ) − c0 (un−1 , un−1 , vh ) + c1 (bn−1 , vh , bn ) + c1 (bn−1 , un , ch ), Then, using (2.5), (2.8), (4.5), (5.6) and the discrete version LBB condition (4.3), we derive (5.7).

G.-D. Zhang et al. / Computers and Mathematics with Applications 68 (2014) 770–788

779

Problem (4.2) subtracting (5.1), (5.2) gives A(uh − un , bh − bn ; vh , ch ) + c1 (bn−1 − bh , vh , bn ) + c1 (bh , vh , bn − bh ) + c1 (bh , uh − un , ch )

+ c1 (bh − bn−1 , un , ch ) + c0 (uh , uh − un−1 , vh ) + c0 (uh − un−1 , un−1 , vh ) + B(ph − pn , rh − r n ; vh , ch ) − B(qh , sh ; uh − un , bh − bn ) = 0,

∀(vh , ch ) ∈ Vh × Ch , (qh , sh ) ∈ Qh × Sh .

(5.10)

Taking vh = uh − u , ch = bh − b , qh = ph − p , sh = rh − r , we obtain n

n

n

n

A(uh − un , bh − bn ; uh − un , bh − bn ) + c1 (bn−1 − bh , uh − un , bn ) + c1 (bh − bn−1 , un , bh − bn )

+ c0 (uh , uh − un−1 , uh − un ) + c0 (uh − un−1 , un−1 , uh − un ) = 0. Using (4.4), (4.5), (2.8), (4.9) and (5.6), we obtain

∥(uh − un , bh − bn )∥ ≤ 3σ ∥(uh − un−1 , bh − bn−1 )∥ ≤ (3σ )n ∥(uh − u0 , bh − b0 )∥, which combines (5.5) and yield (5.8). Next, taking qh = 0, sh = 0 in (5.10), using (2.5), (4.5), (2.8), (4.9), (5.6) and discrete version LBB condition (4.3) give that

 β ∥(ph − p , rh − r )∥ ≤ ∗

n

−1

n

max{Re

+

, Rm } +

3N˜ ∥L∥∗ min{Re−1 , Rm−1 λ4 }

which combines (5.8) to arrive at (5.9).



N˜ ∥L∥∗

−1

min{Re−1 , Rm−1 λ4 }

∥(uh − un , bh − bn )∥

∥(uh − un−1 , bh − bn−1 )∥,



In order to analyze Iteration II, it is convenient to write Iteration II in the following compact form. A(un , bn ; vh , ch ) + C (un−1 , bn−1 ; un , bn ; vh , ch ) + C (un , bn ; un−1 , bn−1 ; vh , ch )

− C (un−1 , bn−1 ; un−1 , bn−1 ; vh , ch ) + B(pn , r n ; v, c) − B(qh , sh ; un , bn ) = ⟨L; vh , ch ⟩, ∀(vh , ch ) ∈ Vh × Ch , (qh , sh ) ∈ Qh × Sh , or, A(un , bn ; vh , ch ) + C (un−1 − un , bn−1 − bn ; un − un−1 , bn − bn−1 ; vh , ch ) + C (un , bn ; un , bn ; vh , ch ) + B(pn , r n ; vh , ch ) − B(qh , sh ; un , bn ) = ⟨L; vh , ch ⟩, ∀(vh , ch ) ∈ Vh × Ch , (qh , sh ) ∈ Qh × Sh .

(5.11)

(5.12)

Now, we give the stability and convergence results of Iteration II. Theorem 5.2. Assume 0 < σ ≤

∥(un , bn )∥ ≤

5 , 11

then, Iteration II is stable and convergent, namely, for n ≥ 0,

4∥L∥∗ 3 min{Re−1 , Rm−1 λ4 }

β ∗ ∥(pn , r n )∥ ≤

11 3

∥L∥∗ +

∥(uh − un , bh − bn )∥ ≤

,

(5.13)

4 max{Re−1 , Rm−1 }∥L∥∗ 3 min{Re−1 , Rm−1 λ4 }



5∥L∥∗ 11 min{Re−1 , Rm−1 λ4 }

β ∗ ∥(ph − pn , rh − r n )∥ ≤



,

15 13

5 max{Re−1 , Rm−1 }∥L∥∗ 11 min{Re−1 , Rm−1 λ4 }

(5.14)

σ

2n −1

3

,

+ ∥L∥∗ 4

(5.15)



15 13

σ

2n −1

.

Proof. We first prove (5.13) by the induction method. It is easy to see that ∥(u0 , b0 )∥ ≤ and (u1 − u0 , b1 − b0 ) satisfies

(5.16) ∥L∥∗ min{Re−1 ,Rm−1 λ4 }

≤

4∥L∥∗ 3 min{Re−1 ,Rm−1 λ4 }

A(u1 − u0 , b1 − b0 ; u1 − u0 , b1 − b0 ) + C (u1 − u0 , b1 − b0 ; u0 , b0 ; u1 − u0 , b1 − b0 )

+ C (u0 , b0 ; u0 , b0 ; u1 − u0 , b1 − b0 ) = 0. Through the above equation and (4.4), (4.7), (5.5), we have

∥(u1 − u0 , b1 − b0 )∥ ≤

5∥L∥∗ 6 min{Re−1 , Rm−1 λ4 }

.

Setting n = 1, (vh , ch ) = (u1 , b1 ) and (qh , sh ) = (p1 , r 1 ) in (5.12), we have A(u1 , b1 ; u1 , b1 ) + C (u0 − u1 , b0 − b1 ; u1 − u0 , b1 − b0 ; u1 , b1 ) = ⟨L; u1 , b1 ⟩.

(5.17)

780

G.-D. Zhang et al. / Computers and Mathematics with Applications 68 (2014) 770–788

Then, using (4.4), (4.7) and (5.17), we obtain

∥(u1 , b1 )∥ ≤

4∥L∥∗ 3 min{Re−1 , Rm−1 λ4 }

Now, we suppose ∥(uJ , bJ )∥ ≤

.

4∥L∥∗ 3 min{Re−1 ,Rm−1 λ4 }

, J = 0, 1, . . . , n − 1, n ≥ 2 and prove ∥(un , bn )∥ ≤

n ≥ 2, the following equation holds

4∥L∥∗ . 3 min{Re−1 ,Rm−1 λ4 }

For

A(un − un−1 , bn − bn−1 ; un − un−1 , bn − bn−1 ) + C (un − un−1 , bn − bn−1 ; un−1 , bn−1 ; un − un−1 , bn − bn−1 )

+ C (un−1 − un−2 , bn−1 − bn−2 ; un−1 − un−2 , bn−1 − bn−2 ; un − un−1 , bn − bn−1 ) = 0. Using (4.4), (4.7), the induction assumption and σ ≤

5 , 11

we obtain

min{Re−1 , Rm−1 λ4 }∥(un − un−1 , bn − bn−1 )∥

≤ N˜ ∥(un−1 , bn−1 )∥∥(un − un−1 , bn − bn−1 )∥ + N˜ ∥(un−1 − un−2 , bn−1 − bn−2 )∥2 4N˜ ∥L∥∗ ≤ ∥(un − un−1 , bn − bn−1 )∥ + N˜ ∥(un−1 − un−2 , bn−1 − bn−2 )∥2 3 min{Re−1 , Rm−1 λ4 } ≤

20 33

min{Re−1 , Rm−1 λ4 }∥(un − un−1 , bn − bn−1 )∥ + N˜ ∥(un−1 − un−2 , bn−1 − bn−2 )∥2 ,

which yields

∥(un − un−1 , bn − bn−1 )∥ ≤

33N˜ 13 min{Re−1 , Rm−1 λ4 }

Thus, by (5.17) and our assumption σ ≤

5 , 11

we derive that

 n

n−1

∥(u − u

,b − b n

n−1

)∥ ≤

2n−1 −1

33N˜

∥(u1 − u0 , b1 − b0 )∥2

13 min{Re−1 , Rm−1 λ4 }

 ≤ ≤

∥(un−1 − un−2 , bn−1 − bn−2 )∥2 .

11 5

σ

2n−1 −1

n−1

5∥L∥∗ 6 min{Re−1 , Rm−1 λ4 }

5∥L∥∗ 6 min{Re−1 , Rm−1 λ4 }

.

(5.18)

Then, setting (vh , ch ) = (un , bn ), (qh , sh ) = (pn , r n ) in (5.12) yields A(un , bn ; un , bn ) + C (un−1 − un , bn−1 − bn ; un − un−1 , bn − bn−1 ; un , bn ) = ⟨L, (un , bn )⟩. Using (4.4), (4.7), (5.18) and σ ≤

5 , 11

we arrive at

∥L∥∗ 25N˜ ∥L∥2∗ + − 1 min{Re , Rm λ4 } 36(min{Re−1 , Rm−1 λ4 })3 4∥L∥∗ ≤ . 3 min{Re−1 , Rm−1 λ4 }

∥(un , bn )∥ ≤

−1

Next, by (5.11), (2.5), (4.7), (5.13) and the discrete version LBB condition (4.3), we derive (5.14). Eq. (4.2) subtracting (5.11) and setting (vh , ch ) = (uh − un , bh − bn ), (qh , sh ) = (ph − pn , rh − r n ) gives A(uh − un , bh − bn ; uh − un , bh − bn ) + C (uh − un−1 , bh − bn−1 ; uh − un−1 , bh − bn−1 ; uh − un , bh − bn )

+ C (uh − un , bh − bn ; un−1 , bn−1 ; uh − un , bh − bn ) = 0,

n ≥ 1,

which combines (4.4), (4.7) and (5.13) derives

∥(uh − un , bh − bn )∥ ≤

33N˜ 13 min{Re−1 , Rm−1 λ4 }

 ≤

∥(uh − un−1 , bh − bn−1 )∥2

33N˜ 13 min{Re−1 , Rm−1 λ4 }

Then, using (5.5), we obtain (5.15).

2n −1 n

∥(uh − u0 , bh − b0 )∥2 ,

n ≥ 1.

G.-D. Zhang et al. / Computers and Mathematics with Applications 68 (2014) 770–788

781

We find (ph − pn , rh − r n ) satisfies the following equation from (4.2) subtracting (5.11) B(pn − ph , r n − rh ; v, c) = A(uh − un , bh − bn ; v, c) + C (uh − un−1 , bh − bn−1 ; uh − un−1 , bh − bn−1 ; v, c)

+ C (uh − un , bh − bn ; un−1 , bn−1 ; v, c) + C (un−1 , bn−1 ; uh − un , bh − bn ; v, c). Using (2.5), (4.3) and (4.7), we have

β ∗ ∥(pn − ph , r n − rh )∥ ≤ max{Re−1 , Rm−1 }∥(uh − un , bh − bn )∥ + N˜ ∥(uh − un−1 , bh − bn−1 )∥2 + 2N˜ ∥(un−1 , bn−1 )∥∥(uh − un , bh − bn )∥. Then, using (5.13) and (5.15), we obtain (5.16).



For analyzing Iteration III, we also write its compact form: A(un , bn ; vh , ch ) + C (un−1 , bn−1 ; un , bn ; vh , ch ) + B(pn , r n ; vh , ch ) − B(qh , sh ; un , bn )

= ⟨L; vh , ch ⟩,

∀(vh , ch ) ∈ Vh × Ch , (qh , sh ) ∈ Qh × Sh .

(5.19)

We state the stability and convergence results of Iteration III in the below. Theorem 5.3. Iteration III is unconditionally stable, namely, for n ≥ 0,

∥(un , bn )∥ ≤

∥L∥∗ min{Re−1 , Rm−1 λ4 }

β ∗ ∥(pn , r n )∥ ≤ ∥L∥∗ +

,

(5.20)

max{Re−1 , Rm−1 }∥L∥∗ −1

min{Re

, Rm−1 λ

4

+

}

N˜ ∥L∥2∗

(min{Re−1 , Rm−1 λ4 })2

.

(5.21)

Furthermore, assume 0 < σ < 1, Iteration III is convergent, namely,

∥(uh − un , bh − bn )∥ ≤ σ n

∥L∥∗ min{Re−1 , Rm−1 λ4 }

β ∗ ∥(ph − pn , rh − r n )∥ ≤ σ n



,

(5.22)

max{Re−1 , Rm−1 }∥L∥∗ min{Re−1 , Rm−1 λ4 }

 + 2∥L∥∗ .

(5.23)

Proof. Choosing (vh , ch ) = (un , bn ), (qh , sh ) = (pn , r n ) in (5.19) and using (4.4), (4.8), we derive (5.20). Next, taking qh = 0, sh = 0 in (5.19), using (2.5), (4.7), (4.3) and (5.20), we get (5.21). (4.2) subtracting (5.19) gives the error equation, for n ≥ 1, A(uh − un , bh − bn ; vh , ch ) + C (uh − un−1 , bh − bn−1 ; uh , bh ; vh , ch )

+ C (un−1 , bn−1 ; uh − un , bh − bn ; vh , ch ) + B(ph − pn , rh − r n ; vh , ch ) − B(qh , sh ; uh − un , bh − bn ) = 0,

∀(vh , ch ) ∈ Vh × Ch , (q, s) ∈ Qh × Sh .

(5.24)

Choosing (vh , ch ) = (uh − u , bh − b ), (qh , sh ) = (ph − p , rh − r ) in (5.24) and using (4.4), (4.7), (4.8) and (4.9), we get for all n ≥ 1, n

n

n

n

∥(uh − un , bh − bn )∥ ≤ σ ∥(uh − un−1 , bh − bn−1 )∥ ≤ σ n ∥(uh − u0 , bh − b0 )∥, which combines (5.5) yield (5.22). Next, choosing qh = 0, sh = 0 in (5.24) gives B(pn − ph , r n − rh ; vh , ch ) = A(uh − un , bh − bn ; vh , ch ) + C (uh − un−1 , bh − bn−1 ; uh , bh ; vh , ch )

+ C (un−1 , bn−1 ; uh − un , bh − bn ; vh , ch ),

∀(vh , ch ) ∈ Vh × Ch .

Then, using (2.5), (4.7), (5.20), (5.22), (4.9) and LBB condition (4.3), we obtain (5.23).



Remark 5.1. In Theorems 5.1, 5.2, we give the sufficient conditions of stability and convergence for Iterations I and II. But, it seems impossible to determine their accurate critical conditions. Anyway, through our above numerical analysis and the later numerical tests, we could draw the conclusion: the condition of stability and convergence of Iteration II is weaker than Iteration I, and the condition of stability and convergence of Iteration III is weaker than Iteration II. Finally, combining Theorem 4.2 with Theorems 5.1–5.3, we obtain the error estimates of Iteration I–III with respect to grid size h and iterative step n. Theorem 5.4. Suppose the conditions of Theorem 4.2 are valid. (i) Assume 0 < σ ≤

1 , 4

Iteration I is stable and convergent, and there holds:

∥(u − un , b − bn )∥ + ∥(p − pn , r − r n )∥ ≤ Chmin{γ ,l,τ ,k} + C (3σ )n .

782

G.-D. Zhang et al. / Computers and Mathematics with Applications 68 (2014) 770–788

(ii) Assume 0 < σ ≤

5 , 11

Iteration II is stable and convergent, and there holds: min{γ ,l,τ ,k}

∥(u − u , b − b )∥ + ∥(p − p , r − r )∥ ≤ Ch n

n

n

n

 +C

15 13

σ

2n −1

.

(iii) Iteration III is unconditional stable, and assume 0 < σ < 1, it is convergent, there holds

∥(u − un , b − bn )∥ + ∥(p − ph , r − rh )∥ ≤ Chmin{γ ,l,τ ,k} + C σ n . 6. Numerical experiments In this section, we present three numerical examples to check the theoretical results developed in the previous sections. The first two verify the convergent rates and compare the CPU cost of the three iterative methods. The last one is a classical MHD model problem, Hartmann flow, which we use to check stability and convergence of the three iterations under different physics parameters. We choose Mini element (P1b , P1 ) for hydrodynamic unknowns (u, p), see Part II Section 4.1 in [31], the lowest order Nédélec’s element for magnetic field b, P1 element for virtual magnet pressure r. Example 1 (MHD Problem with a Smooth Solution). In this example, we consider the MHD equations on (0, 1) × (0, 1) ∈ R2 , the physics parameters Re = Rm = S = 1. Choose source terms f, g and boundary conditions such that the MHD equations has the solution: u1 = x2 (x − 1)2 y(y − 1)(2y − 1), u2 = −y2 (y − 1)2 x(x − 1)(2x − 1), p = (2x − 1)(2y − 1), b1 = sin(π x) cos(π y), b2 = − sin(π y) cos(π x), r = 0. We first solve two small linear equations (5.3) and (5.4) to get the initial value, then solve the MHD equations using the three iterations until the L2 -norm of the difference in successive iteration less than the tolerance Exp = 1.0e−7. The numerical results are shown in Table 1. Seen from Table 1, the formula provides optimal error rates using the three iterations, i.e., ∥∇(u − uh )∥0 and ∥b − bh ∥curl ∼ O(h), the order of ∥p − ph ∥0 approaches O(h) and the Lagrange multiplier rh goes to zero as mesh size h refines. The numerical accuracy of three iterations is almost the same, and Iteration II costs the least CPU time because it has two-order convergent rate about the iterative step n. Example 2 (MHD Problem with a Singular Solution). In this example, we consider the MHD equations on a non-convex L-shaped domain Ω = (−1, 1) × (−1, 1) \ [0, 1] × [0, −1]. We set Re = Rm = 0.1, S = 1, and choose the force terms and the boundary conditions such that the analytical solutions have strong singularities at the re-entrant corner. The singular solutions are described in polar coordinates (ρ, θ ), see [15,21]. u1 (x, y) = ρ λ (1 + λ) sin(θ )ψ(θ ) + cos(θ )ψ ′ (θ ) ,





u2 (x, y) = ρ λ −(1 + λ) cos(θ )ψ(θ ) + sin(θ )ψ ′ (θ ) ,





 ρ λ−1  (1 + λ)2 ψ ′ (θ ) + ψ ′′′ (θ ) , 1−λ    2 2 θ , (b1 (x, y), b2 (x, y)) = ∇ ρ 3 sin p(x, y) = −

3

r =0 where the parameter λ ≈ 0.54448,

ψ(θ ) = sin((1 + λ)θ )

cos(λω) 1+λ

− cos((1 + λ)θ ) − sin((1 − λ)θ )

cos(λω)

2

1−λ

+ cos((1 − λ)θ ).

Note that (u, p) ∈ H1+λ (Ω ) × H λ (Ω ), ∇ · b = 0, curl b = 0 and b ∈ H 3 (Ω ). We compute this problem using the same procedure in example 1, here the tolerance Exp = 1.0e−6. The data are demonstrated in Table 2, from which we see the error rates have optimal convergent order, except the rate of ∥b − bh ∥curl 1 suffers a little loss when h = 64 , and the Lagrange multiplier rh approaches zero as the h refines. Besides, the numerical accuracy of the three iterations is almost the same, and Iteration II consumes the least CPU time.

G.-D. Zhang et al. / Computers and Mathematics with Applications 68 (2014) 770–788

783

Table 1 Three iterations for the MHD problem with smooth solution (Exp = 1.0e−7). h

Type

CPU (s)

∥∇(u − uh )∥0

1/4 1/4 1/4 1/8 1/8 1/8 1/16 1/16 1/16 1/32 1/32 1/32 1/64 1/64 1/64

Iteration I Iteration II Iteration III Iteration I Iteration II Iteration III Iteration I Iteration II Iteration III Iteration I Iteration II Iteration III Iteration I Iteration II Iteration III

0.187 0.219 0.219 0.641 0.61 0.719 3.062 2.562 2.829 11.266 10.235 13.047 46.422 43.188 50.672

0.0228487 0.0228487 0.0228487 0.0100792 0.0100792 0.0100792 0.00481518 0.00481518 0.00481518 0.00236519 0.00236519 0.00236519 0.00117441 0.00117441 0.00117441

Rate

∥p − ph ∥0

1.18 1.18 1.18 1.06 1.06 1.06 1.02 1.02 1.02 1.01 1.01 1.01

0.0367408 0.0367408 0.0367408 0.00976416 0.00976416 0.00976416 0.00271944 0.00271944 0.00271944 0.000807111 0.000807111 0.000807111 0.000255667 0.000255667 0.000255667

Rate

∥b − bh ∥curl

Rate

∥r − rh ∥1

1.91 1.91 1.91 1.84 1.84 1.84 1.75 1.75 1.75 1.65 1.65 1.65

0.825439 0.825439 0.825439 0.417375 0.417375 0.417375 0.209286 0.209286 0.209286 0.104718 0.104718 0.104718 0.0523685 0.0523685 0.0523685

0.98 0.98 0.98 0.99 0.99 0.99 0.99 0.99 0.99 0.99 0.99 0.99

1.23249e−7 1.23249e−7 1.23249e−7 5.67632e−10 5.67631e−10 5.67632e−10 2.33957e−12 2.33757e−12 2.33957e−12 9.39161e−15 2.23926e−13 9.3889e−15 1.34789e−15 8.86329e−13 1.36139e−15

Rate

∥b − bh ∥curl

Rate

∥r − rh ∥0

1.02 1.02 1.03 0.84 0.84 0.84 0.70 0.70 0.71 0.57 0.57 0.56

0.181438 0.181438 0.18144 0.114114 0.114113 0.114113 0.0728583 0.0728582 0.0728579 0.0459007 0.0459007 0.0460386 0.0297594 0.0297594 0.0297594

0.67 0.67 0.67 0.65 0.65 0.65 0.67 0.67 0.66 0.63 0.63 0.63

0.00925069 0.00925069 0.00924446 0.00426083 0.00426083 0.0042389 0.000587599 0.000587599 0.000589569 0.000251998 0.000251998 0.000274508 0.000153217 0.000153217 0.000153217

Table 2 Three iterations for the MHD problem with singular solution (Exp = 1.0e−6). h

Type

CPU (s)

∥u − uh ∥1

1/4 1/4 1/4 1/8 1/8 1/8 1/16 1/16 1/16 1/32 1/32 1/32 1/64 1/64 1/64

Iteration I Iteration II Iteration III Iteration I Iteration II Iteration III Iteration I Iteration II Iteration III Iteration I Iteration II Iteration III Iteration I Iteration II Iteration III

1.925 1.096 1.966 4.182 2.612 4.446 13.151 8.499 13.821 48.242 32.827 52.868 200.719 132.447 213.877

1.86798 1.86801 1.86795 1.22087 1.22092 1.22075 0.762744 0.762744 0.762724 0.519587 0.519587 0.519046 0.357444 0.357444 0.357444

Rate

∥p − ph ∥0

0.61 0.61 0.61 0.67 0.67 0.67 0.55 0.55 0.55 0.54 0.54 0.54

24.679 24.679 24.7601 12.1275 12.1275 12.1065 6.75149 6.75149 6.75029 4.15259 4.15259 4.12037 2.79142 2.79142 2.79142

Fig. 1. Hartmann flow problem: model setting.

Example 3 (Hartmann Flow). Hartmann flow is a classical MHD problem, which describes the flow of a liquid metal through a channel, see [20]. The liquid flows in the x-direction under the influence of a prescribed pressure gradient G, meanwhile, a uniform magnetic field b0 = (0, 1) is applied on the boundaries. The problem model is described in Fig. 1. The Hartmann

784

G.-D. Zhang et al. / Computers and Mathematics with Applications 68 (2014) 770–788

Fig. 2. Hartmann flow in 2D: numerical (points) and analytical (line) u1 along x = 0.5.

Fig. 3. Hartmann flow in 2D: numerical (points) and analytical (line) b1 along x = 0.5.

√ number Ha = u1 = b1 =

ReRmS, the exact solution is given by Re G



Ha tanh(Ha) G S



sinh(Ha y) sinh(Ha)

1−

cosh(Ha y) cosh(Ha)



−y ,



,

u2 = 0,

b2 = 1.

In this test, we set pressure gradient G = 1, and choose the following three cases to compute this problem. Case a : Re = 1, Rm = 0.01, S = 100, Ha = 1. Case b : Re = 50, Rm = 0.02, S = 400, Ha = 20. Case c : Re = 500, Rm = 0.02, S = 1000, Ha = 90. The numerical result of Case a is shown in Table 3. We observe that the three iterations can compute in this case, and the rates of ∥u − uh ∥0 , ∥u − uh ∥1 and ∥b − bh ∥curl all have optimal order. The three iterations have almost the same numerical

G.-D. Zhang et al. / Computers and Mathematics with Applications 68 (2014) 770–788

785

Table 3 Hartmann flow with Re = 1, Rm = 0.01, S = 100, Ha = 1. h

Type

CPU (s)

∥u − uh ∥0

1/2 1/2 1/2 1/4 1/4 1/4 1/8 1/8 1/8 1/16 1/16 1/16 1/32 1/32 1/32

Iteration I Iteration II Iteration III Iteration I Iteration II Iteration III Iteration I Iteration II Iteration III Iteration I Iteration II Iteration III Iteration I Iteration II Iteration III

0.437 0.375 0.343 1.672 1.625 1.891 8.031 7.687 7.75 32.781 33.562 33.061 165.547 128.328 177.125

0.0649506 0.0649506 0.0649495 0.0162761 0.0162761 0.0162761 0.00406323 0.00406323 0.00406323 0.00101215 0.00101215 0.00101215 0.000252866 0.000252866 0.000252866

Rate

∥u − uh ∥1

Rate

∥b − bh ∥curl

Rate

∥r − rh ∥0

1.99 1.99 1.99 2.00 2.00 2.00 2.00 2.00 2.00 2.00 2.00 2.00

0.398687 0.398687 0.398703 0.191058 0.191058 0.191059 0.0931592 0.0931592 0.0931593 0.0461522 0.0461522 0.0461522 0.0230238 0.0230238 0.0230238

1.06 1.06 1.06 1.03 1.03 1.03 1.01 1.01 1.01 1.00 1.00 1.00

0.00173651 0.00173651 0.00173651 0.000901453 0.000901453 0.000901453 0.000446211 0.000446211 0.000446211 0.000221025 0.000221025 0.000221025 0.000109941 0.000109941 0.000109941

0.94 0.94 0.94 1.01 1.01 1.01 1.01 1.01 1.01 1.01 1.01 1.01

3.58927e−12 6.56929e−12 3.58987e−12 1.28936e−11 2.11202e−11 1.28933e−11 3.35182e−11 6.03606e−11 3.35044e−11 9.74578e−11 1.4349e−10 9.74526e−11 2.02328e−10 4.39295e−10 2.02302e−10

∥u − uh ∥1

Rate

∥b − bh ∥curl

Rate

∥r − rh ∥0

Table 4 Hartmann flow with Re = 50, Rm = 0.02, S = 400, Ha = 20. h

Type

1/4 1/4 1/4 1/8 1/8 1/8 1/16 1/16 1/16 1/32 1/32 1/32

Iteration I Iteration II Iteration III Iteration I Iteration II Iteration III Iteration I Iteration II Iteration III Iteration I Iteration II Iteration III

CPU (s) 2.406 6.922 10.859 24.11 47.063 66.906 227.547 279.875

∥u − uh ∥0 N/A 0.904382 0.910733 N/A 0.352477 0.354827 N/A 0.0963823 0.0965139 N/A 0.0258956 0.0259107

Rate

1.36 1.36 1.87 1.88 1.90 1.90

N/A 17.1912 17.794 N/A 11.9432 12.148 N/A 6.33279 6.35607 N/A 3.32836 3.33299

0.52 0.55 0.91 0.93 0.93 0.93

N/A 0.016825 0.0170087 N/A 0.0110284 0.0110628 N/A 0.00571459 0.00571535 N/A 0.00294912 0.00294916

0.95 0.95

N/A 4.73402e−11 2.57895e−11 N/A 1.12455e−10 6.70286e−11 N/A 2.68638e−10 1.94896e−10 N/A 7.59839e−10 4.04619e−10

Rate

∥r − rh ∥0

0.61 0.62 0.95 0.95

Table 5 Hartmann flow with Re = 400, Rm = 0.02, S = 1000, Ha = 90. h

Type

1/16 1/16 1/16 1/32 1/32 1/32 1/48 1/48 1/48 1/60 1/60 1/60

Iteration I Iteration II Iteration III Iteration I Iteration II Iteration III Iteration I Iteration II Iteration III Iteration I Iteration II Iteration III

CPU (s)

45.609

150.172

292.453

417.907

∥u − uh ∥0 N/A N/A 0.413541 N/A N/A 0.153842 N/A N/A 0.0713197 N/A N/A 0.0435406

Rate

1.43

1.90

2.21

∥u − uh ∥1 N/A N/A 39.169 N/A N/A 24.5723 N/A N/A 15.9576 N/A N/A 12.1468

Rate

0.67

1.06

1.22

∥b − bh ∥curl N/A N/A 0.00853968 N/A N/A 0.00462612 N/A N/A 0.00304033 N/A N/A 0.00239276

0.88

1.03

1.07

N/A N/A 2.08289e−10 N/A N/A 4.04476e−10 N/A N/A 6.10873e−10 N/A N/A 8.16397e−10

accuracy, meantime, Iteration II has the best performance in computational time on fine mesh. Figs. 2 and 3 demonstrate the numerical solutions of this case computed on 4540 elements. We see the numerical solutions of the three iterations coincide with exact solution. Table 4 is the data of Case b, from which we observe Iteration I cannot run in this case. The convergence rates of Iterations II and III approach optimal order as mesh size h refines, and the Iteration II still costs less CPU time than Iteration III. Figs. 4, 5 display the numerical solutions of this case computed on 18 582 elements, from which we find the numerical solutions of Iterations II and III coincide with exact solutions. Finally, we observe the data of Case c, see Table 5. In this case, Iterations I and II are all divergent, only Iteration III can work in the high Hartmann number case, and the convergence rates of Iteration III have optimal order. Figs. 6, 7 are the numerical results computed on 65 926 elements, from which we find the numerical solutions of Iteration III almost agree with the analytical solution. The Lagrange multiplier rh all almost equals zero in the above three cases. Figs. 2–7 also present an interesting phenomenon: the profiles of the u1 and b1 flatten and Hartmann boundary layers along the walls appear as the Hartmann number Ha increases, which are the consequences of the electromagnetic coupling effects.

786

G.-D. Zhang et al. / Computers and Mathematics with Applications 68 (2014) 770–788

Fig. 4. Hartmann flow in 2D: numerical (points) and analytical (line) u1 along x = 0.5.

Fig. 5. Hartmann flow in 2D: numerical (points) and analytical (line) b1 along x = 0.5.

Through this numerical test, we see Iteration I can only compute the MHD problem in the case of low Hartmann number, Iteration II can still work in the case of moderate Hartmann number, only Iteration III can run in the high Hartmann number case. This test validates our theoretical analysis about the three iterative methods in Section 5. 7. Conclusion In this article, a finite element approximation formula for MHD equations on general Lipschitz domain proposed in [27] is considered. Under proper conditions, this scheme provides optimal error estimate. Then, for practical numerical computation, three coupled type iterations for this scheme are proposed and analyzed, meanwhile, the different stability conditions of the three coupling iterations are provided. In our future research, we will study the decoupled iterative methods and compare them with the coupling iterations. In the end, we will give the suitable iterative methods for different MHD problems.

G.-D. Zhang et al. / Computers and Mathematics with Applications 68 (2014) 770–788

787

Fig. 6. Hartmann flow in 2D: numerical (points) and analytical (line) u1 along x = 0.5.

Fig. 7. Hartmann flow in 2D: numerical (points) and analytical (line) b1 along x = 0.5.

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14]

A.E. Lifschitz, Magnetohydrodynamics and Spectral Theory, Kluwer, Dordrecht, 1989. E.R. Priest, A.W. Hood, Advances in Solar System Magnetohydrodynamics, Cambridge University Press, Cambridge, 1991. C. Cao, J. Wu, Two regularity criteria for the 3D MHD equations, J. Differential Equations 248 (2010) 2263–2274. D. Chae, Nonexistence of self-similar singularities in the viscous magnetohydrodynamics with zero resistivity, J. Funct. Anal. 254 (2008) 441–453. Q. Chen, C. Miao, Z. Zhang, On the regularity criterion of weak solution for the 3D viscous magneto-hydrodynamics equations, Comm. Math. Phys. 284 (2008) 919–930. C. He, Y. Wang, On the regularity criteria for weak solutions to the magnetohydrodynamic equations, J. Differential Equations 238 (2007) 1–17. C. He, Z. Xin, On the regularity of weak solutions to the magnetohydrodynamic equations, J. Differential Equations 213 (2005) 235–254. C. He, Z. Xin, Partial regularity of suitable weak solutions to the incompressible magnetohydrodynamic equations, J. Funct. Anal. 227 (2005) 113–152. C. Miao, B. Yuan, Well-posedness of the ideal MHD system in critical Besov spaces, Methods Appl. Anal. 13 (2006) 89–106. C. Miao, B. Yuan, B. Zhang, Well-posedness for the incompressible magneto-hydrodynamic system, Math. Methods Appl. Sci. 30 (2007) 961–976. M.E. Schonbek, T.P. Schonbek, E. Söli, Large-time behaviour of solutions to the magnetohydrodynamics equations, Math. Ann. 304 (1996) 717–756. M. Sermange, R. Temam, Some mathematical questions related to the MHD equations, Comm. Pure Appl. Math. 36 (1983) 635–664. J. Wu, Regularity results for weak solutions of the 3D MHD equations, Discrete Contin. Dyn. Syst. 10 (2004) 543–556. J. Wu, Regularity criteria for the generalized MHD equations, Comm. Partial Differential Equations 33 (2008) 285–306.

788

G.-D. Zhang et al. / Computers and Mathematics with Applications 68 (2014) 770–788

[15] S. Badia, R. Codina, R. Planas, On an unconditionally convergent stabilized finite element approximation of resistive magnetohydrodynamics, J. Comput. Phys. 234 (2013) 399–416. [16] Ľ. Baňas, A. Prohl, Convergent finite element discretization of the multi-fluid nonstationary incompressible magnetihydrodynamics equations, Math. Comp. 79 (2010) 1957–1999. [17] N. Ben Salah, A. Soulaimani, W.G. Habashi, A finite element method for magnetohydrodynamics, Comput. Methods Appl. Mech. Engrg. 190 (2001) 5867–5892. [18] R. Codina, N.H. Silva, Stabilized finite element approximation of the stationary magnetohydrodynamics equations, Comput. Mech. 38 (2006) 344–355. [19] J.F. Gerbeau, A stabilized finite element method for the incompressible magnetohydrodynamic equations, Numer. Math. 87 (2000) 83–111. [20] J.F. Gerbeau, C. Le Bris, T. Lelievre, Mathematical Methods for the Magnetohydrodynamics of Liquid Metals, Oxford University Press, USA, 2006. [21] C. Greif, D. Li, D. Schötzau, X. Wei, A mixed finite element method with exactly divergence-free velocities for incompressible magnetohydrodynamics, Comput. Methods Appl. Mech. Engrg. 199 (2010) 2840–2855. [22] M.D. Gunzburger, A.J. Meir, J.P. Peterson, On the existence, uniqueness, and finite element approximation of solutions of the equations of stationary, incompressible magnetohydrodynamics, Math. Comp. 56 (1991) 523–563. [23] U. Hasler, A. Schneebeli, D. Schötzau, Mixed finite element approximation of incompressible MHD problems based on weighted regularization, Appl. Numer. Math. 51 (2004) 19–45. [24] P. Houston, D. Schötzau, X. Wei, A mixed DG method for linearized incompressible magnetohydrodynamics, J. Sci. Comput. 40 (2009) 281–314. [25] A.J. Meir, P.G. Schmidt, Analysis and numerical approximation of a stationary MHD flow problem with nonideal boundary, SIAM J. Numer. Anal. 36 (1999) 1304–1332. [26] A. Prohl, Convergent finite element discretizations of the nonstationary incompressible magnetohydrodynamics system, ESAIM Math. Model. Numer. Anal. 42 (2008) 1065–1087. [27] D. Schötzau, Mixed finite element methods for stationary incompressible magnetohydrodynamics, Numer. Math. 96 (2004) 771–800. [28] P. Monk, Finite Element Methods for Maxwell’s Equations, Oxford University Press, NewYork, 2003. [29] J.C. Nédélec, Mixed finite elements in R3 , Numer. Math. 35 (1980) 315–341. [30] J.C. Nédélec, Éléments finis mixtes incompressibles pour l’ équation de Stokes dans R3 , Numer. Math. 39 (1982) 97–112. [31] V. Girault, P.A. Raviart, Finite Element Method for Navier–Stokes Equations: Theory and Algorithms, Springer-Verlag, Berlin, Heidelberg, 1986. [32] R. Hiptmair, Finite elements in computational electromagnetism, Acta Numer. 11 (2002) 237–339.