- Email: [email protected]
Error analysis of a fractional-step method for magnetohydrodynamics equations
Journal of Computational and Applied Mathematics 313 (2017) 168–184
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Error analysis of a fractional-step method for magnetohydrodynamics equations Rong An ∗ , Can Zhou College of Mathematics and Information Science, Wenzhou University, Wenzhou 325035, PR China
article
info
Article history: Received 11 November 2015 Received in revised form 21 June 2016 MSC: 65M12 76W05
abstract This paper focuses on a fractional-step finite element method for the magnetohydrodynamics problems in three-dimensional bounded domains. It is shown that the proposed fractional-step scheme allows for a discrete energy identity. A rigorous error analysis is presented. We derive the temporal and spatial error estimates of O (∆t + h) for the velocity and the magnetic field in the discrete space l2 (H1 ) ∩ l∞ (L2 ) under the constraint ∆t ≥ Ch. © 2016 Elsevier B.V. All rights reserved.
Keywords: Magnetohydrodynamics equations Fractional-step method Stability Temporal errors Spatial errors
1. Introduction The incompressible magnetohydrodynamics (MHD) problems are used to describe the flow of a viscous, incompressible and electrically conducting fluid, which are governed by the following time-dependent nonlinear coupled problems:
∂u 1 − 1u + (u · ∇)u + ∇ p + Sb × curl b = f, ∂t Re div u = 0, 1 ∂b + curl (curl b) − curl (u × b) = 0, ∂t Rm div b = 0,
(1.1) (1.2) (1.3) (1.4)
for x ∈ Ω and t ∈ (0, T ) with T > 0, where Ω ⊂ R3 is a bounded and simply-connected domain which is either convex or has a C 1,1 boundary ∂ Ω . Re, Rm and S are three positive constants and denote the Reynolds number, the magnetic Reynolds number and the coupling number, respectively. The vector-value function f represents the body forces applied to the fluid. The MHD problems (1.1)–(1.4) couple the incompressible Navier–Stokes equations with Maxwell’s equations. Thus, the unknowns in (1.1)–(1.4) are the fluid velocity u, the pressure p and the magnetic field b. We refer to Hughes [1] and Moreau [2] for the understanding of the physical background of the MHD problems. To study (1.1)–(1.4), the appropriate initial and boundary conditions are needed. For the sake of simplicity, in this paper, we consider the following initial and
∗
Corresponding author. E-mail addresses: [email protected], [email protected] (R. An).
http://dx.doi.org/10.1016/j.cam.2016.09.005 0377-0427/© 2016 Elsevier B.V. All rights reserved.
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
169
boundary conditions: u(x, 0) = u0 , u = 0,
b · n = 0,
b(x, 0) = b0 curl b × n = 0
in Ω ,
(1.5)
on ∂ Ω × [0, T ],
(1.6)
where n denotes the unit outward normal vector on ∂ Ω . It is necessary to require that u0 and b0 satisfy the compatibility conditions div u0 = 0 and div b0 = 0. The numerical methods for the incompressible MHD problems have received much attention in the last decades. We refer to Gerbeau–Bris–Lelièvre [3] for a review of numerical methods. The mixed finite element approximation was first proposed and studied for the stationary MHD problems in [4], where H1 -conforming elements were used to discretize the magnetic field provided that Ω is either convex or has a C 1,1 boundary ∂ Ω . Inspired by the stabilization method for Stokes problems in [5], a stabilization mixed finite element method for the stationary MHD problems was developed by Gerbeau [6]. For the time-dependent MHD problems (1.1)–(1.6), He proposed a linearized semi-implicit Euler scheme in [7], where L2 -unconditional convergence was proved by using the negative norm technique. For the non-convex domain or Lipschitz polyhedra domain of engineering practice, the magnetic field b may have regularity below H1 (Ω ). In this case, the H1 -conforming finite element discretization for b, albeit stable, may not converge to the correct magnetic field. A mixed finite element formulation based on H(curl)-elements (or Nédélec elements) for b was proposed and studied by Schötzau in [8] for the stationary MHD problems. Other different numerical methods can be found in [9–18] and references cited therein. Roughly speaking, the difficulties in solving the MHD problems numerically are mainly of three kinds: the mixed type of the equations; the incompressible condition and the treatment of the pressure; the nonlinearities of the problems, which are very similar to the incompressible Navier–Stokes problems. In the 1960s, Chorin [19] and Temam [20] proposed a projection method for Navier–Stokes problems, which decoupled the velocity and the pressure in the Navier–Stokes problems. The idea of the projection method is first to compute a velocity field without taking into account incompressibility, and then perform a pressure correction, which is a projection back to the subspace of solenoidal (divergence-free) vector fields. However, the drawback is the appearance of the numerical boundary layer due to the incompatibility of the pressure boundary conditions with those of the original Navier–Stokes problems [21,22]. For the time-dependent MHD problems (1.1)–(1.6), inspired by the projection method in [19,20], Prohl in [14] proposed a projection scheme. However, the projection scheme in [14] does not allow for a discrete energy estimate. To avoid using artificial boundary conditions of pressure type, some fractionalstep schemes for the Navier–Stokes problems were introduced and studied in [23,24]. It is a two-step scheme in which the incompressibility and the nonlinearities of the Navier–Stokes problems are split into different steps, and allows the enforcement of the original boundary conditions in all substeps. In this paper, we propose a two-step fractional-step scheme to approximate the solution of the MHD problems (1.1)–(1.4) with the initial and boundary conditions (1.5)–(1.6). We will prove that the proposed fractional-step scheme allows for a discrete energy identity. To state the main results derived, we introduce the following notations. Let X be a Banach space equipped with norm ∥ · ∥X . Let 0 = t0 < t1 < · · · < tN = T be a uniform partition of the time interval [0, T ] with time step 1t = T /N and tn = n1t for 0 ≤ n ≤ N. We denote two discrete norms by
∥u ∥l2 (X ) = 1t n
N
1/2 un 2X
∥ ∥
,
n=1
∥un ∥l∞ (X ) = max ∥un ∥X . 1≤n≤N
It is proved that the time-discrete fractional-step √ scheme provides the temporal error estimates of O (1t ) for the velocity and the magnetic field in l2 (H1 )∩ l∞ (L2 ) and O ( 1t ) for the pressure in l2 (L2 ). For the fully discrete fractional-step scheme, the finite element error estimates of O (1t + h) for the velocity and the magnetic field in l2 (H1 ) ∩ l∞ (L2 ) are obtained under the constraint 1t ≥ Ch. The remainder of this paper is organized as follows: in the next section, we begin with some notations, lay out some assumptions and recall some known inequalities frequently used. The new linearized projection scheme and the main results are presented in Section 3. Meanwhile, the discrete energy identity is derived in Section 3. The proof containing the main results is given in Sections 4 and 5, which is split into several theorems. 2. Mathematical setting For the mathematical setting of the MHD problems (1.1)–(1.4) with the initial and boundary conditions (1.5)–(1.6), we introduce some function spaces and their associated norms. For k ∈ N+ , 1 ≤ s ≤ +∞, let W k,s (Ω ) denote the standard Sobolev space. The norm in W k,s (Ω ) is defined by
∥u∥W k,s =
1/s β s |∇ u| dx
for 1 ≤ s < +∞,
sup |∇ β u|
for s = +∞,
|β|≤k Ω
|β|≤k Ω
170
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
where ∇ is the differential operator with respect to x and
∇β =
∂ |β| β
β
β
∂ x1 1 ∂ x2 2 ∂ x3 3 k,s
for the multi-index β = (β1 , β2 , β3 ), β1 , β2 , β3 ≥ 0 and |β| = β1 + β2 + β3 . We define W0 (Ω ) to be the subspace of W k,s (Ω ) of functions with zero trace on ∂ Ω . When s = 2, we simply use H k (Ω ) to denote W k,2 (Ω ). The dual space of H01 (Ω ) is denoted by H −1 (Ω ). The boldface Sobolev spaces Hk (Ω ), Wk,s (Ω ) and Ls (Ω ) are used to denote the vector Sobolev spaces H k (Ω )3 , W k,s (Ω )3 and Ls (Ω )3 , respectively. For any u, v ∈ L2 (Ω ), (u, v) = Ω u · vdx denotes L2 (Ω ) inner product. Let X be a Banach space. Its dual space is denoted by X ′ . For some T > 0 or T = ∞, denote Bochner spaces Ls (0, T ; X ), 1 ≤ s < +∞, the spaces of measurable functions from the interval [0, T ] into X such that T
∥v(t )∥sX dt < +∞.
0
If s = +∞, the functions in L∞ (0, T ; X ) are required to satisfy sup ∥v(t )∥X < +∞.
t ∈[0,T ]
The symbols C , C1 , C2 , . . . are used to denote the generic positive constants independent of the time step 1t and the mesh size h. Let us denote H = {u ∈ L2 (Ω ), div u = 0 in Ω , u · n = 0 on ∂ Ω }, V0 = {u ∈ V, div u = 0 in Ω },
V = H10 (Ω ),
X = {u ∈ H (Ω ), u · n = 0 on ∂ Ω },
X0 = {u ∈ X, div u = 0 in Ω },
1
M = L20 (Ω ) =
q ∈ L2 (Ω ),
Ω
qdx = 0 .
Thanks to the Poincaré inequality, the norm in V can be defined by
∥u∥V =
Ω
1/2 |∇ u|2 dx .
Introduce two bilinear forms
a2 (u, v) =
Re
1 Rm
1
a1 (u, v) =
Ω
Ω
∇ u : ∇ vdx ∀ u, v ∈ V,
curl u · curl vdx ∀ u, v ∈ X,
d(v, p) =
pdiv vdx ∀ v ∈ V, p ∈ M
Ω
and a trilinear form b(u, v, w) =
Ω
(u · ∇)v · wdx ∀ u, v, w ∈ H1 (Ω ).
It is easy to verify that if div u = 0 in Ω and u · n = 0 on ∂ Ω , the trilinear form b(u, v, w) satisfies the following important property: b(u, v, w) = −b(u, w, v).
(2.1)
With the above notations, given u0 , b0 ∈ H and f ∈ L (0, T ; L (Ω )), the weak variational formulation of the MHD problems (1.1)–(1.6) is defined to find (u, p, b) ∈ L2 (0, T ; V) × L2 (0, T ; M ) × L2 (0, T ; X) such that ∞
2
(ut , v) + a1 (u, v) + b(u, u, v) + S (b × curl b, v) − d(v, p) = (f, v) ∀ v ∈ V,
(2.2)
and
(bt , w) + a2 (b, w) +
1 Rm
(div b, div w) − (u × b, curl w) = 0 ∀ w ∈ X,
(2.3)
where we use the formula:
curl u · vdx = Ω
u · curl vdx +
Ω
∂Ω
n × u · vds.
(2.4)
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
171
Now, we recall Stokes operator A [25]. Let us denote the orthogonal projection operator by PH from L2 (Ω ) onto H. Then Stokes operator A is defined by
∀ u ∈ V0 ∩ H2 (Ω ) := D(A).
Au = −PH 1u
The emphasis of this paper is to prove some temporal and spatial error estimates of the time-discrete scheme (3.1)–(3.3) and the fully-discrete scheme (3.8)–(3.11) stated in Section 3. Then we make the following assumptions. Assumption 1. The initial data u0 , b0 and f satisfy
∥Au0 ∥L2 + ∥b0 ∥H 2 + sup ∥f(t )∥L2 + ∥ft (t )∥L2 ≤ C3 . 0 ≤t ≤T
Assumption 2. The MHD problems (1.1)–(1.6) admit a unique local strong solution (u, p, b) on [0, T ] such that sup 0≤t ≤T
∥Au(t )∥L2 + ∥b(t )∥H 2 + ∥p(t )∥H 1 + ∥ut (t )∥L2 + ∥bt (t )∥L2 T
+ 0
∥∇ ut ∥2L2 + ∥bt ∥2H 1 + ∥utt ∥2H −1 + ∥btt ∥2H −1 + τ (t )∥utt ∥2L2 + τ (t )∥btt ∥2L2 dt ≤ C4 ,
where τ (t ) = min{t , 1} is a weight function vanishing at t = 0. Assumption 3. For given g ∈ L2 (Ω ), the Stokes problems
−1v + ∇π = g,
div v = 0 v=0
in Ω , on ∂ Ω
admit a unique solution (v, π ) ∈ D(A) × H 1 (Ω ) ∩ M satisfying
∥Av∥L2 + ∥∇π ∥L2 ≤ C ∥g∥L2 . Assumption 4. For given g ∈ L2 (Ω ), the magnetic problems curl curl v = g, v · n = 0,
div v = 0,
in Ω ,
curl v × n = 0,
on ∂ Ω
admit a unique solution v ∈ X0 ∩ H2 (Ω ) satisfying ∥v∥H 2 ≤ C ∥g∥L2 . The following known inequalities are frequently used [26]: 1/2
1/2
∀ v ∈ H1 (Ω ),
(2.5)
1/2
1/2
∀ v ∈ H (Ω ),
(2.6)
∥v∥L3 ≤ C ∥v∥L2 ∥v∥H 1 ∥v∥L∞ ≤ C ∥v∥H 1 ∥v∥H 2
2
∥v∥H 2 ≤ C ∥Av∥L2 ∀ v ∈ D(A),
(2.7)
∥v∥H 2 ≤ C ∥curl curl v∥L2 ∀ v ∈ X0 ∩ H (Ω ), ∥v∥H 1 ≤ C (∥curl v∥L2 + ∥div v∥L2 ) ∀ v ∈ X. 2
(2.8) (2.9)
Finally, we recall a discrete version of Gronwall’s inequality established in [27]. Lemma 2.1. Let ak , bk , ck and γk , for integers k ≥ 0, be the nonnegative numbers such that an + 1 t
n k=0
bk ≤ 1 t
n
γk ak + 1t
k=0
n
ck + B
for n ≥ 0.
(2.10)
k=0
Suppose that 1t γk < 1, for all k, and set σk = (1 − 1t γk )−1 . Then an + 1 t
n k=0
bk ≤ exp 1t
n k=0
γk σk
1t
n
ck + B
for n ≥ 0.
(2.11)
k=0
Remark 2.1. If the first sum on the right in (2.10) extends only up to n − 1, then the estimate (2.11) holds for all k > 0 with σk = 1.
172
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
3. Fractional-step scheme Let 0 = t0 < t1 < · · · < tN = T be a uniform partition of the time interval [0, T ] with time step 1t = T /N and tn = n1t , 0 ≤ n ≤ N, where [0, T ] is the maximal time interval such that a local strong solution exists and satisfies Assumption 2. For any sequence of functions {g n }Nn=0 , we define Dt g n+1 = (g n+1 − g n )/1t for 0 ≤ n ≤ N − 1. Start with u0 = u0 and b0 = b0 . The proposed time discrete fractional-step scheme is the following two-step scheme: Step I: The first step of the fractional-step scheme consists of finding an intermediate velocity un+1/2 and bn+1 such that un+1/2 − un
1t bn+1 − bn
1t
+
−
1 Rm
1 Re
1un+1/2 + (un · ∇)un+1/2 + Sbn × curl bn+1 = f(tn+1 ),
curl (curl bn+1 ) − curl (un+1/2 × bn ) = 0,
div bn+1 = 0
(3.1) (3.2)
with the boundary conditions un+1/2 = 0, bn+1 · n = 0 and curl bn+1 × n = 0 on ∂ Ω . Step II: Given un+1/2 from (3.1), the second step of the fractional-step scheme consists of finding un+1 and pn+1 such that un+1 − un+1/2
1t
−
1 Re
∆(un+1 − un+1/2 ) + ∇ pn+1 = 0,
div un+1 = 0
(3.3)
with the boundary condition un+1 = 0 on ∂ Ω . Since the problems (3.1)–(3.3) consist of three linear problems, their solvabilities should follow easily if we can establish the a priori energy estimates for the solutions to these problems. The following energy identity serves this purpose. Theorem 3.1. Let un+1/2 ∈ V, bn+1 ∈ X0 and (un+1 , pn+1 ) ∈ V0 × M be the weak solutions to the problems (3.1)–(3.3). For 0 ≤ m ≤ N − 1, the following energy identity holds:
∥um+1 ∥2L2 + S ∥bm+1 ∥2L2 +
m 2S 1t
Rm
∥curl bn+1 ∥2L2
n =0
m + (∥un+1 − un+1/2 ∥2L2 + ∥un − un+1/2 ∥2L2 + S ∥bn+1 − bn ∥2L2 ) n=0 m 1t (∥∇ un+1 ∥2L2 + ∥∇ un+1/2 ∥2L2 + ∥∇(un+1 − un+1/2 )∥L2 ) +
Re n=0
= ∥u0 ∥2L2 + S ∥b0 ∥2L2 + 21t
m (f(tn+1 ), un+1/2 ).
(3.4)
n =0
Proof. Multiplying (3.1)–(3.3) by 21tun+1/2 , 2S 1tbn+1 and 21tun+1 , respectively, and adding three equations and using the following formula 2(u − v, u) = ∥u∥2L2 + ∥u − v∥2L2 − ∥v∥2L2
∀ u, v ∈ L2 (Ω ),
(3.5)
and the vector relation
(u × curl v, w) = (w × u, curl v) ∀ u, v, w ∈ H1 (Ω ), we complete the proof of (3.4).
(3.6)
As a direct consequence of (2.9), (3.4) and div bn+1 = 0, we can prove that un+1/2 , un+1 and bn+1 are uniformly bounded in l (H1 ), l∞ (L2 ) ∩ l2 (H1 ) and l∞ (L2 ) ∩ l2 (H1 ), respectively. 2
Corollary 3.1. The fractional-step scheme (3.1)–(3.3) is unconditionally stable. Moreover, there exists some C > 0 such that
∥un+1 ∥l∞ (L2 ) + S ∥bn+1 ∥l∞ (L2 ) + +
1 Re
∥∇ un+1/2 ∥l2 (L2 ) +
S Rm
1 Re
∥∇ un+1 ∥l2 (L2 )
∥bn+1 ∥l2 (H 1 ) ≤ C ,
where C depends on Re, Rm, S, u0 , b0 and f. Next, we give the finite element approximations of (3.1)–(3.3). Let Th be a quasi-uniform partition of Ω into tetrahedra of diameters by h with 0 < h < 1. Let C (Ω ) denote the space of continuous function in Ω = Ω ∪ ∂ Ω . For every K ∈ Th and a nonnegative integer r, Pr (K ) denotes the space of the polynomials on K of degree at most r. We use the boldface spaces
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
173
C(Ω ) and Pr (K ) to denote the vector spaces C (Ω )3 and Pr (K )3 . Define the following conforming finite element subspaces of V, M and X for the velocity, the pressure and the magnetic field, respectively, by Vh = {vh ∈ C(Ω ), vh ∈ P2 (K ), ∀ K ∈ Th } ∩ V, Vh = {vh ∈ C(Ω ), vh ∈ P1 (K ), ∀ K ∈ Th } ∩ V, Xh = {vh ∈ C(Ω ), vh ∈ P1 (K ), ∀ K ∈ Th } ∩ X, Mh = {qh ∈ C (Ω ), qh ∈ P1 (K ), ∀ K ∈ Th } ∩ M . In such a way, the approximation spaces Vh and Mh satisfy the discrete inf–sup condition, i.e., there exists some positive constant β > 0 such that [25]
β∥qh ∥L2 ≤
d(vh , qh )
sup
vh ∈Vh /{0}
∥∇ vh ∥L2
.
(3.7)
Under the above notations, the finite element approximations of the fractional-step scheme (3.1)–(3.3) are described as follows. Start with u0h = Rh u0 , b0h = Qh b0 , where the projectors Rh and Qh are defined in Section 5. For 0 ≤ n ≤ N − 1, n+1/2
∈ Vh and bnh+1 ∈ Xh such that
given unh ∈ Vh and bnh ∈ Xh , find uh
n+1/2
uh
− unh
1t
n+1/2
, vh + a1 (uh
n+1/2
, vh ) + b1 (unh , uh
, vh ) + S (bnh × curl bnh+1 , vh ) = (f(tn+1 ), vh ) ∀ vh ∈ Vh , (3.8)
where b1 (u, v, w) =
1 2
b(u, v, w) −
1 2
b(u, w, v) ∀ u, v, w ∈ V,
and
bnh+1 − bnh
1t
, wh + a2 (bnh+1 , wh ) +
1 Rm
n+1/2
(div bnh+1 , div wh ) − (uh
× bnh , curl wh ) = 0 ∀ wh ∈ Xh .
(3.9)
n+1/2
Then given uh
∈ Vh , find (unh+1 , pnh+1 ) ∈ Vh × Mh such that n+1/2
uhn+1 − uh
1t
n+1/2
, vh + a1 (unh+1 − uh
, vh ) − d(vh , phn+1 ) = 0 ∀ vh ∈ Vh , d(uhn+1 , qh ) = 0
(3.10)
∀ qh ∈ Mh .
By using a similar argument for Theorem 3.1, we can prove that the finite element approximation solutions , phn+1 satisfy the following discrete energy estimate.
bnh+1
(3.11) n+1/2 uh
, uhn+1 ,
n+1/2
Theorem 3.2. For 0 ≤ n ≤ N − 1, let uh ∈ Vh , bnh+1 ∈ Xh and (uhn+1 , pnh+1 ) ∈ Vh × Mh be the solutions to the problems (3.8)–(3.10). Then there exists some C > 0 such that
∥unh+1 ∥l∞ (L2 ) + S ∥bnh+1 ∥l∞ (L2 ) +
1 Re
∥∇ unh+1 ∥l2 (L2 ) +
1 Re
n+1/2
∥∇ uh
∥l2 (L2 ) +
S Rm
∥bnh+1 ∥l2 (H 1 ) ≤ C ,
where C depends on Re, Rm, S, u0 , b0 and f. The emphasis of this paper is to show the following error estimates. Theorem 3.3. Suppose that Assumptions 1–4 are satisfied. Then there exists some τ0 > 0 such that when τ < τ0 , the following temporal error estimates hold:
∥u(tn+1 ) − un+1 ∥l∞ (L2 ) +
1 Re
S ∥b(tn+1 ) − bn+1 ∥l∞ (L2 ) +
∥∇ u(tn+1 ) − ∇ un+1 ∥l2 (L2 ) ≤ C 1t , S
∥b(tn+1 ) − bn+1 ∥l2 (H 1 ) ≤ C 1t , √ ∥p(tn+1 ) − pn+1 ∥l2 (L2 ) ≤ C 1t , √ ∥ τ (tn+1 )(∇ u(tn+1 ) − ∇ un+1 )∥l∞ (L2 ) ≤ C 1t , √ ∥ τ (tn+1 )(b(tn+1 ) − bn+1 )∥l∞ (H 1 ) ≤ C 1t , Rm
(3.12) (3.13) (3.14) (3.15) (3.16)
where C depends on Re, Rm and S, and un+1 , pn+1 , bn+1 are the numerical solutions to the problems (3.1)–(3.3) at time-step n + 1, and u(tn+1 ), p(tn+1 ), b(tn+1 ) are the exact solutions evaluated at time tn+1 , and τ (tn+1 ) = min{1, tn+1 } for 0 ≤ n ≤ N − 1.
174
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
Theorem 3.4. Suppose that Assumptions 1–4 are satisfied. If the time step 1t satisfies 1t ≥ Ch for some C > 0, then there exist some τ0 > 0 and h0 > 0 such that when τ < τ0 and h < h0 , the following error estimates hold:
∥u(tn+1 ) − unh+1 ∥l∞ (L2 ) + S ∥b(tn+1 ) − bnh+1 ∥l∞ (L2 ) ≤ C (1t + h), 1 Re
∥∇(u(tn+1 ) − unh+1 )∥l2 (L2 ) +
S Rm
∥b(tn+1 ) − bnh+1 ∥l2 (H 1 ) ≤ C (1t + h),
(3.17) (3.18)
where C depends on Re, Rm and S, and unh+1 , pnh+1 , bnh+1 are the finite element approximation solutions at time-step n + 1 for 0 ≤ n ≤ N − 1. The proof of Theorems 3.3 and 3.4 is split into several lemmas in Sections 4 and 5. 4. Temporal error analysis In this section, we will prove Theorem 3.3. For 0 ≤ n ≤ N − 1, let us denote errors by enu+1/2 = u(tn+1 ) − un+1/2 ,
eun+1 = u(tn+1 ) − un+1 , ebn+1
= b(tn+1 ) − b
n +1
,
enp+1 = p(tn+1 ) − pn+1 .
Taking t = tn+1 in (1.1)–(1.4) yields u(tn+1 ) − u(tn )
−
1t
1 Re
1u(tn+1 ) + (u(tn+1 ) · ∇)u(tn+1 )
+ ∇ p(tn+1 ) + Sb(tn+1 ) × curl b(tn+1 ) = f(tn+1 ) + Rnu+1 , div u(tn+1 ) = 0, b(tn+1 ) − b(tn )
1t
(4.1)
div b(tn+1 ) = 0,
+
1 Rm
(4.2)
curl (curl b(tn+1 )) − curl (u(tn+1 ) × b(tn+1 )) = Rnb+1 ,
(4.3)
where Rnu+1 = −
tn+1
1
1t
(t − tn )utt (t )dt ,
Rnb+1 = −
tn
1
1t
tn+1
(t − tn )btt (t )dt .
tn
From Assumption 2, one has
∥Rnu+1 ∥l2 (H −1 ) + ∥Rnb+1 ∥l2 (H −1 ) ≤ C 1t , ∥ τ (tn+1 )Run+1 ∥l2 (L2 ) + ∥ τ (tn+1 )Rnb+1 ∥l2 (L2 ) ≤ C 1t .
(4.4) (4.5)
Subtracting (3.1) and (3.2) from (4.1) and (4.3), respectively, we obtain n+1/2
eu
− enu
1t
−
1 Re
1enu+1/2 + ((u(tn+1 ) − u(tn )) · ∇)u(tn+1 ) + (enu · ∇)u(tn+1 ) + (un · ∇)eun+1/2
+ S (b(tn+1 ) − b(tn )) × curl b(tn+1 ) + Senb × curl b(tn+1 ) + Sbn × curl ebn+1 + ∇ p(tn+1 ) = Run+1
(4.6)
and ebn+1 − enb
1t
+
1 Rm
curl curl enb+1 − curl (enu+1/2 × bn )
− curl (u(tn+1 ) × enb ) − curl (u(tn+1 ) × (b(tn+1 ) − b(tn ))) = Rbn+1 .
(4.7)
√ First, we show that the fractional-step scheme (3.1)–(3.3) provides the temporal error estimates of O ( 1t ) for the velocity and the magnetic field in l2 (H1 ) ∩ l∞ (L2 ). Lemma 4.1. Suppose that Assumptions 1–2 are satisfied. Then there exists some positive constant C > 0 such that for every 0 ≤ m ≤ N − 1, +1 2 +1 2 ∥em ∥L2 + S ∥em ∥L2 + u b
m S 1t
Rm
n =0
∥curl enb+1 ∥2L2
+
m (∥enu+1/2 − enu ∥2L2 + ∥eun+1/2 − eun+1 ∥2L2 + S ∥ebn+1 − enb ∥2L2 ) n =0
m 1t + (∥∇ enu+1 ∥2L2 + ∥∇ enu+1/2 ∥2L2 + ∥∇(enu+1 − eun+1/2 )∥2L2 ) ≤ C 1t ,
Re n=0
where C depends on Re, Rm and S.
(4.8)
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184 n+1/2
Proof. Testing (4.6) and (4.7) by 21teu to
175
and 2S 1tenb+1 , respectively, using (3.6) and adding two resulting equations lead
∥eun+1/2 ∥2L2 + ∥enu+1/2 − enu ∥2L2 − ∥enu ∥2L2 + S ∥enb+1 ∥2L2 21t 2S 1t + S ∥enb+1 − enb ∥2L2 − S ∥enb ∥2L2 + ∥∇ enu+1/2 ∥2L2 + ∥curl ebn+1 ∥2L2 Re
≤ 21t |b(u(tn+1 ) − u(tn ), u(tn+1 ),
enu+1/2
Rm
) + ( , u(tn+1 ), eun+1/2 )| b enu enb
+ 2S 1t |((b(tn+1 ) − b(tn )) × curl b(tn+1 ) +
× curl b(tn+1 ), eun+1/2 )|
+ 2S 1t |(u(tn+1 ) × enb + u(tn+1 ) × (b(tn+1 ) − b(tn )), curl ebn+1 )| + 21t |(∇ p(tn+1 ), enu+1/2 )| + 21t |(Rnu+1 , enu+1/2 ) + S (Rbn+1 , enb+1 )| = I1 + · · · + I5 .
(4.9)
By using Assumption 2 and Young’s inequality, we have I1 = 21t |b(u(tn+1 ) − u(tn ), u(tn+1 ), enu+1/2 ) + b(enu , u(tn+1 ), enu+1/2 )|
≤ C 1t (∥u(tn+1 ) − u(tn )∥L2 + ∥enu ∥L2 )∥Au(tn+1 )∥L2 ∥∇ eun+1/2 ∥L2 1t ∥∇ eun+1/2 ∥2L2 + C Re1t (∥u(tn+1 ) − u(tn )∥2L2 + ∥enu ∥2L2 ). ≤ 8Re
By a similar method, we have I2 = 2S 1t |((b(tn+1 ) − b(tn )) × curl b(tn+1 ) + enb × curl b(tn+1 ), eun+1/2 )|
≤ CS 1t (∥b(tn+1 ) − b(tn )∥L2 + ∥enb ∥L2 )∥b(tn+1 )∥H 2 ∥∇ eun+1/2 ∥L2 1t ≤ ∥∇ eun+1/2 ∥2L2 + C ReS 2 1t (∥b(tn+1 ) − b(tn )∥2L2 + ∥enb ∥2L2 ) 8Re
and I3 = 2S 1t |(u(tn+1 ) × enb + u(tn+1 ) × (b(tn+1 ) − b(tn )), curl ebn+1 )|
≤ CS 1t (∥b(tn+1 ) − b(tn )∥L2 + ∥enb ∥L2 )∥Au(tn+1 )∥L2 ∥curl enb+1 ∥L2 S 1t ∥curl enb+1 ∥2L2 + CRmS 1t (∥b(tn+1 ) − b(tn )∥2L2 + ∥enb ∥2L2 ). ≤ 2Rm
Due to the fact that div enu = 0 in Ω and enu = 0 on ∂ Ω , we have I4 = 21t |(∇ p(tn+1 ), enu+1/2 )|
= 21t |(∇ p(tn+1 ), enu+1/2 − enu )| ≤ C 1t ∥enu+1/2 − enu ∥L2 ∥∇ p(tn+1 )∥L2 ≤
1 2
∥eun+1/2 − enu ∥2L2 + C (1t )2 .
From (2.9) and Young’s inequality, one has I5 = 21t |(Rnu+1 , enu+1/2 ) + S (Rnb+1 , enb+1 )|
≤ C 1t (∥Rnu+1 ∥H −1 ∥∇ enu+1/2 ∥L2 + S ∥Rnb+1 ∥H −1 ∥curl ebn+1 ∥L2 ) 1t S 1t ≤ ∥∇ eun+1/2 ∥2L2 + ∥curl enb+1 ∥2L2 + C 1t (Re∥Run+1 ∥2H −1 + RmS ∥Rbn+1 ∥2H −1 ). 8Re
2Rm
n+1/2
Using the definitions of enu+1 and eu n+1/2
eun+1 − eu
1t
−
1 Re
, we rewrite (3.3) as
∆(enu+1 − enu+1/2 ) − ∇ pn+1 = 0,
div enu+1 = 0.
(4.10)
Testing (4.10) by 21teun+1 leads to
∥enu+1 ∥2L2 + ∥enu+1 − enu+1/2 ∥2L2 − ∥enu+1/2 ∥2L2 +
1t Re
(∥∇ eun+1 ∥2L2 + ∥∇(eun+1 − eun+1/2 )∥2L2 − ∥∇ enu+1/2 ∥2L2 ) = 0.
(4.11)
Adding (4.9) and (4.11) and combining these estimates for I1 to I5 into the resulting inequality and using discrete Gronwall’s inequality (2.11), we conclude that there exists some positive constant C > 0 depending on Re, Rm, S such that (4.8) holds.
176
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
Remark 4.1. Assumption 2 and (4.8) imply that un+1/2 , un+1 and bn+1 are uniform boundedness in l∞ (H1 ):
∥∇ un+1 ∥l∞ (L2 ) + ∥∇ un+1/2 ∥l∞ (L2 ) + ∥bn+1 ∥l∞ (H 1 ) ≤ C ,
(4.12)
where C depends on Re, Rm and S. Using the above uniform boundedness (4.12), we can improve the orders of the temporal error estimates derived in Lemma 4.1 as follows. Lemma 4.2. Suppose that Assumptions 1–2 are satisfied. Then there exists some positive constant C > 0 such that
∥eun+1 ∥l∞ (L2 ) + S ∥enb+1 ∥l∞ (L2 ) +
1 Re
∥∇ enu+1 ∥l2 (L2 ) + ∥
S Rm
∥enb+1 ∥l2 (H 1 ) ≤ C 1t ,
eun+1
−
enu 2l2 (L2 )
∥
(4.13)
≤ C (1t ) , 3
(4.14)
where C depends on Re, Rm and S. Proof. Taking the sum of (4.6) and (4.10) gives eun+1 − enu
1t
−
1 Re
1enu+1 + ((u(tn+1 ) − u(tn )) · ∇)u(tn+1 ) + (enu · ∇)un+1/2 + (u(tn ) · ∇)eun+1/2
+ S (b(tn+1 ) − b(tn )) × curl b(tn+1 ) + Senb × curl bn+1 + Sb(tn ) × curl ebn+1 + ∇ epn+1 = Run+1 .
(4.15)
Subtracting (3.2) from (4.3) and using u(tn+1 ) × b(tn+1 ) − un+1/2 × bn = u(tn+1 ) × (b(tn+1 ) − b(tn )) + eun+1/2 × b(tn ) + un+1/2 × enb yields ebn+1 − enb
1t
+
1 Rm
curl curl enb+1 − curl (enu+1/2 × b(tn )) − curl (un+1/2 × enb )
− curl (u(tn+1 ) × (b(tn+1 ) − b(tn ))) = Rnb+1
(4.16)
and div ebn+1 = 0. Testing (4.15) and (4.16) by 21tenu+1 and 2S 1tebn+1 , respectively, using (3.6) and adding two resulting equations leads to
∥eun+1 ∥2L2 + ∥enu+1 − enu ∥2L2 − ∥enu ∥2L2 + S ∥enb+1 ∥2L2 2S 1t 21t + S ∥ebn+1 − enb ∥2L2 − S ∥enb ∥2L2 + ∥∇ enu+1 ∥2L2 + ∥curl enb+1 ∥2L2 Re
Rm
≤ 21t |b(u(tn+1 ) − u(tn ), u(tn+1 ), enu+1 )| + 2S 1t |((b(tn+1 ) − b(tn )) × curl b(tn+1 ), eun+1 )| + 2S 1t |(u(tn+1 ) × (b(tn+1 ) − b(tn )), curl enb+1 )| + 21t |(Run+1 , eun+1 ) + S (Rbn+1 , ebn+1 )| + 21t |b(u(tn ), enu+1/2 , enu+1 ) + b(enu , un+1/2 , enu+1 )| + 2S 1t |(un+1/2 × enb , curl enb+1 ) + (enb × curl bn+1 , enu+1 )| + 2S 1t |(b(tn ) × curl enb+1 , enu+1 − enu+1/2 )| = I6 + · · · + I12 .
(4.17)
By a similar method in the proof of Lemma 4.1, we can bound I6 + · · · + I9 as I6 + I7 + I8 + I9
= 21t |b(u(tn+1 ) − u(tn ), u(tn+1 ), enu+1 )| + 2S 1t |((b(tn+1 ) − b(tn )) × curl b(tn+1 ), enu+1 )| + 2S 1t |(u(tn+1 ) × (b(tn+1 ) − b(tn )), curl enb+1 )| + 21t |(Run+1 , eun+1 ) + S (Rbn+1 , ebn+1 )| S 1t 1t ∥∇ eun+1 ∥2L2 + ∥curl enb+1 ∥2L2 + C 1t Re∥u(tn+1 ) − u(tn )∥2L2 + (ReS 2 + RmS )∥b(tn+1 ) − b(tn )∥2L2 ≤ 4Re 4Rm + C 1t Re∥Rnu+1 ∥2H −1 + RmS ∥Rnb+1 ∥2H −1 . From (2.5)–(2.7), we estimate I10 as I10 = 21t |b(u(tn ), enu+1/2 , enu+1 ) + b(enu , un+1/2 , eun+1 )| 1/2
1/2
≤ C 1t ∥Au(tn )∥L2 ∥∇ enu+1 ∥L2 (∥enu+1/2 − enu ∥L2 + ∥enu ∥L2 ) + C 1t ∥∇ un+1/2 ∥L2 ∥∇ eun+1 ∥L2 ∥enu ∥L2 ∥∇ enu ∥L2 1t 1t ≤ ∥∇ enu+1 ∥2L2 + ∥∇ enu ∥2L2 + C Re1t (∥eun+1/2 − enu ∥2L2 + ∥enu ∥2L2 ). 4Re
2Re
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
177
Similarly, from (2.9), I11 is bounded by I11 = 2S 1t |(un+1/2 × enb , curl enb+1 ) + (enb × curl bn+1 , enu+1 )| 1/2
1/2
1/2
1/2
≤ CS 1t ∥∇ un+1/2 ∥L2 ∥curl enb+1 ∥L2 ∥enb ∥L2 ∥curl enb ∥L2 + CS 1t ∥curl bn+1 ∥L2 ∥∇ enu+1 ∥L2 ∥enb ∥L2 ∥curl enb ∥L2 S 1t S 1t 1t ∥∇ eun+1 ∥2L2 + ∥curl enb+1 ∥2L2 + ∥curl enb ∥2L2 + C (Re2 RmS 3 + Rm3 S )1t ∥enb ∥2L2 . ≤ 4Re
4Rm
2Rm
The last term I12 satisfies I12 = 2S 1t |(b(tn ) × curl enb+1 , enu+1 − enu+1/2 )|
≤ CS 1t ∥b(tn )∥H 2 ∥curl enb+1 ∥L2 ∥enu+1 − enu+1/2 ∥L2 S 1t ≤ ∥curl enb+1 ∥2L2 + CRmS 1t ∥enu+1 − enu+1/2 ∥2L2 . 4Rm
Combining these estimates into (4.17) and using the discrete Gronwall’s inequality (2.11), (4.4), (4.8), we complete the proof of (4.13) and (4.14). The temporal error estimate for the pressure in l2 (L2 ) is given in the next lemma by using the inf–sup condition. Lemma 4.3. Suppose that Assumptions 1–2 are satisfied. Then there exists some positive constant C > 0 such that
√ ∥epn+1 ∥l2 (L2 ) ≤ C 1t ,
(4.18)
where C depends on Re, Rm and S. Proof. From (4.15) and inf–sup condition, one has
∥epn+1 ∥L2 ≤ C sup
d(v, enp+1 )
v∈V
≤
∥∇ v∥L2
C
∥en+1 − enu ∥L2 + C ∥∇ enu+1 ∥L2 + C ∥Run+1 ∥H −1 + C (∥∇ enu ∥L2 + ∥∇ eun+1/2 ∥L2 + ∥u(tn+1 ) − u(tn )∥L2 ) 1t u + C (∥curl enb ∥L2 + ∥curl enb+1 ∥L2 + ∥b(tn+1 ) − b(tn )∥L2 ).
Taking the sum of the above inequality from 0 to N − 1 and using (4.4), (4.8) and (4.14), we complete the proof of Lemma 4.3. As a consequence of (4.14) and Assumption 2, we have
∥un+1 − un ∥l∞ (L2 ) + ∥bn+1 − bn ∥l∞ (L2 ) ≤ C 1t .
(4.19)
To give the spatial error estimates for the fully discrete fractional-step scheme (3.8)–(3.11), we need to show the regularity results for un+1/2 , un+1 , pn+1 and bn+1 . Lemma 4.4. Suppose that Assumptions 1–4 are satisfied. Then there exists some τ0 > 0 such that when 1t < τ0 , there holds 1 Re
∥Aun+1 ∥l2 (L2 ) +
1 Re
∥un+1/2 ∥l2 (H 2 ) + ∥∇ pn+1 ∥l2 (L2 ) +
1 Rm
∥bn+1 ∥l2 (H 2 ) ≤ C ,
(4.20)
where C depends on Re, Rm and S. Proof. Taking the sum of (3.1) and (3.3) leads to
−
1 Re
1un+1 + ∇ pn+1 = f(tn+1 ) − Dt un+1 − (un · ∇)un+1/2 − Sbn × curl bn+1 ,
(4.21)
div un+1 = 0.
(4.22)
Using Assumption 3, we get 1 Re
∥Aun+1 ∥L2 + ∥∇ pn+1 ∥L2 ≤ C (∥Dt un+1 ∥L2 + ∥f(tn+1 )∥L2 + S ∥bn × curl bn+1 ∥L2 + ∥(un · ∇)un+1/2 ∥L2 ).
The last two terms in the above inequality can be bounded, respectively, by S ∥bn × curl bn+1 ∥L2 ≤ S ∥enb × curl bn+1 ∥L2 + S ∥b(tn ) × curl bn+1 ∥L2
≤ CS (∥curl enb ∥L2 ∥bn+1 ∥H 2 + ∥curl bn+1 ∥L2 )
(4.23)
178
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
and 1/2
1/2
∥(un · ∇)un+1/2 ∥L2 ≤ C ∥∇ un ∥L2 ∥Aun ∥L2 ∥∇ un+1/2 ∥L2 ≤
1 2Re
∥Aun ∥L2 + C Re∥∇ un ∥L2 ,
where we use (4.12). Similarly, from Assumption 4 and (3.2) and the following formula curl (u × v) = (div v)u − (u · ∇)v + (v · ∇)u − (div u)v,
(4.24)
we have 1 Rm
∥bn+1 ∥H 2 ≤ C (∥Dt bn+1 ∥L2 + ∥(un+1/2 · ∇)bn ∥L2 ) + C (∥(bn · ∇)un+1/2 ∥L2 + ∥(div un+1/2 )bn ∥L2 ) ≤ C (∥Dt bn+1 ∥L2 + ∥∇ bn ∥L3 ∥∇ un+1/2 ∥L2 + ∥bn ∥L∞ ∥∇ un+1/2 ∥L2 ) 1/2
1/2
≤ C (∥Dt bn+1 ∥L2 + ∥curl bn ∥L2 ∥bn ∥H 2 ) ≤ C (∥Dt bn+1 ∥L2 + Rm∥curl bn ∥L2 ) +
1 2Rm
∥bn ∥H 2 .
(4.25)
Combining (4.23) and (4.25), and taking the sum from 0 to N − 1 yield
1t
N −1 1 n =0
≤ C 1t
Re
∥Aun+1 ∥2L2 + ∥∇ pn+1 ∥2L2 + 2
N −1
1
∥bn+1 ∥2H 2 2
Rm
∥Dt un+1 ∥2L2 + ∥Dt bn+1 ∥2L2 + ∥f(tn+1 )∥2L2 + C 1t ∥curl bn+1 ∥2L2 + ∥curl bn ∥2L2 + ∥∇ un ∥2L2
n =0
+ CS 2 1t
N −1
∥curl enb ∥2L2 ∥bn+1 ∥2H 2 +
n =0
N −1 1 1t
4
n=0
∥Aun ∥2L2 + 2
Re
1
n 2 ∥ b ∥ . H2 2
Rm
The last term in the above inequality can be absorbed on the left-hand side. From (4.13) and (4.19), we have
1t
N −1 1 n =0
Let τ0 = 1 Re
Re
1 . 2C6 S 2 Rm2
∥Aun+1 ∥2L2 + ∥∇ pn+1 ∥2L2 + 2
1
∥bn+1 ∥2H 2 2
Rm
≤ C + C6 S 2 (1t )2
N −1
∥bn+1 ∥2H 2 .
n =0
When 1t < τ0 , the above inequality implies that
∥Aun+1 ∥l2 (L2 ) + ∥∇ pn+1 ∥l2 (L2 ) +
1 Rm
∥bn+1 ∥l2 (H 2 ) ≤ C .
(4.26)
By using the regularity result [28] for elliptic problem and (4.26), we can easily get 1 Re
∥un+1/2 ∥l2 (H 2 ) ≤ C ,
if we notice (4.8) and
∥un+1/2 − un ∥L2 ≤ ∥enu+1/2 − enu ∥L2 + ∥u(tn+1 ) − u(tn )∥L2 . We complete the proof of (4.20).
Remark 4.2. From (4.20) and Assumption 2, we have 1 Re
∥Aeun+1 ∥l2 (L2 ) + ∥∇ enp+1 ∥l2 (L2 ) +
1 Rm
∥enb+1 ∥l2 (H 2 ) ≤ C .
(4.27)
Lemma 4.5. Suppose that Assumptions 1–4 are satisfied. Then when τ < τ0 , there holds
√ ∥ τ (tn+1 )∇ eun+1 ∥l∞ (L2 ) + ∥ τ (tn+1 )curl enb+1 ∥l∞ (L2 ) ≤ C 1t , where τ0 is defined in Lemma 4.4, and C depends on Re, Rm and S.
(4.28)
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
179
Proof. Testing (4.15) by −21tAenu+1 gives
∥∇ enu+1 ∥2L2 − ∥∇ enu ∥2L2 + ∥∇(enu+1 − enu )∥2L2 +
1t Re
∥Aenu+1 ∥2L2
≤ C 1t ∥Run+1 ∥2L2 + C 1t ∥((u(tn+1 ) − u(tn )) · ∇)u(tn+1 )∥2L2 + C 1t ∥(enu · ∇)un+1/2 + (u(tn ) · ∇)enu+1/2 ∥2L2 + C 1t ∥(b(tn+1 ) − b(tn )) × curl b(tn+1 )∥2L2 + C 1t ∥enb × curl bn+1 + b(tn ) × curl enb+1 ∥2L2 = C 1t ∥Rnu+1 ∥2L2 + I13 + · · · + I16 .
(4.29)
From Assumption 2, one has I13 = C 1t ∥((u(tn+1 ) − u(tn )) · ∇)u(tn+1 )∥2L2
≤ C 1t ∥Au(tn+1 )∥2L2 ∥∇(u(tn+1 ) − u(tn ))∥2L2 tn+1 ∥∇ ut ∥2L2 dt . ≤ C (1t )2 tn
Using (2.6) and (4.12), we get I14 = C 1t ∥(enu · ∇)un+1/2 + (u(tn ) · ∇)enu+1/2 ∥2L2
≤ C 1t (∥enu ∥2L∞ ∥∇ un+1/2 ∥2L2 + ∥Au(tn )∥2L2 ∥∇ enu+1/2 ∥2L2 ) ≤ C 1t (∥∇ enu ∥L2 ∥Aenu ∥L2 + ∥∇ enu+1/2 ∥2L2 ) ≤
1t 2Re
∥Aenu ∥2L2 + C 1t (∥∇ enu ∥2L2 + ∥∇ enu+1/2 ∥2L2 ).
By a similar argument for the above two inequalities, we have I15 = C 1t ∥(b(tn+1 ) − b(tn )) × curl b(tn+1 )∥2L2
≤ C 1t ∥b(tn+1 )∥2H 2 ∥b(tn+1 ) − b(tn )∥2H 1 tn+1 2 ≤ C (1t ) ∥bt ∥2H 1 dt , tn
and I16 = C 1t ∥enb × curl bn+1 + b(tn ) × curl enb+1 ∥2L2
≤ C 1t (∥enb ∥2L∞ ∥curl bn+1 ∥2L2 + ∥b(tn )∥2H 2 ∥curl enb+1 ∥2L2 ) ≤
C7 1t 2Rm
∥enb ∥2H 2 + C 1t (∥curl enb ∥2L2 + ∥curl enb+1 ∥2L2 ).
Combining these estimates into (4.29) yields
∥∇ eun+1 ∥2L2 − ∥∇ enu ∥2L2 + ≤ C 1t ∥Run+1 ∥2L2 + + C (1t )
2
1t 2Re
tn+1
1t Re
∥Aenu+1 ∥2L2
∥Aenu ∥2L2 + ut 2L2 dt
∥∇ ∥ tn
C7 1 t 2Rm
tn+1
+
∥enb ∥2H 2 + C 1t (∥∇ enu ∥2L2 + ∥∇ eun+1/2 ∥2L2 + ∥curl enb ∥2L2 + ∥curl enb+1 ∥2L2 ) bt 2H 1 dt
∥ ∥ tn
.
(4.30)
Testing (4.7) by 21tcurl curl enb+1 leads to
∥curl ebn+1 ∥2L2 − ∥curl enb ∥2L2 +
C7 1 t Rm
∥enb+1 ∥2H 2
≤ C 1t ∥Rbn+1 ∥2L2 + C 1t ∥curl (enu+1/2 × bn )∥2L2 + C 1t ∥curl (u(tn+1 ) × enb )∥2L2 + C 1t ∥curl (u(tn+1 ) × (b(tn+1 ) − b(tn )))∥2L2 = C 1t ∥Rnb+1 ∥2L2 + I17 + I18 + I19 .
(4.31)
180
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
Using (2.5), (2.6) and (4.24), we obtain I17 = C 1t ∥curl (enu+1/2 × bn )∥2L2
≤ C 1t (∥(enu+1/2 · ∇)bn ∥2L2 + ∥(bn · ∇)enu+1/2 ∥2L2 + ∥(div eun+1/2 )bn ∥2L2 ) ≤ C 1t (∥∇ enu+1/2 ∥2L2 ∥enb ∥H 2 ∥curl enb ∥L2 + ∥∇ eun+1/2 ∥2L2 ∥b(tn )∥2H 2 ) ≤ C (1t )2 ∥enb ∥2H 2 + C ∥curl enb ∥2L2 + C 1t ∥∇ enu+1/2 ∥2L2 . Similarly, we have I18 = C 1t ∥curl (u(tn+1 ) × enb )∥2L2
≤ C 1t (∥(u(tn+1 ) · ∇)enb ∥2L2 + ∥(enb · ∇)u(tn+1 )∥2L2 ) ≤ C 1t ∥curl enb ∥2L2 ∥Au(tn+1 )∥2L2 ≤ C 1t ∥curl enb ∥2L2 , and I19 = C 1t ∥curl (u(tn+1 ) × (b(tn+1 ) − b(tn )))∥2L2
≤ C 1t (∥(u(tn+1 ) · ∇)(b(tn+1 ) − b(tn ))∥2L2 + ∥((b(tn+1 ) − b(tn )) · ∇)u(tn+1 )∥2L2 ) ≤ C 1t ∥b(tn+1 ) − b(tn )∥2H 1 ∥Au(tn+1 )∥2L2 tn+1 ≤ C (1t )2 ∥bt ∥2H 1 dt . tn
Taking the sum of (4.30) and (4.31) and combining the above estimates into the resulting inequality give
∥∇ eun+1 ∥2L2 + ∥curl enb+1 ∥2L2 +
1t Re
∥Aenu+1 ∥2L2 +
C7 1t Rm
∥enb+1 ∥2H 2
≤ ∥∇ enu ∥2L2 + ∥curl enb ∥2L2 + C 1t (∥Rnu+1 ∥2L2 + ∥Rnb+1 ∥2L2 ) +
1t 2Re
∥Aenu ∥2L2 +
C7 1t 2Rm
∥enb ∥2H 2
+ C 1t 2 ∥enb ∥2H 2 + C ∥curl enb ∥2L2 + C 1t (∥∇ enu ∥2L2 + ∥∇ enu+1/2 ∥2L2 + ∥curl enb ∥2L2 + ∥curl ebn+1 ∥2L2 ) tn+1 tn+1 + C (1t )2 ∥∇ ut ∥2L2 dt + ∥bt ∥2H 1 dt . tn
tn
Multiplying it by τ (tn+1 ) and using τ (tn+1 ) ≤ τ (tn ) + 1t and |τ (tn+1 )| ≤ 1, we derive that
τ (tn+1 )∥∇ eun+1 ∥2L2 + τ (tn+1 )∥curl enb+1 ∥2L2 +
1t Re
τ (tn+1 )∥Aenu+1 ∥2L2 +
C7 1 t Rm
τ (tn+1 )∥ebn+1 ∥2H 2
≤ τ (tn )∥∇ enu ∥2L2 + τ (tn )∥curl enb ∥2L2 + 1t (∥∇ enu ∥2L2 + ∥curl enb ∥2L2 ) + C 1t (τ (tn+1 )∥Run+1 ∥2L2 + τ (tn+1 )∥Rbn+1 ∥2L2 ) S 1t n 2 1t τ (tn )∥Aenu ∥2L2 + τ (tn ) ∥eb ∥H 2 + C (1t )2 (∥Aenu ∥2L2 + ∥enb ∥2H 2 ) + 4Re
4Rm
+ C 1t (∥∇ enu ∥2L2 + ∥∇ enu+1/2 ∥2L2 + ∥curl enb ∥2L2 + ∥curl ebn+1 ∥2L2 ) tn+1 tn+1 + C (1t )2 ∥∇ ut ∥2L2 dt + ∥bt ∥2H 1 dt + C ∥curl enb ∥2L2 . tn
Taking the sum of (4.32) from 0 to m ≤ N − 1, and using (4.5), (4.8), (4.13), (4.20), we complete the proof of (4.28). Combining the results in this section, we obtain the temporal error estimates (3.12)–(3.16). 5. Spatial error analysis To show the spatial error estimates, we introduce the projectors (Rh , πh ) : V × M −→ Vh × Mh defined by a1 (Rh u − u, vh ) − d(vh , πh p − p) = 0
∀ vh ∈ Vh , d(Rh u − u, qh ) = 0 ∀ qh ∈ Mh .
In addition, define Kh : V −→ Vh and Qh : X −→ Xh by
(∇(Kh u − u), ∇ vh ) = 0 ∀ vh ∈ Vh , a2 (Qh b − b, wh ) +
(4.32)
tn
1 Rm
(div (Qh b − b), div wh ) = 0 ∀ wh ∈ Xh .
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
181
Then for all (u, p, b) ∈ V ∩ H2 (Ω ) × M ∩ H 1 (Ω ) × X ∩ H2 (Ω ), the following approximation and stability properties hold:
∥u − Rh u∥L2 + h∥∇(u − Rh u)∥L2 ≤ Ch2 ∥u∥H 2 ,
(5.1)
∥u − Kh u∥L2 + h∥∇(u − Kh u)∥L2 ≤ Ch ∥u∥H 2 ,
(5.2)
∥b − Qh b∥L2 + h∥b − Qh b∥H 1 ≤ Ch ∥b∥H 2 , ∥p − πh p∥L2 ≤ Ch∥p∥H 1 , ∥∇ Rh u∥L6 + ∥∇ Kh u∥L6 + ∥Qh b∥W 1,6 + ∥πh p∥L6 ≤ C .
(5.4)
2 2
(5.3) (5.5)
For 0 ≤ n ≤ N − 1, let us denote errors by Enu+1/2 = Kh un+1/2 − un+1/2 ,
Enu = Rh un − un ,
n+1/2
enuh = Rh un − unh ,
euh
n+1/2
= Kh un+1/2 − uh
Enb = Qh bn − bn ,
,
enbh = Qh bn − bnh ,
n+1/2
From (3.1)–(3.3) and (3.8)–(3.11), we obtain that euh
n+1/2
− enuh
euh
1t
n+1/2
, vh + a1 (euh
, vh ) =
n+1/2
Eu
Epn = πh pn − pn , enph = πh pn − pnh .
n +1 , enbh+1 , (euh , enph+1 ) satisfy the following problems:
− Enu
1t
n+1/2
, vh + b1 (unh , uh
, vh )
− b1 (un , un+1/2 , vh ) + S (bnh × curl bnh+1 , vh ) − S (bn × curl bn+1 , vh ) ∀ vh ∈ Vh ,
(5.6)
and
n+1 ebh − enbh
1t
, wh +
(
a2 enbh+1 n+1/2
+ (un+1/2 × bn − uh
, wh ) +
1 Rm
(div
enbh+1
, wh ) =
Enb+1 − Enb
1t
, wh
× bnh , curl wh ) ∀ wh ∈ Xh ,
(5.7)
and
n+1/2
n+1 euh − euh
=
n+1/2
, vh + a1 (enuh+1 − euh
1t
n+1/2
Eun+1 − Eu
1t
n +1 , vh ) − d(vh , eph )
, vh − a1 (Enu+1/2 , vh ) ∀ vh ∈ Vh ,
(5.8)
and n +1 , qh ) = 0, d(euh
∀ qh ∈ Mh .
(5.9)
The main result in this section is to show the following spatial error estimate based on the regularity estimate (4.20). Lemma 5.1. Suppose that Assumptions 1–4 are satisfied. If the time step 1t satisfies 1t ≥ Ch, then there exists some h0 > 0 such that when h < h0 , the following spatial error estimate holds: +1 ∥enuh+1 ∥l∞ (L2 ) + S ∥em bh ∥l∞ (L2 ) +
1 Re
∥∇ enuh+1 ∥l2 (L2 ) +
S Rm
n +1 ∥ebh ∥l2 (H 1 ) ≤ Ch,
(5.10)
where C depends on Re, Rm and S. n+1/2
Proof. Taking vh = 2euh
∥
in (5.6) leads to
n+1/2 2 euh L2
n+1/2
∥ − ∥enuh ∥2L2 + ∥euh
=2
n+1/2
Eu
− enuh ∥2L2
1t
− Enu
1t
n+1/2
, euh
+
2 Re
n+1/2 2 L2
∥∇ euh
∥
n+1/2 n+1/2 n+1/2 + 2 b1 (unh , uh , euh ) − b1 (un , un+1/2 , euh ) n+1/2
+ 2S bnh × curl bnh+1 − bn × curl bn+1 , euh
= J1 + J2 + J3 .
(5.11)
From Hölder inequality and Young’s inequality, J1 satisfies n+1/2
J1 ≤ 2(1t )−1 (∥Enu+1/2 ∥L2 + ∥Enu ∥L2 )∥euh
≤
1 4Re
n+1/2 2 L2
∥∇ euh
∥L2
∥ + Ch4 (1t )−2 (∥un+1/2 ∥2H 2 + ∥un ∥2H 2 ).
182
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184 n+1/2
The definitions of enuh , euh n+1/2
2 b1 (unh , uh
and the trilinear term b1 imply that
n+1/2
, euh
n+1/2
) − b1 (un , un+1/2 , euh
)
n+1/2 n+1/2 n+1/2 = 2 b1 (Enu − enuh , Enu+1/2 , euh ) + b1 (Enu − enuh , un+1/2 , euh ) + b1 (un , Eun+1/2 , euh ) . By using (5.1), (5.2) and (4.20), we estimate the three terms in the above equation as follows: n+1/2
b1 (Enu − enuh , Enu+1/2 , euh
≤ C (∥∇ ∥ + ∥∇ Enu L2
1
≤
) n+1/2
∥ )∥∇ Enu+1/2 ∥L2 ∥∇ euh
enuh L2
∥L2
n+1/2 2 L2
∥ + Ch2 ∥un+1/2 ∥H 2 + C Reh2 ∥∇ enuh ∥2L2 ∥un+1/2 ∥2H 2 ,
∥∇ euh
6Re
and n+1/2
b1 (Enu − enuh , un+1/2 , euh
≤ C (∥ ∥ + ∥ Enu L2
1
≤
∥∇
6Re
) n+1/2
∥ )∥un+1/2 ∥H 2 ∥∇ euh
enuh L2
∥L2
n+1/2 2 euh L2
∥ + Ch2 ∥un+1/2 ∥H 2 + C Re∥enuh ∥2L2 ∥un+1/2 ∥2H 2 ,
and n+1/2
b1 (un , Eun+1/2 , euh
n+1/2
) ≤ C ∥∇ un ∥L2 ∥∇ Enu+1/2 ∥L2 ∥∇ euh ≤
1
∥L2
n+1/2 2 L2
∥ + Ch2 ∥un+1/2 ∥H 2 .
∥∇ euh
6Re
Therefor, J2 satisfies n+1/2
n+1/2
J2 = 2 b1 (unh , uh 1
≤
2Re
n+1/2
) − b1 (un , un+1/2 , euh
, euh
)
n+1/2 2 L2
∥ + Ch2 ∥un+1/2 ∥H 2 + C Re(h2 ∥∇ enuh ∥2L2 + ∥enuh ∥2L2 )∥un+1/2 ∥2H 2 .
∥∇ euh
By a similar method, we have n+1/2
J3 = 2S bnh × curl bnh+1 − bn × curl bn+1 , euh n+1/2
= −2S (bnh × curl enbh+1 , euh + 2S ((Enb − enbh ) × curl ≤ −2S ( +
bnh
1 4Re
× curl
enbh+1
,
n+1/2
) − 2S (enbh × curl Enb+1 , euh
n+1/2 bn+1 euh
,
n+1/2 euh
) + Ch∥
) n+1/2
) + 2S (Enb × curl Enb+1 , euh
enbh 2L2
∥ ∥
n+1/2
) + 2S (bn × curl Enb+1 , euh
)
bn+1 2H 2
∥
n+1/2 2 L2
∥ + Ch2 ∥bn+1 ∥H 2 + C ∥enbh ∥2L2 ∥bn+1 ∥2H 2 ,
∥∇ euh
where we use the inverse inequality ∥vh ∥L3 ≤ Ch−1/6 ∥vh ∥L2 . Combining these estimates into (5.11) yields n+1/2 2 L2
n+1/2
∥ − ∥enuh ∥2L2 + ∥euh
∥euh
1t
n+1/2
≤ −2S (bnh × curl enbh+1 , euh
− enuh ∥2L2
+
1 Re
n+1/2 2 L2
∥∇ euh
∥
) + Ch4 (1t )−2 (∥un+1/2 ∥2H 2 + ∥un ∥2H 2 ) + Ch2 (∥un+1/2 ∥H 2 + ∥bn+1 ∥2H 2 )
+ C Reh2 ∥∇ enuh ∥2L2 ∥un+1/2 ∥2H 2 + C Re∥enuh ∥2L2 ∥un+1/2 ∥2H 2 + C ∥enbh ∥2L2 ∥bn+1 ∥2H 2 .
(5.12)
n +1 Taking vh = 2euh in (5.8) and qh = enph+1 in (5.9) leads to n+1/2 2 L2
n+1 2 ∥euh ∥L2 − ∥euh
n+1/2 2 L2
∥ + ∥enuh+1 − euh
∥
1t
=
n+1/2
Eun+1 − Eu
1t
+
1 Re
n+1/2 2 L2
(∥∇ enuh+1 ∥2L2 − ∥∇ euh
,
enuh+1
− a1 (Enu+1/2 , enuh+1 )
≤ Ch4 (1t )−2 (∥un+1/2 ∥2H 2 + ∥un+1 ∥2H 2 ) +
1 2Re
n+1/2
∥ + ∥∇(enuh+1 − euh
∥∇ enuh+1 ∥2L2 + Ch2 ∥un+1/2 ∥2H 2 ,
)∥2L2 )
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
183
which together with (5.12) yields n +1 2 ∥euh ∥L2 − ∥enuh ∥2L2
1
n+1 2 + ∥∇ euh ∥L2 1t Re n + 1 / 2 ≤ −2S (bnh × curl enbh+1 , euh ) + Ch4 (1t )−2 (∥un+1/2 ∥2H 2 + ∥un ∥2H 2 + ∥un+1 ∥2H 2 )
+ Ch2 (∥un+1/2 ∥H 2 + ∥bn+1 ∥2H 2 ) + C Reh2 ∥∇ enuh ∥2L2 ∥un+1/2 ∥2H 2 + C Re∥enuh ∥2L2 ∥un+1/2 ∥2H 2 + C ∥enbh ∥2L2 ∥bn+1 ∥2H 2 .
(5.13)
Setting wh = 2Senbh+1 in (5.7) leads to n+1 2 S ∥ebh ∥L2 − S ∥enbh ∥2L2 + S ∥enbh+1 − enbh ∥2L2
2SC
∥en+1 ∥2 1 1t Rm bh H n+1/2 n +1 ≤ 2S (un+1/2 × enbh , curl enbh+1 ) − 2S (un+1/2 × Enb , curl enbh+1 ) + 2S (euh × bnh , curl ebh ) +
n +1 − 2S (Enu+1/2 × bn , curl enbh+1 ) + 2S (Enu+1/2 × enbh , curl enbh+1 ) − 2S (Enu+1/2 × Enb , curl ebh )
SC ∥en+1 ∥2 1 + Ch4 (1t )−2 (∥bn ∥2H 2 + ∥bn+1 ∥2H 2 ) 2Rm bh H SC n+1/2 n +1 ∥en+1 ∥2 1 + Ch4 (1t )−2 (∥bn ∥2H 2 + ∥bn+1 ∥2H 2 ) + 2S (euh × bnh , curl ebh ) + C ∥enbh ∥2L2 ∥un+1/2 ∥2H 2 ≤ Rm bh H
+
+ C ∥Eun+1/2 ∥2L2 ∥bn ∥2H 2 + C ∥un+1/2 ∥2H 2 ∥Enb ∥2L2 + C ∥enbh ∥2L3 ∥∇ un+1/2 ∥2L2 + C ∥∇ Enu+1/2 ∥2L2 ∥Enb ∥L2 ∥Enb ∥H 1 .
(5.14)
Under the condition 1t ≥ Ch, adding (5.13) and (5.14), and taking the sum from n = 0 to n = m ≤ N − 1, we get
∥
+1 2 em uh L2
∥ + S∥
∥ + 1t
+1 2 em bh L2
m 1 n =0
≤ Ch2 1t
m
Re
∥∇
enuh+1 2L2
∥ +
SC Rm
∥
enbh+1 2H 1
+ C 1t
∥
m ∥un+1/2 ∥2H 2 + ∥un+1 ∥2H 2 + ∥bn+1 ∥2H 2 + C Reh2 1t ∥∇ enuh ∥2L2 ∥un+1/2 ∥2H 2
n =0 m
n =0
Re∥un+1/2 ∥2H 2 ∥enuh ∥2L2 + ∥un+1/2 ∥2H 2 ∥enbh ∥2L2 + ∥bn+1 ∥2H 2 ∥enbh ∥2L2
n =0
≤ Ch2 + C8 Reh1t
m
∥∇ enuh ∥2L2 + C Re1t
n =0
+ C 1t
m
m n+1/2 2 ∥u ∥H 2 ∥enuh ∥2L2 n=0
∥un+1/2 ∥2H 2 ∥enbh ∥2L2 + ∥bn+1 ∥2H 2 ∥enbh ∥2L2 .
n =0
Let h be sufficiently small such that C8 h < complete the proof of Lemma 5.1.
1 . Re2
Applying the discrete Gronwall’s inequality to the above inequality, we
Combining the temporal error estimates (3.12)–(3.16) with the results in Lemma 5.1, we obtain the temporal–spatial error estimates (3.17) and (3.18). 6. Conclusions In this paper we propose a fractional-step scheme for solving the stationary MHD problems in three-dimensional bounded domains, which is an extension of the previous work for the Navier–Stokes problems studied in [24]. The main feature is that the nonlinearities for velocity and the incompressibility in (1.1)–(1.2) are split into two steps. Compared to the projection scheme [14], the fractional-step scheme studied in this paper allows to enforce the original boundary conditions of the problems in all substeps of the fractional-step scheme, and allows for a discrete energy estimate. We have proved that the time-discrete fractional-step √ scheme provides the temporal error estimates of O (1t ) for the velocity and the magnetic field in l2 (H1 ) ∩ l∞ (L2 ) and O ( 1t ) for the pressure in l2 (L2 ). For the fully discrete FEM scheme, the temporal–spatial error estimates of O (1t + h) for the velocity and the magnetic field in l2 (H1 ) ∩ l∞ (L2 ) are obtained under √the constraint 1t ≥ Ch. However, for the projection scheme, the temporal–spatial error orders derived in [14] are only O ( 1t + h) for the velocity and the magnetic field in l2 (H1 ). Although some temporal–spatial convergence orders are derived in this paper, there are still two questions worthy of study. One is whether the spatial error estimate for the velocity and the magnetic field in l∞ (L2 ) can be improved to
184
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
be optimal. Another is whether there is an optimal error estimate for the pressure in l2 (L2 ). We observed that GuillénGonzález and Redondo-Neble have studied optimal error estimates of the fractional-step scheme for the Navier–Stokes problems in [29,30]. But it is much more complicated to study these optimal error estimates for the MHD problem due to the appearance of the magnetic field equations. In addition, the extension of this fractional-step scheme to higher-order temporal accuracy also an interesting problem, could possibly be achieved using a Crank–Nicolson scheme. However, there needs much higher regularity of the solution to the MHD problems. These objectives are worth studying in our future work.
Acknowledgments The authors would like to thank the anonymous reviewers for their careful reviews and valuable comments. These comments are all valuable and very helpful for improving this manuscript. This manuscript is supported by Zhejiang Provincial Natural Science Foundation with Grant Nos. LY16A010017 and LY14A010020. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30]
W. Hughes, F. Young, The Electromagnetics of Fluids, Wiley, New York, 1966. R. Moreau, Magneto-hydrodynamics, Kluwer Academic Publishers, 1990. J. Gerbeau, C. Le Bris, T. Lelièvre, Mathematical Methods for the Magnetohydrodynamics of Liquid Metals, Oxford University Press, Oxford, 2006. M. Gunzburger, A. Meir, J. Peterson, On the existence, uniqueness, and finite element approximation of solutions of the equations of stationary, incompressible magnetohydrodynamics, Math. Comp. 56 (1991) 523–563. J. Douglas Jr., J. Wang, An absolutely stabilized finite element method for the Stokes problem, Math. Comp. 52 (1989) 495–508. J. Gerbeau, A stabilized finite elemenet method for the incompressible magnetohydrodynamic equations, Numer. Math. 87 (2000) 83–111. Y. He, Unconditional convergence of the Euler semi-implicit scheme for the three-dimensional incompressible MHD equations, IMA J. Numer. Anal. 35 (2015) 767–801. D. Schötzau, Mixed finite element methods for stationary incompressible magneto-hydrodynamics, Numer. Math. 96 (2004) 771–800. S. Badia, R. Codina, R. Planas, On an unconditionally convergent stabilized finite element approximation of resistive magnetohydrodynamics, J. Comput. Phys. 234 (2013) 399–416. S. Badia, R. Planas, J. Gutiérrez-Santacreu, Unconditionally stable operator splitting algorithms for the incompressible magnetohydrodynamics system discretized by a stabilized finite element formulation based on projections, Internat. J. Numer. Methods Engrg. 93 (2013) 302–328. M. Belenli, S. Kaya, L. Rebholz, N. Wilson, A subgrid stabilization finite element method for incompressible magnetohydrodynamics, Int. J. Comput. Math. 90 (2013) 1506–1523. R. Codina, N. Hernández, Approximation of the thermally coupled MHD problem using a stabilized finite element method, J. Comput. Phys. 230 (2011) 1281–1303. 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. A. Prohl, Convergent finite element discretizations of the nonstationary incompressible magnetohydrodynamic system, ESAIM:M2AN 42 (2008) 1065–1087. D. Sondak, J. Shadid, A. Oberai, R. Pawlowski, E. Cyr, T. Smith, A new class of finite element variational multiscale turbulence models for incompressible magnetohydrodynamics, J. Comput. Phys. 295 (2015) 596–616. G. Yuksel, R. Ingram, Numerical analysis of a finite element Crank–Nicolson discretization for MHD flows at small magnetic Reynolds numbers, Int. J. Numer. Anal. Model. 10 (2013) 74–98. G. Yuksel, O. Isik, Numerical analysis of Backward-Euler discretization for simplified magnetohydrodynamic flow, Appl. Math. Model. 39 (2015) 1889–1898. Y. Zhang, Y. Hou, L. Shan, Numerical analysis of the Crank–Nicolson extrapolation time discrete scheme for magnetohydrodynamics flows, Numer. Methods Partial Differential Equations 31 (2015) 2169–2208. A. Chorin, Numerical solution of the Navier–Stokes equations, Math. Comp. 22 (1968) 745–762. R. Temam, Sur l’approximation de la solution des equations de Navier–Stokes par la méthode des pas fractionnaires II, Arch. Ration. Mech. Anal. 33 (1969) 377–385. A. Prohl, On pressure approximation via projection methods for nonstationary incompressible Navier–Stokes equations, SIAM J. Numer. Anal. 47 (2008) 158–180. R. Temam, Remark on the pressure boundary condition for the projection method, Theor. Comput. Fluid Dyn. 3 (1991) 181–184. J. Blasco, R. Codina, A. Huerta, A fractional-step method for the incompressible Navier–Stokes equations related to a predictor-multicorrector algorithm, Internat. J. Numer. Methods Fluids 28 (1997) 1391–1419. J. Blasco, R. Codina, Error estimates for an operator-splitting method for incompressible flows, Appl. Numer. Math. 51 (2004) 1–17. R. Temam, Navier–Stokes Equations, North-Holland Publishing Company, Amsterdam, 1977. R. Adams, Sobolev Spaces, Academic Press, New York, 1975. J. Heywood, R. Rannacher, Finite-element approximation of the nonstationary Navier–Stokes problem Part IV: error analysis for second-order time discretization, SIAM J. Numer. Anal. 27 (1990) 353–384. D. Gilbarg, N.S. Trudinger, Elliptic Partial Differential Equations of Second Order, Springer-Verlag, 1983. F. Guillén-González, M.V. Redondo-Neble, New error estiamtes for a viscosity-splitting scheme in time for the three-dimensional Navier–Stokes equations, IMA J. Numer. Anal. 31 (2011) 556–579. F. Guillén-González, M.V. Redondo-Neble, Spatial error estimates for a finite element viscosity-splitting scheme for the Navier–Stokes equations, Int. J. Numer. Anal. Model. 10 (2013) 826–844.
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Error analysis of a fractional-step method for magnetohydrodynamics equations Rong An ∗ , Can Zhou College of Mathematics and Information Science, Wenzhou University, Wenzhou 325035, PR China
article
info
Article history: Received 11 November 2015 Received in revised form 21 June 2016 MSC: 65M12 76W05
abstract This paper focuses on a fractional-step finite element method for the magnetohydrodynamics problems in three-dimensional bounded domains. It is shown that the proposed fractional-step scheme allows for a discrete energy identity. A rigorous error analysis is presented. We derive the temporal and spatial error estimates of O (∆t + h) for the velocity and the magnetic field in the discrete space l2 (H1 ) ∩ l∞ (L2 ) under the constraint ∆t ≥ Ch. © 2016 Elsevier B.V. All rights reserved.
Keywords: Magnetohydrodynamics equations Fractional-step method Stability Temporal errors Spatial errors
1. Introduction The incompressible magnetohydrodynamics (MHD) problems are used to describe the flow of a viscous, incompressible and electrically conducting fluid, which are governed by the following time-dependent nonlinear coupled problems:
∂u 1 − 1u + (u · ∇)u + ∇ p + Sb × curl b = f, ∂t Re div u = 0, 1 ∂b + curl (curl b) − curl (u × b) = 0, ∂t Rm div b = 0,
(1.1) (1.2) (1.3) (1.4)
for x ∈ Ω and t ∈ (0, T ) with T > 0, where Ω ⊂ R3 is a bounded and simply-connected domain which is either convex or has a C 1,1 boundary ∂ Ω . Re, Rm and S are three positive constants and denote the Reynolds number, the magnetic Reynolds number and the coupling number, respectively. The vector-value function f represents the body forces applied to the fluid. The MHD problems (1.1)–(1.4) couple the incompressible Navier–Stokes equations with Maxwell’s equations. Thus, the unknowns in (1.1)–(1.4) are the fluid velocity u, the pressure p and the magnetic field b. We refer to Hughes [1] and Moreau [2] for the understanding of the physical background of the MHD problems. To study (1.1)–(1.4), the appropriate initial and boundary conditions are needed. For the sake of simplicity, in this paper, we consider the following initial and
∗
Corresponding author. E-mail addresses: [email protected], [email protected] (R. An).
http://dx.doi.org/10.1016/j.cam.2016.09.005 0377-0427/© 2016 Elsevier B.V. All rights reserved.
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
169
boundary conditions: u(x, 0) = u0 , u = 0,
b · n = 0,
b(x, 0) = b0 curl b × n = 0
in Ω ,
(1.5)
on ∂ Ω × [0, T ],
(1.6)
where n denotes the unit outward normal vector on ∂ Ω . It is necessary to require that u0 and b0 satisfy the compatibility conditions div u0 = 0 and div b0 = 0. The numerical methods for the incompressible MHD problems have received much attention in the last decades. We refer to Gerbeau–Bris–Lelièvre [3] for a review of numerical methods. The mixed finite element approximation was first proposed and studied for the stationary MHD problems in [4], where H1 -conforming elements were used to discretize the magnetic field provided that Ω is either convex or has a C 1,1 boundary ∂ Ω . Inspired by the stabilization method for Stokes problems in [5], a stabilization mixed finite element method for the stationary MHD problems was developed by Gerbeau [6]. For the time-dependent MHD problems (1.1)–(1.6), He proposed a linearized semi-implicit Euler scheme in [7], where L2 -unconditional convergence was proved by using the negative norm technique. For the non-convex domain or Lipschitz polyhedra domain of engineering practice, the magnetic field b may have regularity below H1 (Ω ). In this case, the H1 -conforming finite element discretization for b, albeit stable, may not converge to the correct magnetic field. A mixed finite element formulation based on H(curl)-elements (or Nédélec elements) for b was proposed and studied by Schötzau in [8] for the stationary MHD problems. Other different numerical methods can be found in [9–18] and references cited therein. Roughly speaking, the difficulties in solving the MHD problems numerically are mainly of three kinds: the mixed type of the equations; the incompressible condition and the treatment of the pressure; the nonlinearities of the problems, which are very similar to the incompressible Navier–Stokes problems. In the 1960s, Chorin [19] and Temam [20] proposed a projection method for Navier–Stokes problems, which decoupled the velocity and the pressure in the Navier–Stokes problems. The idea of the projection method is first to compute a velocity field without taking into account incompressibility, and then perform a pressure correction, which is a projection back to the subspace of solenoidal (divergence-free) vector fields. However, the drawback is the appearance of the numerical boundary layer due to the incompatibility of the pressure boundary conditions with those of the original Navier–Stokes problems [21,22]. For the time-dependent MHD problems (1.1)–(1.6), inspired by the projection method in [19,20], Prohl in [14] proposed a projection scheme. However, the projection scheme in [14] does not allow for a discrete energy estimate. To avoid using artificial boundary conditions of pressure type, some fractionalstep schemes for the Navier–Stokes problems were introduced and studied in [23,24]. It is a two-step scheme in which the incompressibility and the nonlinearities of the Navier–Stokes problems are split into different steps, and allows the enforcement of the original boundary conditions in all substeps. In this paper, we propose a two-step fractional-step scheme to approximate the solution of the MHD problems (1.1)–(1.4) with the initial and boundary conditions (1.5)–(1.6). We will prove that the proposed fractional-step scheme allows for a discrete energy identity. To state the main results derived, we introduce the following notations. Let X be a Banach space equipped with norm ∥ · ∥X . Let 0 = t0 < t1 < · · · < tN = T be a uniform partition of the time interval [0, T ] with time step 1t = T /N and tn = n1t for 0 ≤ n ≤ N. We denote two discrete norms by
∥u ∥l2 (X ) = 1t n
N
1/2 un 2X
∥ ∥
,
n=1
∥un ∥l∞ (X ) = max ∥un ∥X . 1≤n≤N
It is proved that the time-discrete fractional-step √ scheme provides the temporal error estimates of O (1t ) for the velocity and the magnetic field in l2 (H1 )∩ l∞ (L2 ) and O ( 1t ) for the pressure in l2 (L2 ). For the fully discrete fractional-step scheme, the finite element error estimates of O (1t + h) for the velocity and the magnetic field in l2 (H1 ) ∩ l∞ (L2 ) are obtained under the constraint 1t ≥ Ch. The remainder of this paper is organized as follows: in the next section, we begin with some notations, lay out some assumptions and recall some known inequalities frequently used. The new linearized projection scheme and the main results are presented in Section 3. Meanwhile, the discrete energy identity is derived in Section 3. The proof containing the main results is given in Sections 4 and 5, which is split into several theorems. 2. Mathematical setting For the mathematical setting of the MHD problems (1.1)–(1.4) with the initial and boundary conditions (1.5)–(1.6), we introduce some function spaces and their associated norms. For k ∈ N+ , 1 ≤ s ≤ +∞, let W k,s (Ω ) denote the standard Sobolev space. The norm in W k,s (Ω ) is defined by
∥u∥W k,s =
1/s β s |∇ u| dx
for 1 ≤ s < +∞,
sup |∇ β u|
for s = +∞,
|β|≤k Ω
|β|≤k Ω
170
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
where ∇ is the differential operator with respect to x and
∇β =
∂ |β| β
β
β
∂ x1 1 ∂ x2 2 ∂ x3 3 k,s
for the multi-index β = (β1 , β2 , β3 ), β1 , β2 , β3 ≥ 0 and |β| = β1 + β2 + β3 . We define W0 (Ω ) to be the subspace of W k,s (Ω ) of functions with zero trace on ∂ Ω . When s = 2, we simply use H k (Ω ) to denote W k,2 (Ω ). The dual space of H01 (Ω ) is denoted by H −1 (Ω ). The boldface Sobolev spaces Hk (Ω ), Wk,s (Ω ) and Ls (Ω ) are used to denote the vector Sobolev spaces H k (Ω )3 , W k,s (Ω )3 and Ls (Ω )3 , respectively. For any u, v ∈ L2 (Ω ), (u, v) = Ω u · vdx denotes L2 (Ω ) inner product. Let X be a Banach space. Its dual space is denoted by X ′ . For some T > 0 or T = ∞, denote Bochner spaces Ls (0, T ; X ), 1 ≤ s < +∞, the spaces of measurable functions from the interval [0, T ] into X such that T
∥v(t )∥sX dt < +∞.
0
If s = +∞, the functions in L∞ (0, T ; X ) are required to satisfy sup ∥v(t )∥X < +∞.
t ∈[0,T ]
The symbols C , C1 , C2 , . . . are used to denote the generic positive constants independent of the time step 1t and the mesh size h. Let us denote H = {u ∈ L2 (Ω ), div u = 0 in Ω , u · n = 0 on ∂ Ω }, V0 = {u ∈ V, div u = 0 in Ω },
V = H10 (Ω ),
X = {u ∈ H (Ω ), u · n = 0 on ∂ Ω },
X0 = {u ∈ X, div u = 0 in Ω },
1
M = L20 (Ω ) =
q ∈ L2 (Ω ),
Ω
qdx = 0 .
Thanks to the Poincaré inequality, the norm in V can be defined by
∥u∥V =
Ω
1/2 |∇ u|2 dx .
Introduce two bilinear forms
a2 (u, v) =
Re
1 Rm
1
a1 (u, v) =
Ω
Ω
∇ u : ∇ vdx ∀ u, v ∈ V,
curl u · curl vdx ∀ u, v ∈ X,
d(v, p) =
pdiv vdx ∀ v ∈ V, p ∈ M
Ω
and a trilinear form b(u, v, w) =
Ω
(u · ∇)v · wdx ∀ u, v, w ∈ H1 (Ω ).
It is easy to verify that if div u = 0 in Ω and u · n = 0 on ∂ Ω , the trilinear form b(u, v, w) satisfies the following important property: b(u, v, w) = −b(u, w, v).
(2.1)
With the above notations, given u0 , b0 ∈ H and f ∈ L (0, T ; L (Ω )), the weak variational formulation of the MHD problems (1.1)–(1.6) is defined to find (u, p, b) ∈ L2 (0, T ; V) × L2 (0, T ; M ) × L2 (0, T ; X) such that ∞
2
(ut , v) + a1 (u, v) + b(u, u, v) + S (b × curl b, v) − d(v, p) = (f, v) ∀ v ∈ V,
(2.2)
and
(bt , w) + a2 (b, w) +
1 Rm
(div b, div w) − (u × b, curl w) = 0 ∀ w ∈ X,
(2.3)
where we use the formula:
curl u · vdx = Ω
u · curl vdx +
Ω
∂Ω
n × u · vds.
(2.4)
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
171
Now, we recall Stokes operator A [25]. Let us denote the orthogonal projection operator by PH from L2 (Ω ) onto H. Then Stokes operator A is defined by
∀ u ∈ V0 ∩ H2 (Ω ) := D(A).
Au = −PH 1u
The emphasis of this paper is to prove some temporal and spatial error estimates of the time-discrete scheme (3.1)–(3.3) and the fully-discrete scheme (3.8)–(3.11) stated in Section 3. Then we make the following assumptions. Assumption 1. The initial data u0 , b0 and f satisfy
∥Au0 ∥L2 + ∥b0 ∥H 2 + sup ∥f(t )∥L2 + ∥ft (t )∥L2 ≤ C3 . 0 ≤t ≤T
Assumption 2. The MHD problems (1.1)–(1.6) admit a unique local strong solution (u, p, b) on [0, T ] such that sup 0≤t ≤T
∥Au(t )∥L2 + ∥b(t )∥H 2 + ∥p(t )∥H 1 + ∥ut (t )∥L2 + ∥bt (t )∥L2 T
+ 0
∥∇ ut ∥2L2 + ∥bt ∥2H 1 + ∥utt ∥2H −1 + ∥btt ∥2H −1 + τ (t )∥utt ∥2L2 + τ (t )∥btt ∥2L2 dt ≤ C4 ,
where τ (t ) = min{t , 1} is a weight function vanishing at t = 0. Assumption 3. For given g ∈ L2 (Ω ), the Stokes problems
−1v + ∇π = g,
div v = 0 v=0
in Ω , on ∂ Ω
admit a unique solution (v, π ) ∈ D(A) × H 1 (Ω ) ∩ M satisfying
∥Av∥L2 + ∥∇π ∥L2 ≤ C ∥g∥L2 . Assumption 4. For given g ∈ L2 (Ω ), the magnetic problems curl curl v = g, v · n = 0,
div v = 0,
in Ω ,
curl v × n = 0,
on ∂ Ω
admit a unique solution v ∈ X0 ∩ H2 (Ω ) satisfying ∥v∥H 2 ≤ C ∥g∥L2 . The following known inequalities are frequently used [26]: 1/2
1/2
∀ v ∈ H1 (Ω ),
(2.5)
1/2
1/2
∀ v ∈ H (Ω ),
(2.6)
∥v∥L3 ≤ C ∥v∥L2 ∥v∥H 1 ∥v∥L∞ ≤ C ∥v∥H 1 ∥v∥H 2
2
∥v∥H 2 ≤ C ∥Av∥L2 ∀ v ∈ D(A),
(2.7)
∥v∥H 2 ≤ C ∥curl curl v∥L2 ∀ v ∈ X0 ∩ H (Ω ), ∥v∥H 1 ≤ C (∥curl v∥L2 + ∥div v∥L2 ) ∀ v ∈ X. 2
(2.8) (2.9)
Finally, we recall a discrete version of Gronwall’s inequality established in [27]. Lemma 2.1. Let ak , bk , ck and γk , for integers k ≥ 0, be the nonnegative numbers such that an + 1 t
n k=0
bk ≤ 1 t
n
γk ak + 1t
k=0
n
ck + B
for n ≥ 0.
(2.10)
k=0
Suppose that 1t γk < 1, for all k, and set σk = (1 − 1t γk )−1 . Then an + 1 t
n k=0
bk ≤ exp 1t
n k=0
γk σk
1t
n
ck + B
for n ≥ 0.
(2.11)
k=0
Remark 2.1. If the first sum on the right in (2.10) extends only up to n − 1, then the estimate (2.11) holds for all k > 0 with σk = 1.
172
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
3. Fractional-step scheme Let 0 = t0 < t1 < · · · < tN = T be a uniform partition of the time interval [0, T ] with time step 1t = T /N and tn = n1t , 0 ≤ n ≤ N, where [0, T ] is the maximal time interval such that a local strong solution exists and satisfies Assumption 2. For any sequence of functions {g n }Nn=0 , we define Dt g n+1 = (g n+1 − g n )/1t for 0 ≤ n ≤ N − 1. Start with u0 = u0 and b0 = b0 . The proposed time discrete fractional-step scheme is the following two-step scheme: Step I: The first step of the fractional-step scheme consists of finding an intermediate velocity un+1/2 and bn+1 such that un+1/2 − un
1t bn+1 − bn
1t
+
−
1 Rm
1 Re
1un+1/2 + (un · ∇)un+1/2 + Sbn × curl bn+1 = f(tn+1 ),
curl (curl bn+1 ) − curl (un+1/2 × bn ) = 0,
div bn+1 = 0
(3.1) (3.2)
with the boundary conditions un+1/2 = 0, bn+1 · n = 0 and curl bn+1 × n = 0 on ∂ Ω . Step II: Given un+1/2 from (3.1), the second step of the fractional-step scheme consists of finding un+1 and pn+1 such that un+1 − un+1/2
1t
−
1 Re
∆(un+1 − un+1/2 ) + ∇ pn+1 = 0,
div un+1 = 0
(3.3)
with the boundary condition un+1 = 0 on ∂ Ω . Since the problems (3.1)–(3.3) consist of three linear problems, their solvabilities should follow easily if we can establish the a priori energy estimates for the solutions to these problems. The following energy identity serves this purpose. Theorem 3.1. Let un+1/2 ∈ V, bn+1 ∈ X0 and (un+1 , pn+1 ) ∈ V0 × M be the weak solutions to the problems (3.1)–(3.3). For 0 ≤ m ≤ N − 1, the following energy identity holds:
∥um+1 ∥2L2 + S ∥bm+1 ∥2L2 +
m 2S 1t
Rm
∥curl bn+1 ∥2L2
n =0
m + (∥un+1 − un+1/2 ∥2L2 + ∥un − un+1/2 ∥2L2 + S ∥bn+1 − bn ∥2L2 ) n=0 m 1t (∥∇ un+1 ∥2L2 + ∥∇ un+1/2 ∥2L2 + ∥∇(un+1 − un+1/2 )∥L2 ) +
Re n=0
= ∥u0 ∥2L2 + S ∥b0 ∥2L2 + 21t
m (f(tn+1 ), un+1/2 ).
(3.4)
n =0
Proof. Multiplying (3.1)–(3.3) by 21tun+1/2 , 2S 1tbn+1 and 21tun+1 , respectively, and adding three equations and using the following formula 2(u − v, u) = ∥u∥2L2 + ∥u − v∥2L2 − ∥v∥2L2
∀ u, v ∈ L2 (Ω ),
(3.5)
and the vector relation
(u × curl v, w) = (w × u, curl v) ∀ u, v, w ∈ H1 (Ω ), we complete the proof of (3.4).
(3.6)
As a direct consequence of (2.9), (3.4) and div bn+1 = 0, we can prove that un+1/2 , un+1 and bn+1 are uniformly bounded in l (H1 ), l∞ (L2 ) ∩ l2 (H1 ) and l∞ (L2 ) ∩ l2 (H1 ), respectively. 2
Corollary 3.1. The fractional-step scheme (3.1)–(3.3) is unconditionally stable. Moreover, there exists some C > 0 such that
∥un+1 ∥l∞ (L2 ) + S ∥bn+1 ∥l∞ (L2 ) + +
1 Re
∥∇ un+1/2 ∥l2 (L2 ) +
S Rm
1 Re
∥∇ un+1 ∥l2 (L2 )
∥bn+1 ∥l2 (H 1 ) ≤ C ,
where C depends on Re, Rm, S, u0 , b0 and f. Next, we give the finite element approximations of (3.1)–(3.3). Let Th be a quasi-uniform partition of Ω into tetrahedra of diameters by h with 0 < h < 1. Let C (Ω ) denote the space of continuous function in Ω = Ω ∪ ∂ Ω . For every K ∈ Th and a nonnegative integer r, Pr (K ) denotes the space of the polynomials on K of degree at most r. We use the boldface spaces
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
173
C(Ω ) and Pr (K ) to denote the vector spaces C (Ω )3 and Pr (K )3 . Define the following conforming finite element subspaces of V, M and X for the velocity, the pressure and the magnetic field, respectively, by Vh = {vh ∈ C(Ω ), vh ∈ P2 (K ), ∀ K ∈ Th } ∩ V, Vh = {vh ∈ C(Ω ), vh ∈ P1 (K ), ∀ K ∈ Th } ∩ V, Xh = {vh ∈ C(Ω ), vh ∈ P1 (K ), ∀ K ∈ Th } ∩ X, Mh = {qh ∈ C (Ω ), qh ∈ P1 (K ), ∀ K ∈ Th } ∩ M . In such a way, the approximation spaces Vh and Mh satisfy the discrete inf–sup condition, i.e., there exists some positive constant β > 0 such that [25]
β∥qh ∥L2 ≤
d(vh , qh )
sup
vh ∈Vh /{0}
∥∇ vh ∥L2
.
(3.7)
Under the above notations, the finite element approximations of the fractional-step scheme (3.1)–(3.3) are described as follows. Start with u0h = Rh u0 , b0h = Qh b0 , where the projectors Rh and Qh are defined in Section 5. For 0 ≤ n ≤ N − 1, n+1/2
∈ Vh and bnh+1 ∈ Xh such that
given unh ∈ Vh and bnh ∈ Xh , find uh
n+1/2
uh
− unh
1t
n+1/2
, vh + a1 (uh
n+1/2
, vh ) + b1 (unh , uh
, vh ) + S (bnh × curl bnh+1 , vh ) = (f(tn+1 ), vh ) ∀ vh ∈ Vh , (3.8)
where b1 (u, v, w) =
1 2
b(u, v, w) −
1 2
b(u, w, v) ∀ u, v, w ∈ V,
and
bnh+1 − bnh
1t
, wh + a2 (bnh+1 , wh ) +
1 Rm
n+1/2
(div bnh+1 , div wh ) − (uh
× bnh , curl wh ) = 0 ∀ wh ∈ Xh .
(3.9)
n+1/2
Then given uh
∈ Vh , find (unh+1 , pnh+1 ) ∈ Vh × Mh such that n+1/2
uhn+1 − uh
1t
n+1/2
, vh + a1 (unh+1 − uh
, vh ) − d(vh , phn+1 ) = 0 ∀ vh ∈ Vh , d(uhn+1 , qh ) = 0
(3.10)
∀ qh ∈ Mh .
By using a similar argument for Theorem 3.1, we can prove that the finite element approximation solutions , phn+1 satisfy the following discrete energy estimate.
bnh+1
(3.11) n+1/2 uh
, uhn+1 ,
n+1/2
Theorem 3.2. For 0 ≤ n ≤ N − 1, let uh ∈ Vh , bnh+1 ∈ Xh and (uhn+1 , pnh+1 ) ∈ Vh × Mh be the solutions to the problems (3.8)–(3.10). Then there exists some C > 0 such that
∥unh+1 ∥l∞ (L2 ) + S ∥bnh+1 ∥l∞ (L2 ) +
1 Re
∥∇ unh+1 ∥l2 (L2 ) +
1 Re
n+1/2
∥∇ uh
∥l2 (L2 ) +
S Rm
∥bnh+1 ∥l2 (H 1 ) ≤ C ,
where C depends on Re, Rm, S, u0 , b0 and f. The emphasis of this paper is to show the following error estimates. Theorem 3.3. Suppose that Assumptions 1–4 are satisfied. Then there exists some τ0 > 0 such that when τ < τ0 , the following temporal error estimates hold:
∥u(tn+1 ) − un+1 ∥l∞ (L2 ) +
1 Re
S ∥b(tn+1 ) − bn+1 ∥l∞ (L2 ) +
∥∇ u(tn+1 ) − ∇ un+1 ∥l2 (L2 ) ≤ C 1t , S
∥b(tn+1 ) − bn+1 ∥l2 (H 1 ) ≤ C 1t , √ ∥p(tn+1 ) − pn+1 ∥l2 (L2 ) ≤ C 1t , √ ∥ τ (tn+1 )(∇ u(tn+1 ) − ∇ un+1 )∥l∞ (L2 ) ≤ C 1t , √ ∥ τ (tn+1 )(b(tn+1 ) − bn+1 )∥l∞ (H 1 ) ≤ C 1t , Rm
(3.12) (3.13) (3.14) (3.15) (3.16)
where C depends on Re, Rm and S, and un+1 , pn+1 , bn+1 are the numerical solutions to the problems (3.1)–(3.3) at time-step n + 1, and u(tn+1 ), p(tn+1 ), b(tn+1 ) are the exact solutions evaluated at time tn+1 , and τ (tn+1 ) = min{1, tn+1 } for 0 ≤ n ≤ N − 1.
174
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
Theorem 3.4. Suppose that Assumptions 1–4 are satisfied. If the time step 1t satisfies 1t ≥ Ch for some C > 0, then there exist some τ0 > 0 and h0 > 0 such that when τ < τ0 and h < h0 , the following error estimates hold:
∥u(tn+1 ) − unh+1 ∥l∞ (L2 ) + S ∥b(tn+1 ) − bnh+1 ∥l∞ (L2 ) ≤ C (1t + h), 1 Re
∥∇(u(tn+1 ) − unh+1 )∥l2 (L2 ) +
S Rm
∥b(tn+1 ) − bnh+1 ∥l2 (H 1 ) ≤ C (1t + h),
(3.17) (3.18)
where C depends on Re, Rm and S, and unh+1 , pnh+1 , bnh+1 are the finite element approximation solutions at time-step n + 1 for 0 ≤ n ≤ N − 1. The proof of Theorems 3.3 and 3.4 is split into several lemmas in Sections 4 and 5. 4. Temporal error analysis In this section, we will prove Theorem 3.3. For 0 ≤ n ≤ N − 1, let us denote errors by enu+1/2 = u(tn+1 ) − un+1/2 ,
eun+1 = u(tn+1 ) − un+1 , ebn+1
= b(tn+1 ) − b
n +1
,
enp+1 = p(tn+1 ) − pn+1 .
Taking t = tn+1 in (1.1)–(1.4) yields u(tn+1 ) − u(tn )
−
1t
1 Re
1u(tn+1 ) + (u(tn+1 ) · ∇)u(tn+1 )
+ ∇ p(tn+1 ) + Sb(tn+1 ) × curl b(tn+1 ) = f(tn+1 ) + Rnu+1 , div u(tn+1 ) = 0, b(tn+1 ) − b(tn )
1t
(4.1)
div b(tn+1 ) = 0,
+
1 Rm
(4.2)
curl (curl b(tn+1 )) − curl (u(tn+1 ) × b(tn+1 )) = Rnb+1 ,
(4.3)
where Rnu+1 = −
tn+1
1
1t
(t − tn )utt (t )dt ,
Rnb+1 = −
tn
1
1t
tn+1
(t − tn )btt (t )dt .
tn
From Assumption 2, one has
∥Rnu+1 ∥l2 (H −1 ) + ∥Rnb+1 ∥l2 (H −1 ) ≤ C 1t , ∥ τ (tn+1 )Run+1 ∥l2 (L2 ) + ∥ τ (tn+1 )Rnb+1 ∥l2 (L2 ) ≤ C 1t .
(4.4) (4.5)
Subtracting (3.1) and (3.2) from (4.1) and (4.3), respectively, we obtain n+1/2
eu
− enu
1t
−
1 Re
1enu+1/2 + ((u(tn+1 ) − u(tn )) · ∇)u(tn+1 ) + (enu · ∇)u(tn+1 ) + (un · ∇)eun+1/2
+ S (b(tn+1 ) − b(tn )) × curl b(tn+1 ) + Senb × curl b(tn+1 ) + Sbn × curl ebn+1 + ∇ p(tn+1 ) = Run+1
(4.6)
and ebn+1 − enb
1t
+
1 Rm
curl curl enb+1 − curl (enu+1/2 × bn )
− curl (u(tn+1 ) × enb ) − curl (u(tn+1 ) × (b(tn+1 ) − b(tn ))) = Rbn+1 .
(4.7)
√ First, we show that the fractional-step scheme (3.1)–(3.3) provides the temporal error estimates of O ( 1t ) for the velocity and the magnetic field in l2 (H1 ) ∩ l∞ (L2 ). Lemma 4.1. Suppose that Assumptions 1–2 are satisfied. Then there exists some positive constant C > 0 such that for every 0 ≤ m ≤ N − 1, +1 2 +1 2 ∥em ∥L2 + S ∥em ∥L2 + u b
m S 1t
Rm
n =0
∥curl enb+1 ∥2L2
+
m (∥enu+1/2 − enu ∥2L2 + ∥eun+1/2 − eun+1 ∥2L2 + S ∥ebn+1 − enb ∥2L2 ) n =0
m 1t + (∥∇ enu+1 ∥2L2 + ∥∇ enu+1/2 ∥2L2 + ∥∇(enu+1 − eun+1/2 )∥2L2 ) ≤ C 1t ,
Re n=0
where C depends on Re, Rm and S.
(4.8)
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184 n+1/2
Proof. Testing (4.6) and (4.7) by 21teu to
175
and 2S 1tenb+1 , respectively, using (3.6) and adding two resulting equations lead
∥eun+1/2 ∥2L2 + ∥enu+1/2 − enu ∥2L2 − ∥enu ∥2L2 + S ∥enb+1 ∥2L2 21t 2S 1t + S ∥enb+1 − enb ∥2L2 − S ∥enb ∥2L2 + ∥∇ enu+1/2 ∥2L2 + ∥curl ebn+1 ∥2L2 Re
≤ 21t |b(u(tn+1 ) − u(tn ), u(tn+1 ),
enu+1/2
Rm
) + ( , u(tn+1 ), eun+1/2 )| b enu enb
+ 2S 1t |((b(tn+1 ) − b(tn )) × curl b(tn+1 ) +
× curl b(tn+1 ), eun+1/2 )|
+ 2S 1t |(u(tn+1 ) × enb + u(tn+1 ) × (b(tn+1 ) − b(tn )), curl ebn+1 )| + 21t |(∇ p(tn+1 ), enu+1/2 )| + 21t |(Rnu+1 , enu+1/2 ) + S (Rbn+1 , enb+1 )| = I1 + · · · + I5 .
(4.9)
By using Assumption 2 and Young’s inequality, we have I1 = 21t |b(u(tn+1 ) − u(tn ), u(tn+1 ), enu+1/2 ) + b(enu , u(tn+1 ), enu+1/2 )|
≤ C 1t (∥u(tn+1 ) − u(tn )∥L2 + ∥enu ∥L2 )∥Au(tn+1 )∥L2 ∥∇ eun+1/2 ∥L2 1t ∥∇ eun+1/2 ∥2L2 + C Re1t (∥u(tn+1 ) − u(tn )∥2L2 + ∥enu ∥2L2 ). ≤ 8Re
By a similar method, we have I2 = 2S 1t |((b(tn+1 ) − b(tn )) × curl b(tn+1 ) + enb × curl b(tn+1 ), eun+1/2 )|
≤ CS 1t (∥b(tn+1 ) − b(tn )∥L2 + ∥enb ∥L2 )∥b(tn+1 )∥H 2 ∥∇ eun+1/2 ∥L2 1t ≤ ∥∇ eun+1/2 ∥2L2 + C ReS 2 1t (∥b(tn+1 ) − b(tn )∥2L2 + ∥enb ∥2L2 ) 8Re
and I3 = 2S 1t |(u(tn+1 ) × enb + u(tn+1 ) × (b(tn+1 ) − b(tn )), curl ebn+1 )|
≤ CS 1t (∥b(tn+1 ) − b(tn )∥L2 + ∥enb ∥L2 )∥Au(tn+1 )∥L2 ∥curl enb+1 ∥L2 S 1t ∥curl enb+1 ∥2L2 + CRmS 1t (∥b(tn+1 ) − b(tn )∥2L2 + ∥enb ∥2L2 ). ≤ 2Rm
Due to the fact that div enu = 0 in Ω and enu = 0 on ∂ Ω , we have I4 = 21t |(∇ p(tn+1 ), enu+1/2 )|
= 21t |(∇ p(tn+1 ), enu+1/2 − enu )| ≤ C 1t ∥enu+1/2 − enu ∥L2 ∥∇ p(tn+1 )∥L2 ≤
1 2
∥eun+1/2 − enu ∥2L2 + C (1t )2 .
From (2.9) and Young’s inequality, one has I5 = 21t |(Rnu+1 , enu+1/2 ) + S (Rnb+1 , enb+1 )|
≤ C 1t (∥Rnu+1 ∥H −1 ∥∇ enu+1/2 ∥L2 + S ∥Rnb+1 ∥H −1 ∥curl ebn+1 ∥L2 ) 1t S 1t ≤ ∥∇ eun+1/2 ∥2L2 + ∥curl enb+1 ∥2L2 + C 1t (Re∥Run+1 ∥2H −1 + RmS ∥Rbn+1 ∥2H −1 ). 8Re
2Rm
n+1/2
Using the definitions of enu+1 and eu n+1/2
eun+1 − eu
1t
−
1 Re
, we rewrite (3.3) as
∆(enu+1 − enu+1/2 ) − ∇ pn+1 = 0,
div enu+1 = 0.
(4.10)
Testing (4.10) by 21teun+1 leads to
∥enu+1 ∥2L2 + ∥enu+1 − enu+1/2 ∥2L2 − ∥enu+1/2 ∥2L2 +
1t Re
(∥∇ eun+1 ∥2L2 + ∥∇(eun+1 − eun+1/2 )∥2L2 − ∥∇ enu+1/2 ∥2L2 ) = 0.
(4.11)
Adding (4.9) and (4.11) and combining these estimates for I1 to I5 into the resulting inequality and using discrete Gronwall’s inequality (2.11), we conclude that there exists some positive constant C > 0 depending on Re, Rm, S such that (4.8) holds.
176
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
Remark 4.1. Assumption 2 and (4.8) imply that un+1/2 , un+1 and bn+1 are uniform boundedness in l∞ (H1 ):
∥∇ un+1 ∥l∞ (L2 ) + ∥∇ un+1/2 ∥l∞ (L2 ) + ∥bn+1 ∥l∞ (H 1 ) ≤ C ,
(4.12)
where C depends on Re, Rm and S. Using the above uniform boundedness (4.12), we can improve the orders of the temporal error estimates derived in Lemma 4.1 as follows. Lemma 4.2. Suppose that Assumptions 1–2 are satisfied. Then there exists some positive constant C > 0 such that
∥eun+1 ∥l∞ (L2 ) + S ∥enb+1 ∥l∞ (L2 ) +
1 Re
∥∇ enu+1 ∥l2 (L2 ) + ∥
S Rm
∥enb+1 ∥l2 (H 1 ) ≤ C 1t ,
eun+1
−
enu 2l2 (L2 )
∥
(4.13)
≤ C (1t ) , 3
(4.14)
where C depends on Re, Rm and S. Proof. Taking the sum of (4.6) and (4.10) gives eun+1 − enu
1t
−
1 Re
1enu+1 + ((u(tn+1 ) − u(tn )) · ∇)u(tn+1 ) + (enu · ∇)un+1/2 + (u(tn ) · ∇)eun+1/2
+ S (b(tn+1 ) − b(tn )) × curl b(tn+1 ) + Senb × curl bn+1 + Sb(tn ) × curl ebn+1 + ∇ epn+1 = Run+1 .
(4.15)
Subtracting (3.2) from (4.3) and using u(tn+1 ) × b(tn+1 ) − un+1/2 × bn = u(tn+1 ) × (b(tn+1 ) − b(tn )) + eun+1/2 × b(tn ) + un+1/2 × enb yields ebn+1 − enb
1t
+
1 Rm
curl curl enb+1 − curl (enu+1/2 × b(tn )) − curl (un+1/2 × enb )
− curl (u(tn+1 ) × (b(tn+1 ) − b(tn ))) = Rnb+1
(4.16)
and div ebn+1 = 0. Testing (4.15) and (4.16) by 21tenu+1 and 2S 1tebn+1 , respectively, using (3.6) and adding two resulting equations leads to
∥eun+1 ∥2L2 + ∥enu+1 − enu ∥2L2 − ∥enu ∥2L2 + S ∥enb+1 ∥2L2 2S 1t 21t + S ∥ebn+1 − enb ∥2L2 − S ∥enb ∥2L2 + ∥∇ enu+1 ∥2L2 + ∥curl enb+1 ∥2L2 Re
Rm
≤ 21t |b(u(tn+1 ) − u(tn ), u(tn+1 ), enu+1 )| + 2S 1t |((b(tn+1 ) − b(tn )) × curl b(tn+1 ), eun+1 )| + 2S 1t |(u(tn+1 ) × (b(tn+1 ) − b(tn )), curl enb+1 )| + 21t |(Run+1 , eun+1 ) + S (Rbn+1 , ebn+1 )| + 21t |b(u(tn ), enu+1/2 , enu+1 ) + b(enu , un+1/2 , enu+1 )| + 2S 1t |(un+1/2 × enb , curl enb+1 ) + (enb × curl bn+1 , enu+1 )| + 2S 1t |(b(tn ) × curl enb+1 , enu+1 − enu+1/2 )| = I6 + · · · + I12 .
(4.17)
By a similar method in the proof of Lemma 4.1, we can bound I6 + · · · + I9 as I6 + I7 + I8 + I9
= 21t |b(u(tn+1 ) − u(tn ), u(tn+1 ), enu+1 )| + 2S 1t |((b(tn+1 ) − b(tn )) × curl b(tn+1 ), enu+1 )| + 2S 1t |(u(tn+1 ) × (b(tn+1 ) − b(tn )), curl enb+1 )| + 21t |(Run+1 , eun+1 ) + S (Rbn+1 , ebn+1 )| S 1t 1t ∥∇ eun+1 ∥2L2 + ∥curl enb+1 ∥2L2 + C 1t Re∥u(tn+1 ) − u(tn )∥2L2 + (ReS 2 + RmS )∥b(tn+1 ) − b(tn )∥2L2 ≤ 4Re 4Rm + C 1t Re∥Rnu+1 ∥2H −1 + RmS ∥Rnb+1 ∥2H −1 . From (2.5)–(2.7), we estimate I10 as I10 = 21t |b(u(tn ), enu+1/2 , enu+1 ) + b(enu , un+1/2 , eun+1 )| 1/2
1/2
≤ C 1t ∥Au(tn )∥L2 ∥∇ enu+1 ∥L2 (∥enu+1/2 − enu ∥L2 + ∥enu ∥L2 ) + C 1t ∥∇ un+1/2 ∥L2 ∥∇ eun+1 ∥L2 ∥enu ∥L2 ∥∇ enu ∥L2 1t 1t ≤ ∥∇ enu+1 ∥2L2 + ∥∇ enu ∥2L2 + C Re1t (∥eun+1/2 − enu ∥2L2 + ∥enu ∥2L2 ). 4Re
2Re
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
177
Similarly, from (2.9), I11 is bounded by I11 = 2S 1t |(un+1/2 × enb , curl enb+1 ) + (enb × curl bn+1 , enu+1 )| 1/2
1/2
1/2
1/2
≤ CS 1t ∥∇ un+1/2 ∥L2 ∥curl enb+1 ∥L2 ∥enb ∥L2 ∥curl enb ∥L2 + CS 1t ∥curl bn+1 ∥L2 ∥∇ enu+1 ∥L2 ∥enb ∥L2 ∥curl enb ∥L2 S 1t S 1t 1t ∥∇ eun+1 ∥2L2 + ∥curl enb+1 ∥2L2 + ∥curl enb ∥2L2 + C (Re2 RmS 3 + Rm3 S )1t ∥enb ∥2L2 . ≤ 4Re
4Rm
2Rm
The last term I12 satisfies I12 = 2S 1t |(b(tn ) × curl enb+1 , enu+1 − enu+1/2 )|
≤ CS 1t ∥b(tn )∥H 2 ∥curl enb+1 ∥L2 ∥enu+1 − enu+1/2 ∥L2 S 1t ≤ ∥curl enb+1 ∥2L2 + CRmS 1t ∥enu+1 − enu+1/2 ∥2L2 . 4Rm
Combining these estimates into (4.17) and using the discrete Gronwall’s inequality (2.11), (4.4), (4.8), we complete the proof of (4.13) and (4.14). The temporal error estimate for the pressure in l2 (L2 ) is given in the next lemma by using the inf–sup condition. Lemma 4.3. Suppose that Assumptions 1–2 are satisfied. Then there exists some positive constant C > 0 such that
√ ∥epn+1 ∥l2 (L2 ) ≤ C 1t ,
(4.18)
where C depends on Re, Rm and S. Proof. From (4.15) and inf–sup condition, one has
∥epn+1 ∥L2 ≤ C sup
d(v, enp+1 )
v∈V
≤
∥∇ v∥L2
C
∥en+1 − enu ∥L2 + C ∥∇ enu+1 ∥L2 + C ∥Run+1 ∥H −1 + C (∥∇ enu ∥L2 + ∥∇ eun+1/2 ∥L2 + ∥u(tn+1 ) − u(tn )∥L2 ) 1t u + C (∥curl enb ∥L2 + ∥curl enb+1 ∥L2 + ∥b(tn+1 ) − b(tn )∥L2 ).
Taking the sum of the above inequality from 0 to N − 1 and using (4.4), (4.8) and (4.14), we complete the proof of Lemma 4.3. As a consequence of (4.14) and Assumption 2, we have
∥un+1 − un ∥l∞ (L2 ) + ∥bn+1 − bn ∥l∞ (L2 ) ≤ C 1t .
(4.19)
To give the spatial error estimates for the fully discrete fractional-step scheme (3.8)–(3.11), we need to show the regularity results for un+1/2 , un+1 , pn+1 and bn+1 . Lemma 4.4. Suppose that Assumptions 1–4 are satisfied. Then there exists some τ0 > 0 such that when 1t < τ0 , there holds 1 Re
∥Aun+1 ∥l2 (L2 ) +
1 Re
∥un+1/2 ∥l2 (H 2 ) + ∥∇ pn+1 ∥l2 (L2 ) +
1 Rm
∥bn+1 ∥l2 (H 2 ) ≤ C ,
(4.20)
where C depends on Re, Rm and S. Proof. Taking the sum of (3.1) and (3.3) leads to
−
1 Re
1un+1 + ∇ pn+1 = f(tn+1 ) − Dt un+1 − (un · ∇)un+1/2 − Sbn × curl bn+1 ,
(4.21)
div un+1 = 0.
(4.22)
Using Assumption 3, we get 1 Re
∥Aun+1 ∥L2 + ∥∇ pn+1 ∥L2 ≤ C (∥Dt un+1 ∥L2 + ∥f(tn+1 )∥L2 + S ∥bn × curl bn+1 ∥L2 + ∥(un · ∇)un+1/2 ∥L2 ).
The last two terms in the above inequality can be bounded, respectively, by S ∥bn × curl bn+1 ∥L2 ≤ S ∥enb × curl bn+1 ∥L2 + S ∥b(tn ) × curl bn+1 ∥L2
≤ CS (∥curl enb ∥L2 ∥bn+1 ∥H 2 + ∥curl bn+1 ∥L2 )
(4.23)
178
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
and 1/2
1/2
∥(un · ∇)un+1/2 ∥L2 ≤ C ∥∇ un ∥L2 ∥Aun ∥L2 ∥∇ un+1/2 ∥L2 ≤
1 2Re
∥Aun ∥L2 + C Re∥∇ un ∥L2 ,
where we use (4.12). Similarly, from Assumption 4 and (3.2) and the following formula curl (u × v) = (div v)u − (u · ∇)v + (v · ∇)u − (div u)v,
(4.24)
we have 1 Rm
∥bn+1 ∥H 2 ≤ C (∥Dt bn+1 ∥L2 + ∥(un+1/2 · ∇)bn ∥L2 ) + C (∥(bn · ∇)un+1/2 ∥L2 + ∥(div un+1/2 )bn ∥L2 ) ≤ C (∥Dt bn+1 ∥L2 + ∥∇ bn ∥L3 ∥∇ un+1/2 ∥L2 + ∥bn ∥L∞ ∥∇ un+1/2 ∥L2 ) 1/2
1/2
≤ C (∥Dt bn+1 ∥L2 + ∥curl bn ∥L2 ∥bn ∥H 2 ) ≤ C (∥Dt bn+1 ∥L2 + Rm∥curl bn ∥L2 ) +
1 2Rm
∥bn ∥H 2 .
(4.25)
Combining (4.23) and (4.25), and taking the sum from 0 to N − 1 yield
1t
N −1 1 n =0
≤ C 1t
Re
∥Aun+1 ∥2L2 + ∥∇ pn+1 ∥2L2 + 2
N −1
1
∥bn+1 ∥2H 2 2
Rm
∥Dt un+1 ∥2L2 + ∥Dt bn+1 ∥2L2 + ∥f(tn+1 )∥2L2 + C 1t ∥curl bn+1 ∥2L2 + ∥curl bn ∥2L2 + ∥∇ un ∥2L2
n =0
+ CS 2 1t
N −1
∥curl enb ∥2L2 ∥bn+1 ∥2H 2 +
n =0
N −1 1 1t
4
n=0
∥Aun ∥2L2 + 2
Re
1
n 2 ∥ b ∥ . H2 2
Rm
The last term in the above inequality can be absorbed on the left-hand side. From (4.13) and (4.19), we have
1t
N −1 1 n =0
Let τ0 = 1 Re
Re
1 . 2C6 S 2 Rm2
∥Aun+1 ∥2L2 + ∥∇ pn+1 ∥2L2 + 2
1
∥bn+1 ∥2H 2 2
Rm
≤ C + C6 S 2 (1t )2
N −1
∥bn+1 ∥2H 2 .
n =0
When 1t < τ0 , the above inequality implies that
∥Aun+1 ∥l2 (L2 ) + ∥∇ pn+1 ∥l2 (L2 ) +
1 Rm
∥bn+1 ∥l2 (H 2 ) ≤ C .
(4.26)
By using the regularity result [28] for elliptic problem and (4.26), we can easily get 1 Re
∥un+1/2 ∥l2 (H 2 ) ≤ C ,
if we notice (4.8) and
∥un+1/2 − un ∥L2 ≤ ∥enu+1/2 − enu ∥L2 + ∥u(tn+1 ) − u(tn )∥L2 . We complete the proof of (4.20).
Remark 4.2. From (4.20) and Assumption 2, we have 1 Re
∥Aeun+1 ∥l2 (L2 ) + ∥∇ enp+1 ∥l2 (L2 ) +
1 Rm
∥enb+1 ∥l2 (H 2 ) ≤ C .
(4.27)
Lemma 4.5. Suppose that Assumptions 1–4 are satisfied. Then when τ < τ0 , there holds
√ ∥ τ (tn+1 )∇ eun+1 ∥l∞ (L2 ) + ∥ τ (tn+1 )curl enb+1 ∥l∞ (L2 ) ≤ C 1t , where τ0 is defined in Lemma 4.4, and C depends on Re, Rm and S.
(4.28)
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
179
Proof. Testing (4.15) by −21tAenu+1 gives
∥∇ enu+1 ∥2L2 − ∥∇ enu ∥2L2 + ∥∇(enu+1 − enu )∥2L2 +
1t Re
∥Aenu+1 ∥2L2
≤ C 1t ∥Run+1 ∥2L2 + C 1t ∥((u(tn+1 ) − u(tn )) · ∇)u(tn+1 )∥2L2 + C 1t ∥(enu · ∇)un+1/2 + (u(tn ) · ∇)enu+1/2 ∥2L2 + C 1t ∥(b(tn+1 ) − b(tn )) × curl b(tn+1 )∥2L2 + C 1t ∥enb × curl bn+1 + b(tn ) × curl enb+1 ∥2L2 = C 1t ∥Rnu+1 ∥2L2 + I13 + · · · + I16 .
(4.29)
From Assumption 2, one has I13 = C 1t ∥((u(tn+1 ) − u(tn )) · ∇)u(tn+1 )∥2L2
≤ C 1t ∥Au(tn+1 )∥2L2 ∥∇(u(tn+1 ) − u(tn ))∥2L2 tn+1 ∥∇ ut ∥2L2 dt . ≤ C (1t )2 tn
Using (2.6) and (4.12), we get I14 = C 1t ∥(enu · ∇)un+1/2 + (u(tn ) · ∇)enu+1/2 ∥2L2
≤ C 1t (∥enu ∥2L∞ ∥∇ un+1/2 ∥2L2 + ∥Au(tn )∥2L2 ∥∇ enu+1/2 ∥2L2 ) ≤ C 1t (∥∇ enu ∥L2 ∥Aenu ∥L2 + ∥∇ enu+1/2 ∥2L2 ) ≤
1t 2Re
∥Aenu ∥2L2 + C 1t (∥∇ enu ∥2L2 + ∥∇ enu+1/2 ∥2L2 ).
By a similar argument for the above two inequalities, we have I15 = C 1t ∥(b(tn+1 ) − b(tn )) × curl b(tn+1 )∥2L2
≤ C 1t ∥b(tn+1 )∥2H 2 ∥b(tn+1 ) − b(tn )∥2H 1 tn+1 2 ≤ C (1t ) ∥bt ∥2H 1 dt , tn
and I16 = C 1t ∥enb × curl bn+1 + b(tn ) × curl enb+1 ∥2L2
≤ C 1t (∥enb ∥2L∞ ∥curl bn+1 ∥2L2 + ∥b(tn )∥2H 2 ∥curl enb+1 ∥2L2 ) ≤
C7 1t 2Rm
∥enb ∥2H 2 + C 1t (∥curl enb ∥2L2 + ∥curl enb+1 ∥2L2 ).
Combining these estimates into (4.29) yields
∥∇ eun+1 ∥2L2 − ∥∇ enu ∥2L2 + ≤ C 1t ∥Run+1 ∥2L2 + + C (1t )
2
1t 2Re
tn+1
1t Re
∥Aenu+1 ∥2L2
∥Aenu ∥2L2 + ut 2L2 dt
∥∇ ∥ tn
C7 1 t 2Rm
tn+1
+
∥enb ∥2H 2 + C 1t (∥∇ enu ∥2L2 + ∥∇ eun+1/2 ∥2L2 + ∥curl enb ∥2L2 + ∥curl enb+1 ∥2L2 ) bt 2H 1 dt
∥ ∥ tn
.
(4.30)
Testing (4.7) by 21tcurl curl enb+1 leads to
∥curl ebn+1 ∥2L2 − ∥curl enb ∥2L2 +
C7 1 t Rm
∥enb+1 ∥2H 2
≤ C 1t ∥Rbn+1 ∥2L2 + C 1t ∥curl (enu+1/2 × bn )∥2L2 + C 1t ∥curl (u(tn+1 ) × enb )∥2L2 + C 1t ∥curl (u(tn+1 ) × (b(tn+1 ) − b(tn )))∥2L2 = C 1t ∥Rnb+1 ∥2L2 + I17 + I18 + I19 .
(4.31)
180
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
Using (2.5), (2.6) and (4.24), we obtain I17 = C 1t ∥curl (enu+1/2 × bn )∥2L2
≤ C 1t (∥(enu+1/2 · ∇)bn ∥2L2 + ∥(bn · ∇)enu+1/2 ∥2L2 + ∥(div eun+1/2 )bn ∥2L2 ) ≤ C 1t (∥∇ enu+1/2 ∥2L2 ∥enb ∥H 2 ∥curl enb ∥L2 + ∥∇ eun+1/2 ∥2L2 ∥b(tn )∥2H 2 ) ≤ C (1t )2 ∥enb ∥2H 2 + C ∥curl enb ∥2L2 + C 1t ∥∇ enu+1/2 ∥2L2 . Similarly, we have I18 = C 1t ∥curl (u(tn+1 ) × enb )∥2L2
≤ C 1t (∥(u(tn+1 ) · ∇)enb ∥2L2 + ∥(enb · ∇)u(tn+1 )∥2L2 ) ≤ C 1t ∥curl enb ∥2L2 ∥Au(tn+1 )∥2L2 ≤ C 1t ∥curl enb ∥2L2 , and I19 = C 1t ∥curl (u(tn+1 ) × (b(tn+1 ) − b(tn )))∥2L2
≤ C 1t (∥(u(tn+1 ) · ∇)(b(tn+1 ) − b(tn ))∥2L2 + ∥((b(tn+1 ) − b(tn )) · ∇)u(tn+1 )∥2L2 ) ≤ C 1t ∥b(tn+1 ) − b(tn )∥2H 1 ∥Au(tn+1 )∥2L2 tn+1 ≤ C (1t )2 ∥bt ∥2H 1 dt . tn
Taking the sum of (4.30) and (4.31) and combining the above estimates into the resulting inequality give
∥∇ eun+1 ∥2L2 + ∥curl enb+1 ∥2L2 +
1t Re
∥Aenu+1 ∥2L2 +
C7 1t Rm
∥enb+1 ∥2H 2
≤ ∥∇ enu ∥2L2 + ∥curl enb ∥2L2 + C 1t (∥Rnu+1 ∥2L2 + ∥Rnb+1 ∥2L2 ) +
1t 2Re
∥Aenu ∥2L2 +
C7 1t 2Rm
∥enb ∥2H 2
+ C 1t 2 ∥enb ∥2H 2 + C ∥curl enb ∥2L2 + C 1t (∥∇ enu ∥2L2 + ∥∇ enu+1/2 ∥2L2 + ∥curl enb ∥2L2 + ∥curl ebn+1 ∥2L2 ) tn+1 tn+1 + C (1t )2 ∥∇ ut ∥2L2 dt + ∥bt ∥2H 1 dt . tn
tn
Multiplying it by τ (tn+1 ) and using τ (tn+1 ) ≤ τ (tn ) + 1t and |τ (tn+1 )| ≤ 1, we derive that
τ (tn+1 )∥∇ eun+1 ∥2L2 + τ (tn+1 )∥curl enb+1 ∥2L2 +
1t Re
τ (tn+1 )∥Aenu+1 ∥2L2 +
C7 1 t Rm
τ (tn+1 )∥ebn+1 ∥2H 2
≤ τ (tn )∥∇ enu ∥2L2 + τ (tn )∥curl enb ∥2L2 + 1t (∥∇ enu ∥2L2 + ∥curl enb ∥2L2 ) + C 1t (τ (tn+1 )∥Run+1 ∥2L2 + τ (tn+1 )∥Rbn+1 ∥2L2 ) S 1t n 2 1t τ (tn )∥Aenu ∥2L2 + τ (tn ) ∥eb ∥H 2 + C (1t )2 (∥Aenu ∥2L2 + ∥enb ∥2H 2 ) + 4Re
4Rm
+ C 1t (∥∇ enu ∥2L2 + ∥∇ enu+1/2 ∥2L2 + ∥curl enb ∥2L2 + ∥curl ebn+1 ∥2L2 ) tn+1 tn+1 + C (1t )2 ∥∇ ut ∥2L2 dt + ∥bt ∥2H 1 dt + C ∥curl enb ∥2L2 . tn
Taking the sum of (4.32) from 0 to m ≤ N − 1, and using (4.5), (4.8), (4.13), (4.20), we complete the proof of (4.28). Combining the results in this section, we obtain the temporal error estimates (3.12)–(3.16). 5. Spatial error analysis To show the spatial error estimates, we introduce the projectors (Rh , πh ) : V × M −→ Vh × Mh defined by a1 (Rh u − u, vh ) − d(vh , πh p − p) = 0
∀ vh ∈ Vh , d(Rh u − u, qh ) = 0 ∀ qh ∈ Mh .
In addition, define Kh : V −→ Vh and Qh : X −→ Xh by
(∇(Kh u − u), ∇ vh ) = 0 ∀ vh ∈ Vh , a2 (Qh b − b, wh ) +
(4.32)
tn
1 Rm
(div (Qh b − b), div wh ) = 0 ∀ wh ∈ Xh .
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
181
Then for all (u, p, b) ∈ V ∩ H2 (Ω ) × M ∩ H 1 (Ω ) × X ∩ H2 (Ω ), the following approximation and stability properties hold:
∥u − Rh u∥L2 + h∥∇(u − Rh u)∥L2 ≤ Ch2 ∥u∥H 2 ,
(5.1)
∥u − Kh u∥L2 + h∥∇(u − Kh u)∥L2 ≤ Ch ∥u∥H 2 ,
(5.2)
∥b − Qh b∥L2 + h∥b − Qh b∥H 1 ≤ Ch ∥b∥H 2 , ∥p − πh p∥L2 ≤ Ch∥p∥H 1 , ∥∇ Rh u∥L6 + ∥∇ Kh u∥L6 + ∥Qh b∥W 1,6 + ∥πh p∥L6 ≤ C .
(5.4)
2 2
(5.3) (5.5)
For 0 ≤ n ≤ N − 1, let us denote errors by Enu+1/2 = Kh un+1/2 − un+1/2 ,
Enu = Rh un − un ,
n+1/2
enuh = Rh un − unh ,
euh
n+1/2
= Kh un+1/2 − uh
Enb = Qh bn − bn ,
,
enbh = Qh bn − bnh ,
n+1/2
From (3.1)–(3.3) and (3.8)–(3.11), we obtain that euh
n+1/2
− enuh
euh
1t
n+1/2
, vh + a1 (euh
, vh ) =
n+1/2
Eu
Epn = πh pn − pn , enph = πh pn − pnh .
n +1 , enbh+1 , (euh , enph+1 ) satisfy the following problems:
− Enu
1t
n+1/2
, vh + b1 (unh , uh
, vh )
− b1 (un , un+1/2 , vh ) + S (bnh × curl bnh+1 , vh ) − S (bn × curl bn+1 , vh ) ∀ vh ∈ Vh ,
(5.6)
and
n+1 ebh − enbh
1t
, wh +
(
a2 enbh+1 n+1/2
+ (un+1/2 × bn − uh
, wh ) +
1 Rm
(div
enbh+1
, wh ) =
Enb+1 − Enb
1t
, wh
× bnh , curl wh ) ∀ wh ∈ Xh ,
(5.7)
and
n+1/2
n+1 euh − euh
=
n+1/2
, vh + a1 (enuh+1 − euh
1t
n+1/2
Eun+1 − Eu
1t
n +1 , vh ) − d(vh , eph )
, vh − a1 (Enu+1/2 , vh ) ∀ vh ∈ Vh ,
(5.8)
and n +1 , qh ) = 0, d(euh
∀ qh ∈ Mh .
(5.9)
The main result in this section is to show the following spatial error estimate based on the regularity estimate (4.20). Lemma 5.1. Suppose that Assumptions 1–4 are satisfied. If the time step 1t satisfies 1t ≥ Ch, then there exists some h0 > 0 such that when h < h0 , the following spatial error estimate holds: +1 ∥enuh+1 ∥l∞ (L2 ) + S ∥em bh ∥l∞ (L2 ) +
1 Re
∥∇ enuh+1 ∥l2 (L2 ) +
S Rm
n +1 ∥ebh ∥l2 (H 1 ) ≤ Ch,
(5.10)
where C depends on Re, Rm and S. n+1/2
Proof. Taking vh = 2euh
∥
in (5.6) leads to
n+1/2 2 euh L2
n+1/2
∥ − ∥enuh ∥2L2 + ∥euh
=2
n+1/2
Eu
− enuh ∥2L2
1t
− Enu
1t
n+1/2
, euh
+
2 Re
n+1/2 2 L2
∥∇ euh
∥
n+1/2 n+1/2 n+1/2 + 2 b1 (unh , uh , euh ) − b1 (un , un+1/2 , euh ) n+1/2
+ 2S bnh × curl bnh+1 − bn × curl bn+1 , euh
= J1 + J2 + J3 .
(5.11)
From Hölder inequality and Young’s inequality, J1 satisfies n+1/2
J1 ≤ 2(1t )−1 (∥Enu+1/2 ∥L2 + ∥Enu ∥L2 )∥euh
≤
1 4Re
n+1/2 2 L2
∥∇ euh
∥L2
∥ + Ch4 (1t )−2 (∥un+1/2 ∥2H 2 + ∥un ∥2H 2 ).
182
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184 n+1/2
The definitions of enuh , euh n+1/2
2 b1 (unh , uh
and the trilinear term b1 imply that
n+1/2
, euh
n+1/2
) − b1 (un , un+1/2 , euh
)
n+1/2 n+1/2 n+1/2 = 2 b1 (Enu − enuh , Enu+1/2 , euh ) + b1 (Enu − enuh , un+1/2 , euh ) + b1 (un , Eun+1/2 , euh ) . By using (5.1), (5.2) and (4.20), we estimate the three terms in the above equation as follows: n+1/2
b1 (Enu − enuh , Enu+1/2 , euh
≤ C (∥∇ ∥ + ∥∇ Enu L2
1
≤
) n+1/2
∥ )∥∇ Enu+1/2 ∥L2 ∥∇ euh
enuh L2
∥L2
n+1/2 2 L2
∥ + Ch2 ∥un+1/2 ∥H 2 + C Reh2 ∥∇ enuh ∥2L2 ∥un+1/2 ∥2H 2 ,
∥∇ euh
6Re
and n+1/2
b1 (Enu − enuh , un+1/2 , euh
≤ C (∥ ∥ + ∥ Enu L2
1
≤
∥∇
6Re
) n+1/2
∥ )∥un+1/2 ∥H 2 ∥∇ euh
enuh L2
∥L2
n+1/2 2 euh L2
∥ + Ch2 ∥un+1/2 ∥H 2 + C Re∥enuh ∥2L2 ∥un+1/2 ∥2H 2 ,
and n+1/2
b1 (un , Eun+1/2 , euh
n+1/2
) ≤ C ∥∇ un ∥L2 ∥∇ Enu+1/2 ∥L2 ∥∇ euh ≤
1
∥L2
n+1/2 2 L2
∥ + Ch2 ∥un+1/2 ∥H 2 .
∥∇ euh
6Re
Therefor, J2 satisfies n+1/2
n+1/2
J2 = 2 b1 (unh , uh 1
≤
2Re
n+1/2
) − b1 (un , un+1/2 , euh
, euh
)
n+1/2 2 L2
∥ + Ch2 ∥un+1/2 ∥H 2 + C Re(h2 ∥∇ enuh ∥2L2 + ∥enuh ∥2L2 )∥un+1/2 ∥2H 2 .
∥∇ euh
By a similar method, we have n+1/2
J3 = 2S bnh × curl bnh+1 − bn × curl bn+1 , euh n+1/2
= −2S (bnh × curl enbh+1 , euh + 2S ((Enb − enbh ) × curl ≤ −2S ( +
bnh
1 4Re
× curl
enbh+1
,
n+1/2
) − 2S (enbh × curl Enb+1 , euh
n+1/2 bn+1 euh
,
n+1/2 euh
) + Ch∥
) n+1/2
) + 2S (Enb × curl Enb+1 , euh
enbh 2L2
∥ ∥
n+1/2
) + 2S (bn × curl Enb+1 , euh
)
bn+1 2H 2
∥
n+1/2 2 L2
∥ + Ch2 ∥bn+1 ∥H 2 + C ∥enbh ∥2L2 ∥bn+1 ∥2H 2 ,
∥∇ euh
where we use the inverse inequality ∥vh ∥L3 ≤ Ch−1/6 ∥vh ∥L2 . Combining these estimates into (5.11) yields n+1/2 2 L2
n+1/2
∥ − ∥enuh ∥2L2 + ∥euh
∥euh
1t
n+1/2
≤ −2S (bnh × curl enbh+1 , euh
− enuh ∥2L2
+
1 Re
n+1/2 2 L2
∥∇ euh
∥
) + Ch4 (1t )−2 (∥un+1/2 ∥2H 2 + ∥un ∥2H 2 ) + Ch2 (∥un+1/2 ∥H 2 + ∥bn+1 ∥2H 2 )
+ C Reh2 ∥∇ enuh ∥2L2 ∥un+1/2 ∥2H 2 + C Re∥enuh ∥2L2 ∥un+1/2 ∥2H 2 + C ∥enbh ∥2L2 ∥bn+1 ∥2H 2 .
(5.12)
n +1 Taking vh = 2euh in (5.8) and qh = enph+1 in (5.9) leads to n+1/2 2 L2
n+1 2 ∥euh ∥L2 − ∥euh
n+1/2 2 L2
∥ + ∥enuh+1 − euh
∥
1t
=
n+1/2
Eun+1 − Eu
1t
+
1 Re
n+1/2 2 L2
(∥∇ enuh+1 ∥2L2 − ∥∇ euh
,
enuh+1
− a1 (Enu+1/2 , enuh+1 )
≤ Ch4 (1t )−2 (∥un+1/2 ∥2H 2 + ∥un+1 ∥2H 2 ) +
1 2Re
n+1/2
∥ + ∥∇(enuh+1 − euh
∥∇ enuh+1 ∥2L2 + Ch2 ∥un+1/2 ∥2H 2 ,
)∥2L2 )
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
183
which together with (5.12) yields n +1 2 ∥euh ∥L2 − ∥enuh ∥2L2
1
n+1 2 + ∥∇ euh ∥L2 1t Re n + 1 / 2 ≤ −2S (bnh × curl enbh+1 , euh ) + Ch4 (1t )−2 (∥un+1/2 ∥2H 2 + ∥un ∥2H 2 + ∥un+1 ∥2H 2 )
+ Ch2 (∥un+1/2 ∥H 2 + ∥bn+1 ∥2H 2 ) + C Reh2 ∥∇ enuh ∥2L2 ∥un+1/2 ∥2H 2 + C Re∥enuh ∥2L2 ∥un+1/2 ∥2H 2 + C ∥enbh ∥2L2 ∥bn+1 ∥2H 2 .
(5.13)
Setting wh = 2Senbh+1 in (5.7) leads to n+1 2 S ∥ebh ∥L2 − S ∥enbh ∥2L2 + S ∥enbh+1 − enbh ∥2L2
2SC
∥en+1 ∥2 1 1t Rm bh H n+1/2 n +1 ≤ 2S (un+1/2 × enbh , curl enbh+1 ) − 2S (un+1/2 × Enb , curl enbh+1 ) + 2S (euh × bnh , curl ebh ) +
n +1 − 2S (Enu+1/2 × bn , curl enbh+1 ) + 2S (Enu+1/2 × enbh , curl enbh+1 ) − 2S (Enu+1/2 × Enb , curl ebh )
SC ∥en+1 ∥2 1 + Ch4 (1t )−2 (∥bn ∥2H 2 + ∥bn+1 ∥2H 2 ) 2Rm bh H SC n+1/2 n +1 ∥en+1 ∥2 1 + Ch4 (1t )−2 (∥bn ∥2H 2 + ∥bn+1 ∥2H 2 ) + 2S (euh × bnh , curl ebh ) + C ∥enbh ∥2L2 ∥un+1/2 ∥2H 2 ≤ Rm bh H
+
+ C ∥Eun+1/2 ∥2L2 ∥bn ∥2H 2 + C ∥un+1/2 ∥2H 2 ∥Enb ∥2L2 + C ∥enbh ∥2L3 ∥∇ un+1/2 ∥2L2 + C ∥∇ Enu+1/2 ∥2L2 ∥Enb ∥L2 ∥Enb ∥H 1 .
(5.14)
Under the condition 1t ≥ Ch, adding (5.13) and (5.14), and taking the sum from n = 0 to n = m ≤ N − 1, we get
∥
+1 2 em uh L2
∥ + S∥
∥ + 1t
+1 2 em bh L2
m 1 n =0
≤ Ch2 1t
m
Re
∥∇
enuh+1 2L2
∥ +
SC Rm
∥
enbh+1 2H 1
+ C 1t
∥
m ∥un+1/2 ∥2H 2 + ∥un+1 ∥2H 2 + ∥bn+1 ∥2H 2 + C Reh2 1t ∥∇ enuh ∥2L2 ∥un+1/2 ∥2H 2
n =0 m
n =0
Re∥un+1/2 ∥2H 2 ∥enuh ∥2L2 + ∥un+1/2 ∥2H 2 ∥enbh ∥2L2 + ∥bn+1 ∥2H 2 ∥enbh ∥2L2
n =0
≤ Ch2 + C8 Reh1t
m
∥∇ enuh ∥2L2 + C Re1t
n =0
+ C 1t
m
m n+1/2 2 ∥u ∥H 2 ∥enuh ∥2L2 n=0
∥un+1/2 ∥2H 2 ∥enbh ∥2L2 + ∥bn+1 ∥2H 2 ∥enbh ∥2L2 .
n =0
Let h be sufficiently small such that C8 h < complete the proof of Lemma 5.1.
1 . Re2
Applying the discrete Gronwall’s inequality to the above inequality, we
Combining the temporal error estimates (3.12)–(3.16) with the results in Lemma 5.1, we obtain the temporal–spatial error estimates (3.17) and (3.18). 6. Conclusions In this paper we propose a fractional-step scheme for solving the stationary MHD problems in three-dimensional bounded domains, which is an extension of the previous work for the Navier–Stokes problems studied in [24]. The main feature is that the nonlinearities for velocity and the incompressibility in (1.1)–(1.2) are split into two steps. Compared to the projection scheme [14], the fractional-step scheme studied in this paper allows to enforce the original boundary conditions of the problems in all substeps of the fractional-step scheme, and allows for a discrete energy estimate. We have proved that the time-discrete fractional-step √ scheme provides the temporal error estimates of O (1t ) for the velocity and the magnetic field in l2 (H1 ) ∩ l∞ (L2 ) and O ( 1t ) for the pressure in l2 (L2 ). For the fully discrete FEM scheme, the temporal–spatial error estimates of O (1t + h) for the velocity and the magnetic field in l2 (H1 ) ∩ l∞ (L2 ) are obtained under √the constraint 1t ≥ Ch. However, for the projection scheme, the temporal–spatial error orders derived in [14] are only O ( 1t + h) for the velocity and the magnetic field in l2 (H1 ). Although some temporal–spatial convergence orders are derived in this paper, there are still two questions worthy of study. One is whether the spatial error estimate for the velocity and the magnetic field in l∞ (L2 ) can be improved to
184
R. An, C. Zhou / Journal of Computational and Applied Mathematics 313 (2017) 168–184
be optimal. Another is whether there is an optimal error estimate for the pressure in l2 (L2 ). We observed that GuillénGonzález and Redondo-Neble have studied optimal error estimates of the fractional-step scheme for the Navier–Stokes problems in [29,30]. But it is much more complicated to study these optimal error estimates for the MHD problem due to the appearance of the magnetic field equations. In addition, the extension of this fractional-step scheme to higher-order temporal accuracy also an interesting problem, could possibly be achieved using a Crank–Nicolson scheme. However, there needs much higher regularity of the solution to the MHD problems. These objectives are worth studying in our future work.
Acknowledgments The authors would like to thank the anonymous reviewers for their careful reviews and valuable comments. These comments are all valuable and very helpful for improving this manuscript. This manuscript is supported by Zhejiang Provincial Natural Science Foundation with Grant Nos. LY16A010017 and LY14A010020. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30]
W. Hughes, F. Young, The Electromagnetics of Fluids, Wiley, New York, 1966. R. Moreau, Magneto-hydrodynamics, Kluwer Academic Publishers, 1990. J. Gerbeau, C. Le Bris, T. Lelièvre, Mathematical Methods for the Magnetohydrodynamics of Liquid Metals, Oxford University Press, Oxford, 2006. M. Gunzburger, A. Meir, J. Peterson, On the existence, uniqueness, and finite element approximation of solutions of the equations of stationary, incompressible magnetohydrodynamics, Math. Comp. 56 (1991) 523–563. J. Douglas Jr., J. Wang, An absolutely stabilized finite element method for the Stokes problem, Math. Comp. 52 (1989) 495–508. J. Gerbeau, A stabilized finite elemenet method for the incompressible magnetohydrodynamic equations, Numer. Math. 87 (2000) 83–111. Y. He, Unconditional convergence of the Euler semi-implicit scheme for the three-dimensional incompressible MHD equations, IMA J. Numer. Anal. 35 (2015) 767–801. D. Schötzau, Mixed finite element methods for stationary incompressible magneto-hydrodynamics, Numer. Math. 96 (2004) 771–800. S. Badia, R. Codina, R. Planas, On an unconditionally convergent stabilized finite element approximation of resistive magnetohydrodynamics, J. Comput. Phys. 234 (2013) 399–416. S. Badia, R. Planas, J. Gutiérrez-Santacreu, Unconditionally stable operator splitting algorithms for the incompressible magnetohydrodynamics system discretized by a stabilized finite element formulation based on projections, Internat. J. Numer. Methods Engrg. 93 (2013) 302–328. M. Belenli, S. Kaya, L. Rebholz, N. Wilson, A subgrid stabilization finite element method for incompressible magnetohydrodynamics, Int. J. Comput. Math. 90 (2013) 1506–1523. R. Codina, N. Hernández, Approximation of the thermally coupled MHD problem using a stabilized finite element method, J. Comput. Phys. 230 (2011) 1281–1303. 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. A. Prohl, Convergent finite element discretizations of the nonstationary incompressible magnetohydrodynamic system, ESAIM:M2AN 42 (2008) 1065–1087. D. Sondak, J. Shadid, A. Oberai, R. Pawlowski, E. Cyr, T. Smith, A new class of finite element variational multiscale turbulence models for incompressible magnetohydrodynamics, J. Comput. Phys. 295 (2015) 596–616. G. Yuksel, R. Ingram, Numerical analysis of a finite element Crank–Nicolson discretization for MHD flows at small magnetic Reynolds numbers, Int. J. Numer. Anal. Model. 10 (2013) 74–98. G. Yuksel, O. Isik, Numerical analysis of Backward-Euler discretization for simplified magnetohydrodynamic flow, Appl. Math. Model. 39 (2015) 1889–1898. Y. Zhang, Y. Hou, L. Shan, Numerical analysis of the Crank–Nicolson extrapolation time discrete scheme for magnetohydrodynamics flows, Numer. Methods Partial Differential Equations 31 (2015) 2169–2208. A. Chorin, Numerical solution of the Navier–Stokes equations, Math. Comp. 22 (1968) 745–762. R. Temam, Sur l’approximation de la solution des equations de Navier–Stokes par la méthode des pas fractionnaires II, Arch. Ration. Mech. Anal. 33 (1969) 377–385. A. Prohl, On pressure approximation via projection methods for nonstationary incompressible Navier–Stokes equations, SIAM J. Numer. Anal. 47 (2008) 158–180. R. Temam, Remark on the pressure boundary condition for the projection method, Theor. Comput. Fluid Dyn. 3 (1991) 181–184. J. Blasco, R. Codina, A. Huerta, A fractional-step method for the incompressible Navier–Stokes equations related to a predictor-multicorrector algorithm, Internat. J. Numer. Methods Fluids 28 (1997) 1391–1419. J. Blasco, R. Codina, Error estimates for an operator-splitting method for incompressible flows, Appl. Numer. Math. 51 (2004) 1–17. R. Temam, Navier–Stokes Equations, North-Holland Publishing Company, Amsterdam, 1977. R. Adams, Sobolev Spaces, Academic Press, New York, 1975. J. Heywood, R. Rannacher, Finite-element approximation of the nonstationary Navier–Stokes problem Part IV: error analysis for second-order time discretization, SIAM J. Numer. Anal. 27 (1990) 353–384. D. Gilbarg, N.S. Trudinger, Elliptic Partial Differential Equations of Second Order, Springer-Verlag, 1983. F. Guillén-González, M.V. Redondo-Neble, New error estiamtes for a viscosity-splitting scheme in time for the three-dimensional Navier–Stokes equations, IMA J. Numer. Anal. 31 (2011) 556–579. F. Guillén-González, M.V. Redondo-Neble, Spatial error estimates for a finite element viscosity-splitting scheme for the Navier–Stokes equations, Int. J. Numer. Anal. Model. 10 (2013) 826–844.