- Email: [email protected]
Error analysis of first-order projection method for time-dependent magnetohydrodynamics equations
Accepted Manuscript Error analysis of first-order projection method for time-dependent magnetohydrodynamics equations
Rong An, Yuan Li
PII: DOI: Reference:
S0168-9274(16)30204-5 http://dx.doi.org/10.1016/j.apnum.2016.10.010 APNUM 3116
To appear in:
Applied Numerical Mathematics
Received date: Revised date: Accepted date:
25 August 2015 20 September 2016 2 October 2016
Please cite this article in press as: R. An, Y. Li, Error analysis of first-order projection method for time-dependent magnetohydrodynamics equations, Appl. Numer. Math. (2016), http://dx.doi.org/10.1016/j.apnum.2016.10.010
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Error analysis of first-order projection method for timedependent magnetohydrodynamics equations Rong An∗
Yuan Li
College of Mathematics and Information Science, Wenzhou University, Wenzhou 325025, P.R.China.
Abstract. This paper focuses on a linearized fully discrete projection scheme for time-dependent magnetohydrodynamics equations in three-dimensional bounded domain. It is shown that the proposed projection scheme allows for a discrete energy inequality and is unconditionally stable. In addition, we present a rigorous analysis for the rates of convergence. AMS subject classifications: 65M12, 76W05 Key words: Magnetohydrodynamics equations, Projection method, Stability, Rate of convergence, Temporal errors.
1
Introduction
The incompressible magnetohydrodynamics (MHD) equations are used to describe the flow of a viscous, incompressible and electrically conducting fluid. For the understanding of the physical background of the MHD equations, we refer to Hughes [10] and Moreau [11]. Let Ω ⊂ R3 be a bounded and simply-connected domain which is either convex or has a C1,1 boundary ∂Ω. We consider the non-dimensional MHD equations in the primitive variable formulation ∂u 1 − Δu + (u · ∇)u + ∇ p + Sb × curl b = f, ∂t Re div u = 0, div b = 0, ∂b 1 + curl (curl b) − curl (u × b) = 0, ∂t Rm
(1.1) (1.2) (1.3)
for x ∈ Ω and t ∈ (0, T ) with T > 0, where 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 equations (1.1)-(1.3) couple the incompressible Navier-Stokes ∗ Corresponding
author. Email: [email protected], [email protected] (Rong An) 1
c
2
R.An
equations with Maxwell’s equations. Thus, the unknowns in (1.1)-(1.3) are the fluid velocity u, the pressure p and the magnetic field b. The MHD equations should be completed by the appropriate initial and boundary conditions. For the sake of simplicity, we consider the following initial and boundary conditions: u( x, 0) = u0 , u = 0,
b · n = 0,
b( x, 0) = b0 curl b × n = 0
in Ω,
(1.4)
on ∂Ω × [0, T ],
(1.5)
where n denotes the unit outward normal vector on ∂Ω. The initial vector functions u0 and b0 satisfy the compatibility condition div u0 = 0 and div b0 = 0. It was observed that testing (1.1) and (1.3) by u and Sb, respectively, and adding the resulting equations leads to the following energy identity: 1 S 1 d |u|2 + S|b|2 dx + |∇u|2 + |curl b|2 dx = f · udx ∀ t > 0. 2 dt Ω Re Rm Ω Ω Thus, for any prescribed initial data (u0 , b0 ) ∈ L2 (Ω) × L2 (Ω), the MHD problem (1.1)-(1.5) exists the global weak solutions. On the regularities of the weak solutions, Sermange-Temam in [14] established the existences of local strong solution with large initial data and global strong solution with small initial data. For the numerical analysis of the MHD problem, the mixed finite element approximations were first proposed and studied for the stationary MHD problem in [7], where H1 - conforming elements were used to discretize the magnetic field provided that Ω is either convex or has a C1,1 boundary ∂Ω. Inspired by the absolutely stable methods for Stokes problem in [3], a stabilized mixed finite element method for stationary MHD problem was developed by Gerbeau [4]. For the time-dependent MHD equations (1.1)-(1.3), recently, He proposed a linearized semi-implicit Euler scheme in [8], where L2 -unconditional convergence of this scheme was proved by using the negative norm technique. Other fully discrete Crank-Nicolson schemes were studied in [19, 20]. 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 corresponding magnetic field. A mixed finite element formulation based on H(curl)¨ elements (or N´ed´elec elements) for b was proposed and studied by Schotzau in [13] for the stationary MHD problem. For a review of various numerical methods for the MHD equations, we refer to Gerbeau-Bris-Leli`evre [5]. In the present work, we will consider the three-dimensional time-dependent MHD equations (1.1)-(1.3) with the initial and boundary conditions (1.4) and (1.5), and propose a linearized fully discrete scheme based on a modified Chorin’s projection scheme for Navier-Stokes equations in [15]. The projection methods were first proposed by Chorin [2] and Temam [17], and have been further developed in various directions. The main advantage of the projection methods 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. 2
3
Projection method for MHD equations
For a review of projection methods for the Navier-Stokes equations, we refer to [6]. To state our main results derived in this paper, we recall the following terminology [15]: Definition 1.1. Let X be a Banach space equipped with norm || · || X and f : [0, T ] −→ X is continuous. Let 0 = t0 < t1 < · · · < t N = T be a uniform partition of the time interval [0, T ] with time step Δt = T/N and tn = nΔt for 0 ≤ n ≤ N. We say f Δt is a weakly order α approximation of f in X if there exists a constant C independent of Δt such that Δt
N
∑ || f Δt (tn ) − f (tn )||2X ≤ C(Δt)2α ;
n =0
and we say f Δt is a strongly order α approximation of f in X if there exists a constant C independent of Δt and n such that max || f Δt (tn ) − f (tn )||2X ≤ C (Δt)2α .
0≤ n ≤ N
To our best knowledge, Prohl in [12] first proposed a projection scheme for the MHD problem. Unfortunately, the projection scheme in [12] does not allow for a discrete energy estimate. Moreover, it was proved that this scheme provided the weakly order 12 approximations of the velocity field and the magnetic field in H1 (Ω). In this paper, we propose a new linearized projection scheme, which allows for a discrete energy inequality, thus, and is unconditionally stable. We will prove that this new projection scheme provides the weakly first-order approximations of the velocity field and the magnetic field in H1 (Ω), and the strongly first-order approximations of the velocity field and the magnetic field in L2 (Ω) under the regularity assumptions 1 and 2 below. To obtain the optimal convergence order of the pressure, we need some further regularity assumption as like in [16], under which we can derive the weakly first-order approximation of the pressure in L2 (Ω). This paper is organized as follows: in the next section, we begin with some notation and lay out some regularity 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 inequality is derived in Section 3. The proofs of main results are given in Sections 4 and 5, which are split into several theorems.
2 Mathematical setting For the mathematical setting of the MHD equations (1.1)-(1.3) with the initial and boundary conditions (1.4)-(1.5), we introduce some function spaces and their associated norms. For k ∈ N + , 1 ≤ p ≤ +∞, let W k,p (Ω) denote the standard Sobolev 3
4
R.An
space. The norm in W k,p (Ω) is defined by
||u||W k,p =
⎧ ⎪ ⎪ ⎨ ( ⎪ ⎪ ⎩
∑
| β|≤k Ω
∑
|∇ β u| p dx )1/p
for 1 ≤ p < +∞,
sup |∇ β u|
for p = +∞,
| β|≤k Ω
where ∇ is the differential operator with respect to x and
∇β =
∂| β| β
β
β
∂x1 1 ∂x2 2 ∂x3 3
for the multi-index β = ( β 1 , β 2 , β 3 ) and | β| = β 1 + β 2 + β 3 with β 1 , β 2 , β 3 ≥ 0. We k,p define W0 (Ω) to be the subspace of W k,p (Ω) of functions with zero trace on ∂Ω. When p = 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,p (Ω) and L p (Ω) are used to denote the vector Sobolev spaces H k (Ω)3 , W k,p (Ω)3 and L p (Ω)3 , respectively. In particular, (·, ·) denotes the L2 (Ω) inner product. Let X be a Banach space. For some T > 0 or T = ∞, denote Bochner spaces L p (0, T; X ), 1 ≤ p < +∞ the spaces of measurable functions from the interval [0, T ] into X such that T 0
p
||u(t)|| X dt < +∞.
If p = +∞, the functions in L∞ (0, T; X ) are required to satisfy sup ||u(t)|| X < +∞.
t∈[0,T ]
The symbols C, C1 , C2 , · · · are used to denote the generic positive constants independent of the time step Δt. We define H = {u ∈ L2 (Ω), div u = 0 in Ω, u · n = 0 on ∂Ω}, V = H10 (Ω),
V0 = {u ∈ V, div u = 0 in Ω},
X = {u ∈ H1 (Ω), u · n = 0 on ∂Ω}, M=
L20 (Ω)
X0 = {u ∈ X, div u = 0 in Ω},
= { q ∈ L ( Ω ), 2
Ω
qdx = 0}.
Thanks to Poincar`e inequality, the norm in V is defined by
||u||V =
|∇u| dx 2
Ω
4
1/2 .
5
Projection method for MHD equations
Define the bilinear forms
a2 (u, v) =
1 Rm
1 a1 (u, v) = Re
Ω
∇u : ∇vdx
∀ u, v ∈ V,
(curl u · curl v + div udiv v)dx
∀ u, v ∈ X,
Ω
d(v, p) =
Ω
∀ v ∈ V, p ∈ M
pdiv vdx
and the trilinear form
1 (div u)v · wdx b(u, v, w) = (u · ∇)v · wdx + 2 Ω Ω 1 1 = (u · ∇)v · wdx − (u · ∇)w · vdx 2 Ω 2 Ω
∀ u, v, w ∈ V.
It is obvious that b(u, v, w) satisfies b(u, v, w) = −b(u, w, v)
∀ u, v, w ∈ V.
(2.1)
With the above notations, given u0 , b0 ∈ H and f ∈ L∞ (0, T; L2 (Ω)), the weak variational formulation of the MHD problem (1.1)-(1.5) is defined to find (u, p, b) ∈ L2 (0, T; V) × L2 (0, T; M) × L2 (0, T; X) such that
(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) − (u × b, curl w) = 0
∀ w ∈ X,
(2.3)
where we use the following formula: Ω
curl u · vdx =
Ω
u · curl vdx +
∂Ω
n × u · vds.
(2.4)
Corresponding to (1.1), we recall Stokes operator A [18]. Let us denote the orthogonal projection operator from L2 (Ω) onto H by PH . Then Stokes operator A is defined by Au = −PH Δu ∀ u ∈ D( A) = V0 ∩ H2 (Ω). The emphasis of this paper is to prove the temporal error estimates of the projection scheme (3.1)-(3.3) stated in Section 3. Then we make the following assumptions on the initial data and the regularities of the solution. Assumption 1: The initial data u0 , b0 and f satisfy
|| Au0 || L2 + ||b0 || H2 + sup {||f(t)|| L2 + ||ft (t)|| L2 } ≤ C. 0≤ t ≤ T
5
(2.5)
6
R.An
Assumption 2: The MHD problem (1.1)-(1.5) admits a unique local strong solution (u, p, b) on (0, T ) such that sup {|| Au(t)|| L2 + ||b(t)|| H2 + || p(t)|| H1 + ||ut (t)|| L2 + ||bt (t)|| L2 }
0≤ t ≤ T
+
T 0
||ut ||2H1 + ||bt ||2H1 + || pt ||2H1 + ||utt ||2L2 + ||btt ||2L2 dt ≤ C. (2.6)
The following known inequalities are frequently used [1, 18]:
||PH v|| H1 ≤ C ||v|| H1 ∀ v ∈ H1 (Ω), ||v|| H1 ≤ C (||curl v|| L2 + ||div v|| L2 ) ∀ v ∈ X, ||v|| Lr ≤ C ||v|| H1 (2 ≤ r ≤ 6) ∀ v ∈ H (Ω), 1
||v|| L3 ≤ ||v|| L∞ ≤
C ||v||1/2 ||v||1/2 L2 H1 1/2 1/2 C ||v|| H1 ||v|| H2
||v|| H2 ≤ C || Av|| L2
(2.7) (2.8) (2.9)
∀ v ∈ H ( Ω ),
(2.10)
∀ v ∈ H (Ω), i = 1, 2,
(2.11)
1 2
∀ v ∈ D( A ).
(2.12)
Finally, we recall a discrete version of Gronwall inequality established in [9]. Lemma 2.1. Let ak , bk , ck and γk , for integers k ≥ 0, be the nonnegative numbers such that an + τ
n
n
n
k =0
k =0
k =0
∑ bk ≤ τ ∑ γk a k + τ ∑ c k + B
for n ≥ 0.
(2.13)
Suppose that τγk < 1, for all k, and set σk = (1 − τγk )−1 . Then an + τ
n
∑
k =0
bk ≤ exp(τ
n
∑
k =0
γk σk )(τ
n
∑ ck + B)
k =0
for n ≥ 0.
(2.14)
Remark 2.1. If the first sum on the right in (2.13) extends only up to n − 1, then estimate (2.14) holds for all k > 0 with σk = 1.
3
A projection scheme and main results
Let 0 = t0 < t1 < · · · < t N = T be a uniform partition of the time interval [0, T ] with time step Δt = T/N and tn = nΔt, 0 ≤ n ≤ N, where [0, T ] is the maximal time interval such that a local strong solution exists. For 0 ≤ n ≤ N, let us denote un = u(tn ), pn = p(tn ), bn = b(tn ) and fn = f(tn ). For any sequence of functions { gn }nN=0 , we define Dt gn+1 = ( gn+1 − gn )/Δt for 0 ≤ n ≤ N − 1. We start with 0 = u0 , B0 = b0 and an arbitrary P0 ∈ M ∩ H 1 (Ω). We solve (U n , Bn ) and U0 = U (Un , Pn ) for 1 ≤ n ≤ N, respectively, by n − Un −1 U 1 n n + ∇ Pn−1 + SBn−1 × curl Bn = fn , − ΔU + (Un−1 · ∇)U Δt Re 1 n n × Bn−1 ) = 0, div Bn = 0 Dt B + curl (curl Bn ) − curl (U Rm 6
(3.1) (3.2)
7
Projection method for MHD equations
n = 0, Bn · n = 0 and curl Bn × n = 0 on ∂Ω, and with the boundary condition U n
Un − U + 2∇( Pn − Pn−1 ) = 0, Δt
div Un = 0,
(3.3)
with the boundary condition Un · n = 0 on ∂Ω. A reformulation of the problem (3.3) is 1 n div U Δt
2ΔPn = 2ΔPn−1 +
(3.4)
with the boundary condition ∇( Pn − Pn−1 ) · n = 0 on ∂Ω, and n − 2Δt∇( Pn − Pn−1 ). Un = U
(3.5)
Since the problems (3.1)-(3.3) consist of two linear equations, its solvability should follow easily if we can establish a priori energy estimates for its solutions. The following discrete energy identity just serves this purpose. n , Bn ) and (Un , Pn ) be the solution of (3.1)-(3.3). Then Theorem 3.1. For 1 ≤ n ≤ N, let (U the following energy identity holds:
||Um ||2L2 + S||Bm ||2L2 + 2(Δt)2 ||∇ Pm ||2L2 + Δt +
m
∑ (||U
n
n =1
m
2
∑ ( Re ||∇U
n =1
n
|| L2 +
2S ||curl Bn ||2L2 ) Rm
1 n − Un−1 ||2L2 + ||U − Un ||2L2 + S||Bn − Bn−1 ||2L2 ) 2
= ||u0 ||2L2 + S||b0 ||2L2 +2(Δt)2 ||∇ P0 ||2L2 + 2Δt
m
∑ (fn , U
n =1
n
) (3.6)
for every 1 ≤ m ≤ N. n
and 2SΔtBn , respectively, and adding the Proof. Testing (3.1) and (3.2) by 2ΔtU resulting equations yield n ||2 2 + ||U n − Un−1 ||2 2 − ||Un−1 ||2 2 + S(||Bn ||2 2 − ||Bn−1 ||2 2 + ||Bn − Bn−1 ||2 2 ) ||U L L L L L L 2Δt 2SΔt n n n n n 2 n − 1 ) = 2Δt(f , U ). (3.7) || L2 + ||∇U ||curl B || L2 + 2Δt(∇ P , U + Re Rm Here we use the following identity:
(u × curl v, w) = (w × u, curl v) Testing (3.5) by Un and lead to
∀ u, v, w ∈ H1 (Ω).
(3.8)
1 n n (U + U ), respectively, and adding the resulting equations 2
n ||2 2 + 1 ||Un − U n ||2 2 = Δt(∇( Pn−1 − Pn ), U n ). ||Un ||2L2 − ||U L L 2 7
(3.9)
8
R.An
The sum of (3.7) and (3.9) yields n − Un−1 ||2 2 − ||Un−1 ||2 2 + 1 ||Un − U n ||2 2 ||Un ||2L2 + ||U L L L 2 2Δt n || L2 + 2SΔt ||curl Bn ||2 2 + S(||Bn ||2L2 − ||Bn−1 ||2L2 + ||Bn − Bn−1 ||2L2 ) + ||∇U L Re Rm n n n n −1 n = 2Δt(f , U ) − Δt(∇( P + P ), U ). (3.10) n = Un + 2Δt∇( Pn − Pn−1 ) and div Un = 0, we have Observing U n ) = −2(Δt)2 (||∇ Pn ||2 2 − ||∇ Pn−1 ||2 2 ). −Δt(∇( Pn + Pn−1 ), U L L Putting the above equation into (3.10) and taking the sum for n from 1 to m ≤ N, we complete the proof of (3.6). As a consequence of (3.6), we have the following result. Corollary 3.1. The projection scheme (3.1)-(3.3) is unconditionally stable. Moreover, there exists some C > 0 such that
||Um ||2L2 + S||Bm ||2L2 + 2(Δt)2 ||∇ Pm ||2L2 + Δt m
m
1
∑ ( Re ||∇U
n =1
n
|| L2 +
S ||curl Bn ||2L2 ) Rm
− Un−1 ||2 2 + 1 ||U n − Un ||2 2 + S||Bn − Bn−1 ||2 2 ) ≤ C (3.11) + ∑ (||U L L L 2 n =1 n
for every 1 ≤ m ≤ N. Let us denote enu = un − Un ,
n, enu = un − U
enb = bn − Bn ,
P n = Pn−1 + 2( Pn − Pn−1 ),
enp = pn − Pn ,
enp = pn − P n .
The main results derived in this paper are presented in the following two theorems. Theorem 3.2. Suppose that Assumptions 1-2 are satisfied. Then there exists some positive constant C > 0 such that sup
||um − Um ||2L2 + ||bm − Bm ||2L2 ≤ C (Δt)2 ,
(3.12)
+ ||bn − Bn ||2H1 + Δt|| pn − P n ||2L2 ) ≤ C (Δt)2 .
(3.13)
1≤ m ≤ N
Δt
N
∑ (||∇(un − Un )||2L
n =1
2
8
9
Projection method for MHD equations
From (3.13), we can see that the weakly order 12 approximation of the pressure in L2 (Ω) is derived. However, this convergence order is not optimal. To improve it, we need some further regularity assumptions as like in [16]. Now, we assume that sup (||∇ut (t)|| L2 + ||curl bt (t)|| L2 + ||utt (t)|| L2 + ||btt (t)|| L2 + || pt (t)|| H1 )
0≤ t ≤ T
+
T 0
(|| Aut ||2L2 + ||bt ||2H2 + ||uttt ||2L2 + ||bttt ||2L2 + || ptt ||2H1 )dt ≤ C. (3.14)
Then we can prove the wekly first-order approximation of the pressure in L2 (Ω). Theorem 3.3. Suppose that Assumptions 1-2 and (3.14) are satisfied. If
||∇( P0 − p0 )|| L2 ≤ CΔt, then there exists some C > 0 such that Δt
4
N
∑ || pn − P n ||2L
n =1
2
≤ C (Δt)2 .
(3.15)
Proof of Theorem 3.2
For 1 ≤ n ≤ N, taking t = tn in (1.1)-(1.3), we obtain 1 un − un −1 − Δun + (un · ∇)un + ∇ pn + Sbn × curl bn = fn − Rnu , Δt Re
(4.1)
and bn − bn −1 1 − curl (curl bn ) − curl (un × bn ) = − Rnb , Δt Rm
(4.2)
where Rnu =
1 Δt
tn t n −1
(s − tn−1 )∂tt u(s)ds,
Rnb =
1 Δt
tn t n −1
(s − tn−1 )∂tt b(s)ds.
Subtracting (3.1) and (3.2) from (4.1) and (4.2), respectively, leads to 1 enu − enu−1 n − (un · ∇)un + ∇( Pn−1 − pn ) − Δ en = (Un−1 · ∇)U Δt Re u + S(Bn−1 × curl Bn − bn × curl bn ) − Rnu , (4.3) and enb − enb −1 1 n × Bn −1 ) − R n , + curl (curl enb ) = curl (un × bn ) − curl (U b Δt Rm 9
(4.4)
10
R.An
An alternative to (3.5) is enu = 2Δt∇( Pn − Pn−1 ). enu −
(4.5)
n
− (un · ∇)un and S(Bn−1 × curl Bn − bn × curl bn ) in (4.3) can The terms (Un−1 · ∇)U be rewritten as n − (un · ∇)un (Un−1 · ∇)U
(4.6)
enu − (eun−1 · ∇)un − (un−1 · ∇) enu − ((un − un−1 ) · ∇)un = (enu−1 · ∇) and S(Bn−1 × curl Bn − bn × curl bn )
(4.7)
= −S(Bn−1 × curl enb + ebn−1 × curl bn + (bn − bn−1 ) × curl bn ). n
× Bn−1 ) in (4.4) can be rewritten as The term curl (un × bn ) − curl (U n × Bn −1 ) curl (un × bn ) − curl (U
(4.8)
enu × Bn−1 ) + curl (un × (bn − bn−1 )) + curl (un × ebn−1 ). = curl (
Theorem 4.1. Suppose that Assumptions 1-2 are satisfied. Then there exists some positive constant C1 > 0 such that 2 m 2 2 m 2 ||em u || L2 + ||eb || L2 + ( Δt ) ||∇ e p || L2 + Δt
+ Proof. Testing (4.3) by
2Δt enu
m
∑ (||∇ enu ||2L
n =1
m
∑ (|| enu − enu−1 ||2L
n =1
2
2
+ ||∇enu ||2L2 + ||curl enb ||2L2 )
+ ||enu − enu ||2L2 ) ≤ C1 (Δt)2 . (4.9)
and using (u − v, 2u) = |u|2 − |v|2 + |u − v|2 , we have
2Δt enu ||2L2 ||∇ Re ≤ 2Δt|b(un − un−1 , un , enu ) + b(enu−1 , un , enu )| + 2Δt|( Rnu , enu )| enu ||2L2 − ||eun−1 ||2L2 + || enu − eun−1 ||2L2 + || n −1
n
+ 2SΔt(B × curl B − b × curl b = I1 + I2 + I3 + I4 . n
n
, enu )
+ 2Δt(∇( P
n −1
(4.10)
−p
n
), enu )
From (2.6) and the Young’s inequality, we bound I1 as enu || L2 I1 ≤ CΔt(||un − un−1 || L2 + ||enu−1 || L2 )|| Aun || L2 ||∇ Δt enu ||2L2 + CΔt(||un − un−1 ||2L2 + ||eun−1 ||2L2 ). ≤ ||∇ 4Re It is easy to show that enu || L2 ≤ I2 ≤ CΔt|| Rnu || L2 ||∇
Δt enu ||2L2 + CΔt|| Rnu ||2L2 . ||∇ 4Re
10
11
Projection method for MHD equations
¨ inequality, Poincar´e inequality and the From (2.9), (4.7), H1 (Ω) → L6 (Ω), Holder Young’s inequality, I3 is bounded by I3 = −2SΔt(Bn−1 × curl enb , enu ) − 2SΔt(enb −1 × curl bn , enu )
− 2SΔt((bn − bn−1 ) × curl bn , enu )
≤ −2SΔt(Bn−1 × curl enb , enu ) + CΔt(||enb −1 || L2 + ||bn − bn−1 || L2 )||curl bn || L3 ||∇ enu || L2 Δt enu ||2L2 − 2SΔt(Bn−1 × curl enb , ≤ ||∇ enu ) 4Re + CΔt(||enb −1 ||2L2 + ||bn − bn−1 ||2L2 ). Combining these estimates into (4.10) leads to Δt enu ||2L2 ||∇ Re ≤ CΔt(||ebn−1 ||2L2 + ||bn − bn−1 ||2L2 ) − 2SΔt(Bn−1 × curl enb , enu ) enu ||2L2 − ||enu−1 ||2L2 + || enu − enu−1 ||2L2 + ||
(4.11)
+ CΔt(||un − un−1 ||2L2 + ||eun−1 ||2L2 + || Run ||2L2 ) + 2Δt(∇( Pn−1 − pn ), enu ).
1 Testing (4.5) by enu and (enu + enu ), respectively, we obtain 2 1 enu ||2L2 + ||enu − enu ||2L2 ) = 0, (||enu ||2L2 − || 2 1 enu ||2L2 ) = Δt(∇( Pn − Pn−1 ), (||enu ||2L2 − || enu ), 2 where we use div enu = 0 in Ω and enu = 0 on ∂Ω. Taking the sum of the above identities leads to 1 enu ||2L2 + ||enu − enu ||2L2 = Δt(∇( Pn − Pn−1 ), ||enu ||2L2 − || enu ), 2 which together with (4.11) implies that 1 Δt enu − enu−1 ||2L2 + ||enu − enu ||2L2 + enu ||2L2 ||enu ||2L2 − ||enu−1 ||2L2 + || ||∇ 2 Re ≤ CΔt(||enb −1 ||2L2 + ||bn − bn−1 ||2L2 ) − 2SΔt(Bn−1 × curl enb , enu )
+ CΔt(||un − un−1 ||2L2 + ||eun−1 ||2L2 + || Rnu ||2L2 )
− Δt(∇(enp + enp−1 + pn − pn−1 ), enu ). From (4.5) again, one has
enu = enu + 2Δt∇( Pn−1 − Pn )
= enu + 2Δt∇(enp − enp−1 − ( pn − pn−1 )). 11
(4.12)
12
R.An
Then we have
− Δt(∇(enp + enp−1 + pn − pn−1 ), enu ) = −2(Δt)2 (∇(enp + enp−1 + pn − pn−1 ), ∇(enp − enp−1 − ( pn − pn−1 )) = 2(Δt)2 (||∇enp−1 ||2L2 − ||∇enp ||2L2 ) + 2(Δt)2 ||∇( pn − pn−1 )||2L2 + 4(Δt)2 (∇enp−1 , ∇( pn − pn−1 )) ≤ 2(Δt)2 (||∇enp−1 ||2L2 − ||∇enp ||2L2 ) + 2(Δt)2 ||∇( pn − pn−1 )||2L2 + 2Δt||∇( pn − pn−1 )||2L2 + 2(Δt)3 ||∇enp−1 ||2L2 . Substituting the above inequality into (4.12) we obtain
||enu ||2L2 − ||enu−1 ||2L2 + 2(Δt)2 (||∇enp ||2L2 − ||∇enp−1 ||2L2 ) 1 Δt + || enu − enu−1 ||2L2 + ||enu − enu ||2L2 + enu ||2L2 ||∇ 2 Re ≤ CΔt(||ebn−1 ||2L2 + ||bn − bn−1 ||2L2 ) − 2SΔt(Bn−1 × curl enb , enu ) + CΔt(||un − un−1 ||2L2 + ||enu−1 ||2L2 + || Rnu ||2L2 )
+ 4Δt||∇( pn − pn−1 )||2L2 + 2(Δt)3 ||∇enp−1 ||2L2 . Using (4.8) and testing (4.4) by 2SΔtenb leads to 2SΔt ||curl enb ||2L2 Rm enu × Bn−1 , curl enb ) + (un × (bn − bn−1 ), curl enb ) = 2SΔt(( S||enb ||2L2 − S||ebn−1 ||2L2 + S||enb − enb −1 ||2L2 +
+ (un × enb −1 , curl enb )) − 2SΔt( Rbn , enb )
≤ 2SΔt( enu × Bn−1 , curl enb ) + CSΔt|| Rbn || L2 ||curl enb ||2L2
+ CSΔt(||bn − bn−1 || L2 + ||enb −1 || L2 )||un || H2 ||curl enb || L2 SΔt enu × Bn−1 , curl enb ) ≤ ||curl enb ||2L2 + 2SΔt( Rm + CΔt(||bn − bn−1 ||2L2 + ||enb −1 ||2L2 + || Rnb ||2L2 ). From (3.8), the sum of the above two inequalities yields
||enu ||2L2 − ||eun−1 ||2L2 + S||enb ||2L2 − S||ebn−1 ||2L2 + 2(Δt)2 (||∇enp ||2L2 − ||∇enp−1 ||2L2 ) 1 Δt 2SΔt enu − enu−1 ||2L2 + ||enu − enu ||2L2 + enu ||2L2 + ||∇ ||curl enb ||2L2 + || 2 Re Rm ≤ CΔt(||ebn−1 ||2L2 + ||enu−1 ||2L2 + (Δt)2 ||∇enp−1 ||2L2 ) + CΔt(|| Rnu ||2L2 + || Rnb ||2L2 ) + CΔt(||un − un−1 ||2L2 + ||bn − bn−1 ||2L2 + ||∇( pn − pn−1 )||2L2 ). Taking the sum of the above inequality for n from 1 to m ≤ N and using the discrete Gronwall inequality (2.14), we conclude that there exists some positive constant C1 > 12
13
Projection method for MHD equations
0 such that 2 m 2 2 m 2 ||em u || L2 + ||eb || L2 + ( Δt ) ||∇ e p || L2 + Δt
+
m
∑ (||∇ enu ||2L
n =1
m
∑ (|| enu − enu−1 ||2L
n =1
2
2
+ ||curl enb ||2L2 )
enu ||2L2 ) ≤ C1 (Δt)2 (4.13) + ||enu −
by using N
∑ (||un − un−1 ||2L
n =1
2
+ ||bn − bn−1 ||2L2 + ||∇( pn − pn−1 )||2L2 ) ≤ CΔt, N
∑ (|| Rnu ||2L
2
n =1
+ || Rnb ||2L2 ) ≤ CΔt.
(4.14) (4.15)
Applying the projector PH to (4.5) and using (2.7), (4.13) implies the inequality (4.9) and completes the proof of Theorem 4.1. Theorem 4.1 implies that the projection scheme (3.1)-(3.3) provides the weakly and strongly first-order approximations of the velocity field and the magnetic field in H1 (Ω) and L2 (Ω), respectively. Based on the above theorem, the rate of convergence for the pressure is stated in the following theorem. Theorem 4.2. Suppose that Assumptions 1-2 are satisfied. Then there exists some positive constant C2 > 0 such that Δt
N
∑ || enp ||2L
n =1
2
≤ C2 Δt.
(4.16)
Proof. Taking the sum of (3.1) and (3.3) gives Un − Un −1 1 n n + ∇ P n + SBn−1 × curl Bn = fn . − ΔU + (Un−1 · ∇)U Δt Re Subtracting (4.17) from (4.1) yields
(4.17)
1 enu − enu−1 n − (un · ∇)un − ∇ enp − Δ en = (Un−1 · ∇)U Δt Re u + SBn−1 × curl Bn − Sbn × curl bn − Rnu . (4.18) From (4.18) and inf-sup condition [18], one has enp || L2 ≤ C sup ||
v∈V,v =0
≤
d(v, enp )
||∇v|| L2
C n enu || L2 + C || Rnu || L2 ||e − enu−1 || L2 + C ||∇ Δt u enu || L2 + ||∇eun−1 || L2 + ||∇ enu || L2 + ||un − un−1 || L2 ) + C (||∇enu−1 || L2 ||∇
+ C (||ebn−1 || H1 ||curl enb || L2 + ||curl enb || L2 + ||ebn−1 || L2 + ||bn − bn−1 || L2 ). 13
14
R.An
Taking the sum of the above inequality from 1 to N and using (3.11), (4.9), (4.14) and (4.15), we have Δt
N
∑ || enp ||2L
2
n =1
≤ CΔt + C (Δt)2
N
∑ (||enb −1 ||2H
1
n =1
+ ||∇eun−1 ||2L2 ),
which together with (3.11) and (4.9) completes the proof of Theorem 4.2.
5
Proof of Theorem 3.3
As a consequence of (4.9) and Assumption 2, we have
|| Dt Un || L2 + || Dt Bn || L2 ≤ C
for 1 ≤ n ≤ N.
(5.1)
n , Un , Bn ) can be derived in next lemma. Then, the regularities of (U Lemma 5.1. Suppose that Assumptions 1-2 are satisfied. Then there exists some τ0 > 0 such that when Δt < τ0 , there holds Δt
N
∑ (||Un ||2H
n =1
2
n ||2 2 + ||Bn ||2 2 ) ≤ C. + ||U H H
(5.2)
Proof. From (3.1) and the regularity of elliptic problem, one has n
|| H2 ≤ C (|| Dt Un || L2 + (Δt)−1 || enu − enu || L2 + ||fn || L2 ||U n
|| L2 ). + ||∇ Pn−1 || L2 + ||Bn−1 × curl Bn || L2 + ||(Un−1 · ∇)U The last three terms in the above inequality can be bounded, respectively, by
||∇ Pn−1 || L2 ≤ ||∇enp−1 || L2 + ||∇ pn−1 || L2 and
||Bn−1 × curl Bn || L2 ≤ ||enb −1 × curl Bn || L2 + ||bn−1 × curl Bn || L2 ≤ C (||ebn−1 || H1 ||Bn || H2 + ||Bn || H1 )
and n
n
n
|| L2 ≤ ||(eun−1 · ∇)U || L2 + ||(un−1 · ∇)U || L2 ||(Un−1 · ∇)U n n || H2 + ||∇U || L2 ). ≤ C (||∇eun−1 || L2 ||U Similarly, from (3.2) and the following formulation curl (u × v) = (div v)u − (u · ∇)v + (v · ∇)u − (div u)v, 14
(5.3)
15
Projection method for MHD equations
we have n · ∇)Bn−1 || L2 + ||(Bn−1 · ∇)U n || L2 + ||(div U n )Bn−1 || L2 ) ||Bn || H2 ≤ C (|| Dt Bn || L2 + ||(U n || H2 + ||∇U n || L2 ). ≤ C (|| Dt Bn || L2 + ||en−1 || H1 ||U b
Taking the sum of the above inequalities leads to N
∑ (||U
Δt
+ || P n ||2H1 + ||Bn ||2H2 )
∑ (|| Dt Un ||2L
2
∑ (||∇enp−1 ||2L
2
≤ Δt +
n =1 N
n 2 || H2
n =1 N
n =1
+ Δt
enu − enu ||2L2 ) + || Dt Bn ||2L2 + (Δt)−2 || n ||2 2 ) + ||∇ pn−1 ||2L2 + ||fn ||2L2 + ||Bn ||2H1 + ||∇U L
N
∑ (||ebn−1 ||2H ||U 1
n =1
≤ C + C3 (Δt)2
N
∑ (||U
n =1
n 2 || H2
n 2 || H2
n
||2 2 ) + ||enb −1 ||2H1 ||Bn ||2H2 + ||∇eun−1 ||2L2 ||U H
+ ||Bn ||2H2 ).
Let τ0 < 1/2C3 . Then when Δt < τ0 , we complete the proof of Lemma 4.2.
Let us denote enu − eun−1 , r np = enp − enp−1 , bn = enb − enb −1 , un = enu − eun−1 , un =
∀ 1 ≤ n ≤ N.
Then (4.9) implies
||un ||2L2 + ||bn ||2L2 + Δt
m
∑ (||∇un ||2L
n =1
2
un ||2L2 + ||curl bn ||2L2 ) ≤ C (Δt)2 . + ||∇
(5.4)
Lemma 5.2. Suppose that Assumptions 1-2 and (3.14) are satisfied. If
||∇( P0 − p0 )|| L2 ≤ CΔt, then there exists some C4 > 0 such that
||enu − enu−1 || L2 ≤ C4 (Δt)2 ,
∀ 1 ≤ n ≤ N.
(5.5)
Proof. Taking the difference of (4.4) with two consecutive indices, we have bn − bn−1 1 + curl (curl bn ) = T n − T n−1 − ( Rnb − Rbn−1 ), Δt Rm 15
(5.6)
16
R.An
where enu × Bn−1 ) + curl (un × (bn − bn−1 )) + curl (un × ebn−1 ). T n = curl ( By a simple calculation, we find that T n − T n−1 in (5.6) is equal to enu−1 × (bn−1 − bn−2 )) − curl ( eun−1 × bn−1 ) T n − T n−1 = curl ( n × Bn−1 ) + curl (
+ curl ((un − un−1 ) × (bn − bn−1 )) + curl (un−1 × (bn − 2bn−1 + bn−2 )) + curl ((un − un−1 ) × enb −1 ) + curl (un−1 × bn−1 ). (5.7)
Testing (5.6) by 2SΔtbn yields 2SΔt ||curl bn ||2L2 Rm = 2SΔt(T n − T n−1 , bn ) − 2SΔt( Rnb − Rbn−1 , bn ). S||bn ||2L2 − S||bn−1 ||2L2 + S||bn − bn−1 ||2L2 +
(5.8)
From (5.7), (3.14) and the Young’s inequality, we get 2SΔt(T n − T n−1 , bn ) un × Bn−1 , curl bn ) + 2SΔt( enu−1 × (bn−1 − bn−2 ), curl bn ) = 2SΔt( enu−1 × bn−1 , curl bn ) + 2SΔt((un − un−1 ) × (bn − bn−1 ), curl bn ) − 2SΔt(
+ 2SΔt(un−1 × (bn − 2bn−1 + bn−2 ), curl bn )
+ 2SΔt((un − un−1 ) × ebn−1 , curl bn ) + 2SΔt(un−1 × bn−1 , curl bn ) SΔt ≤ 2SΔt( ||curl bn ||2L2 un × Bn−1 , curl bn ) + 2Rm eun−1 ||2L2 + CΔt||bn−1 − bn−2 || L2 ||curl (bn−1 − bn−2 )|| L2 ||∇ n −1 2 || H2 )||bn−1 ||2L2 + CΔt|| Aun−1 ||2L2 ||bn + CΔt||un − un−1 || L2 ||∇(un − un−1 )|| L2 ||curl (bn − bn−1 )||2L2 + CΔt||un − un−1 || L2 ||∇(un − un−1 )|| L2 ||curl ebn−1 ||2L2
+ CΔt(|| Aun−1 ||2L2 + ||U
un × Bn−1 , curl bn ) + ≤ 2SΔt(
SΔt ||curl bn ||2L2 + C (Δt)5 2Rm
enu−1 ||2L2 + ||curl enb −1 ||2L2 ) + CΔt(|| Aun−1 ||2L2 + ||U + C (Δt)3 (||∇
− 2bn−1 + bn−2 ||2L2
n −1 2 || H2 )||bn−1 ||2L2 ,
where we use
||bn − 2bn−1 + bn−2 ||2L2 = ||
tn t n −1
(tn+1 − t)btt (t)dt −
4 ≤ (Δt)4 ||btt ||2L2 . 3 16
t n −1 t n −2
(tn+1 − t)btt (t)dt||2L2
17
Projection method for MHD equations
In addition, by Taylor formula with the integral residue, we have
|| Rnb − Rbn−1 ||2L2 = ||
1 n (b − 2bn−1 + bn−2 ) − (bnt − btn−1 )||2L2 Δt
≤ 2(Δt)
3
tn
t n −2
||bttt (t)||2L2 dt.
Then 2SΔt( Rnb − Rnb −1 , bn ) ≤
SΔt ||curl bn ||2L2 + C (Δt)4 2Rm
tn t n −2
||bttt (t)||2L2 dt.
Combining these estimates into (5.8) leads to S||bn ||2L2 − S||bn−1 ||2L2 + S||bn − bn−1 ||2L2 +
SΔt ||curl bn ||2L2 Rm
un × Bn−1 , curl bn ) + C (Δt)5 + C (Δt)4 ≤ 2SΔt( enu−1 ||2L2 + ||curl enb −1 ||2L2 ) + C (Δt)3 (||∇
tn
t n −2
||bttt (t)||2L2 dt
(5.9)
n−1 ||2 2 )||n−1 ||2 2 . + CΔt(|| Aun−1 ||2L2 + ||U H L b
Similarly, taking the difference of (4.18) with two consecutive indices, we have 1 un − un−1 − Δ n = −∇r np−1 − 2∇(r np − r np−1 ) Δt Re u + ∇( pn − 2pn−1 + pn−2 ) + (X n − X n−1 ) − (W n − W n−1 ) − ( Rnu − Run−1 ), (5.10) where enu − (eun−1 · ∇)un − (un−1 · ∇) enu − ((un − un−1 ) · ∇)un , X n = (enu−1 · ∇) W n = S(Bn−1 × curl enb + ebn−1 × curl bn + (bn − bn−1 ) × curl bn ). Thanks to the proof of Lemma A2 in [16], we find Δt un ||2L2 ||∇ Re + (Δt)2 (||∇r np ||2L2 − ||∇r np−1 ||2L2 + ||∇(r np − r np−1 )||2L2 )
||un ||2L2 − ||un−1 ||2L2 + ||un − un−1 ||2L2 +
≤ 2Δt(X n − X n−1 , un ) − 2Δt(W n − W n−1 , un ) + C (Δt)3 ||∇r np−1 ||2L2 + C (Δt)4
tn
t n −2
(5.11)
(||uttt (t)||2L2 + ||∇ ptt (t)||2L2 )dt.
Since X n − X n−1 = (enu−1 · ∇) n + (un−1 · ∇) eun−1 − (eun−1 · ∇)(un − un−1 ) un − ((un−1 − un−2 ) · ∇) − (un−1 · ∇)un−1 − (un−1 · ∇) eun−1
− ((un − 2un−1 + un−2 ) · ∇)un + ((un−1 − un−2 ) · ∇)(un − un−1 ), 17
18
R.An
and W n − W n−1 = SBn−1 × curl bn − Sbn−1 × curl enb −1 + S(bn−1 − bn−2 ) × curl ebn−1
+ Senb −1 × curl (bn − bn−1 ) + Sbn−1 × curl bn−1
+ S(bn − bn−1 ) × curl (bn − bn−1 ) + S(bn − 2bn−1 + bn ) × curl bn−1 , the terms 2Δt(X n − X n−1 , un ) and −2Δt(W n − W n−1 , un ) in (5.11) can be bounded as follows: Δt n−1 ||2 2 ) 2Δt(X n − X n−1 , un ||2L2 + CΔt||un−1 ||2L2 (|| Aun−1 ||2L2 + ||U un ) ≤ ||∇ H 4Re n −1 2 n n −1 2 n −1 2 + CΔt||eu || L2 || A(u − u )|| L2 + CΔt||u || L2 || Aun−1 ||2L2 eun−1 ||2L2 + CΔt||un − un−1 || L2 ||∇(un − un−1 )|| L2 ||∇
+ CΔt||un − 2un−1 + un−2 ||2L2 || Aun ||2L2
+ CΔt||un−1 − un−2 || L2 ||∇(un−1 − un−2 )|| L2 ||∇(un − un−1 )||2L2 Δt n−1 ||2 2 ) un ||2L2 + CΔt||un−1 ||2L2 (|| Aun−1 ||2L2 + ||U ≤ ||∇ H 4Re + C (Δt)4
tn
t n −1
eun−1 ||2L2 + C (Δt)5 , || Aut (t)||2L2 dt + C (Δt)3 ||∇
and
−2Δt(W n − W n−1 , un ) ≤
Δt ||∇ un ) un ||2L2 − 2SΔt(Bn−1 × curl bn , 4Re + CΔt||bn−1 ||2L2 (||bn−1 ||2H2 + ||Bn−1 ||2H2 )
+ CΔt||bn−1 − bn−2 || L2 ||curl (bn−1 − bn−2 )|| L2 ||curl enb −1 ||2L2 + CΔt||enb −1 ||2L2 ||bn − bn−1 ||2H2 + CΔt||bn−1 ||2L2 ||bn−1 ||2H2 + CΔt||bn − bn−1 || L2 ||curl (bn−1 − bn−2 )||3L2
+ CΔt||bn − 2bn−1 + bn−2 ||2L2 ||bn−1 ||2H2 Δt ≤ ||∇ un ) un ||2L2 − 2SΔt(Bn−1 × curl bn , 4Re + CΔt||bn−1 ||2L2 (||bn−1 ||2H2 + ||Bn−1 ||2H2 ) + C (Δt)5 + C (Δt)3 ||curl ebn−1 ||2L2 + C (Δt)4
tn
t n −1
||bt (t)||2H2 dt,
un , un ) = b(un−1 , un , un ) = 0 and where we use b(eun−1 ,
|| A(un − un−1 )||2L2 ≤ Δt ||bn − bn−1 )||2H2 ≤ Δt
tn
t n −1 tn t n −1
|| Aut (t)||2L2 dt, ||bt (t)||2H2 dt,
4 ||un − 2un−1 + un−2 ||2L2 ≤ (Δt)4 ||utt ||2L2 . 3 18
19
Projection method for MHD equations
Combining the above two estimates into (5.11) yields Δt ||∇ un ||2L2 2Re + (Δt)2 (||∇r np ||2L2 − ||∇r np−1 ||2L2 + ||∇(r np − r np−1 )||2L2 )
||un ||2L2 − ||un−1 ||2L2 + ||un − un−1 ||2L2 +
≤ C (Δt)3 ||∇r np−1 ||2L2 − 2SΔt(Bn−1 × curl bn , un ) + C (Δt)5 + C (Δt)
4
+ C (Δt)4
tn
t n −2 tn t n −1
(||uttt (t)||2L2 + ||∇ ptt (t)||2L2 )dt
(5.12)
(|| Aut (t)||2L2 + ||bt (t)||2H2 )dt
eun−1 ||2L2 + ||curl ebn−1 ||2L2 ) + C (Δt)3 (||∇
n −1 2 || H2 ) ||Bn−1 ||2H2 ).
+ CΔt||un−1 ||2L2 (|| Aun−1 ||2L2 + ||U
+ CΔt||bn−1 ||2L2 (||bn−1 ||2H2 +
Taking the sum of (5.9) and (5.12), we get Δt un ||2L2 ||∇ 2Re SΔt + S||bn ||2L2 − S||bn−1 ||2L2 + S||bn − bn−1 ||2L2 + ||curl bn ||2L2 Rm + (Δt)2 (||∇r np ||2L2 − ||∇r np−1 ||2L2 + ||∇(r np − r np−1 )||2L2 )
||un ||2L2 − ||un−1 ||2L2 + ||un − un−1 ||2L2 +
≤ C (Δt)3 ||∇r np−1 ||2L2 + C (Δt)5 + C (Δt)4 + C (Δt)4
tn
t n −2 tn t n −1
(||uttt (t)||2L2 + ||bttt (t)||2L2 + ||∇ ptt (t)||2L2 )dt (|| Aut (t)||2L2 + ||bt (t)||2H2 )dt
enu−1 ||2L2 + ||curl enb −1 ||2L2 ) + C (Δt)3 (||∇ + CΔt(||un−1 ||2L2 + ||bn−1 ||2L2 )(|| Aun−1 ||2L2 + ||U
+ CΔt||bn−1 ||2L2 (||bn−1 ||2H2 + ||Bn−1 ||2H2 ).
n −1 2 || H2 )
Summing up the above inequality for n from 2 to m ≤ N, using (4.9), (3.14) and the discrete Gronwall inequality, we obtain
||um ||2L2 ≤ C (Δt)4 + C ||u1 || L2 + C ||b1 ||2L2 + C (Δt)2 ||∇r1p ||2L2 = C (Δt)4 + C ||e1u || L2 + C ||e1b ||2L2 + C (Δt)2 ||∇r1p ||2L2 , where we use e0u = e0b = 0. To complete the proof of (5.5), we need to show
||e1u ||2L2 + ||e1b ||2L2 + (Δt)2 ||∇r1p ||2L2 ≤ C (Δt)4 . 19
(5.13)
20
R.An
Let n = 1 in (4.3) and (4.4). It follows from (4.6)-(4.8) and (5.3) that e1u −
Δt 1 e1u − Δt((u1 − u0 ) · ∇)u1 Δ e = Δt∇( P0 − p1 ) − Δt(u0 · ∇) Re u − SΔt(b0 × curl e1b + (b1 − b0 ) × curl b1 ) − ΔtR1u , (5.14)
and e1b +
Δt e1u × b0 ) − ΔtR1b curl (curl e1b ) = Δtcurl ( Rm + Δt((b1 − b0 ) · ∇)u1 − Δt(u1 · ∇)(b1 − b0 ). (5.15)
Testing (5.14) and (5.15) by e1u and Se1b , taking the sum of them, we find e1u ||2L2 + ||e1b ||2L2 ≤ C (Δt)2 (|| R1u ||2L2 + || R1b ||2L2 + ||∇( P0 − p0 )||2L2 + ||∇( p1 − p0 )||2L2 ) ||
+ C (Δt)2 (||∇(u1 − u0 )||2L2 + ||curl (b1 − b0 )||2L2 ) ≤ C (Δt)4 . e1u , then Since e1u = P H
e1u || L2 ≤ C (Δt)2 . ||e1u || L2 ≤ ||
On the other hand, in terms of (3.5), we have Δt||∇r1p || L2 ≤
1 e1u || L2 ) + Δt||∇( p1 − p0 )|| L2 ≤ C (Δt)2 . (||e1u || L2 + || 2
Therefore, we derive (5) and the complete the proof of (5.5).
From (5.5) and the proof of Theorem 4.2, it is easy to show Δt
N
∑ || enp ||2L
2
n =1
≤ C (Δt)2 ,
which implies that (3.15) in Theorem 3.3.
Acknowledgments The authors would like to thank the anonymous reviewers for their careful reviews and valuable comments to improve the quality of the original manuscript. This work is supported by Zhejiang Provincial Natural Science Foundation with Grant No. LY16A010017 and LY14A010020.
References [1] R.Adams, Sobolev Spaces, Academic Press, New York, 1975.
20
21
Projection method for MHD equations
[2] A.Chorin, Numerical solution of the Navier-Stokes equations, Math. Comp., 22(1968), 745-762. [3] J.Douglas Jr., J.Wang, An absolutely stabilized finite element method for the Stokes problem, Math. Comp., 52(1989), 495-508. [4] J.Gerbeau, A stabilized finite elemenet method for the incompressible magnetohydrodynamic equations, Numer. Math., 87(2000), 83-111. [5] J.Gerbeau, C.Le Bris, T.Leli`evre, Mathematical Methods for the Magnetohydrodynamics of Liquid Metals, Oxford: Oxford University Press, 2006. [6] J.Guermond, P.Minev, J.Shen, An overview of projection methods for incompressible flows, Comput. Methods Appl. Mech. Engrg., 195(2006), 6011-6045. [7] M.Gunzburger, A.Meir, J.Peterson, On the existence, uniqueness, and finite element approximation of solutions of the equations of stationary, incompresible magnetohydrodynamics, Math. Comp., 56(1991), 523-563. [8] Y.He, Unconditional convergence of the Euler semi-implicit scheme for the threedimensional incompressible MHD equations, IMA J. Numer. Anal., 35(2015),767-801. [9] J.Heywood, R.Rannacher, Finite-element approximation of the nonstationary NavierStokes problem Part IV: error analysis for second-order time discretization, SIAM J. Numer. Anal., 27(1990), 353-384. [10] W.Hughes, F.Young, The electromagnetics of fluids, Wiley, New York, 1966. [11] R.Moreau, Magneto-hydrodynamics, Kluwer Academic Publishers, 1990. [12] A.Prohl, Convergent finite element discretizations of the nonstationary incompressible magnetohydrodynamic system, ESAIM:M2AN, 42(2008), 1065-1087. ¨ [13] D.Schotzau, Mixed finite element methods for stationary incompressible magnetohydrodynamics, Numer. Math., 96(2004), 771-800. [14] M.Sermange, R.Temam, Some mathematical questions related to the MHD equations, Comm. Pure Appl. Math., 36(1983), 635-664. [15] J.Shen, On error estimates of projection methods for Navier-Stokes equations: first-order schemes, SIAM J. Numer. Anal., 29(1992):57-77. [16] J.Shen, Remarks on the pressure error estimates for the projection methods, Numer.Math., 67(1994), 513-520. [17] R.Temam, Sur l’approximation de la solution des equations de Navier-Stokes par la m´ethode des pas fractionnaires II, Arch. Rational Mech. Anal., 33(1969), 377-385. [18] R.Temam, Navier-Stokes Equations, North-Holland, Amsterdam, 1977. [19] 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. [20] Y.Zhang, Y.Hou, L.Shan, Numerical analysis of the Crank-Nicolson extrapolation time discrete scheme for magnetohydrodynamics flows, Numer. Methods Part. Diff. Equa., 31(2015), 2169-2208.
21
Rong An, Yuan Li
PII: DOI: Reference:
S0168-9274(16)30204-5 http://dx.doi.org/10.1016/j.apnum.2016.10.010 APNUM 3116
To appear in:
Applied Numerical Mathematics
Received date: Revised date: Accepted date:
25 August 2015 20 September 2016 2 October 2016
Please cite this article in press as: R. An, Y. Li, Error analysis of first-order projection method for time-dependent magnetohydrodynamics equations, Appl. Numer. Math. (2016), http://dx.doi.org/10.1016/j.apnum.2016.10.010
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Error analysis of first-order projection method for timedependent magnetohydrodynamics equations Rong An∗
Yuan Li
College of Mathematics and Information Science, Wenzhou University, Wenzhou 325025, P.R.China.
Abstract. This paper focuses on a linearized fully discrete projection scheme for time-dependent magnetohydrodynamics equations in three-dimensional bounded domain. It is shown that the proposed projection scheme allows for a discrete energy inequality and is unconditionally stable. In addition, we present a rigorous analysis for the rates of convergence. AMS subject classifications: 65M12, 76W05 Key words: Magnetohydrodynamics equations, Projection method, Stability, Rate of convergence, Temporal errors.
1
Introduction
The incompressible magnetohydrodynamics (MHD) equations are used to describe the flow of a viscous, incompressible and electrically conducting fluid. For the understanding of the physical background of the MHD equations, we refer to Hughes [10] and Moreau [11]. Let Ω ⊂ R3 be a bounded and simply-connected domain which is either convex or has a C1,1 boundary ∂Ω. We consider the non-dimensional MHD equations in the primitive variable formulation ∂u 1 − Δu + (u · ∇)u + ∇ p + Sb × curl b = f, ∂t Re div u = 0, div b = 0, ∂b 1 + curl (curl b) − curl (u × b) = 0, ∂t Rm
(1.1) (1.2) (1.3)
for x ∈ Ω and t ∈ (0, T ) with T > 0, where 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 equations (1.1)-(1.3) couple the incompressible Navier-Stokes ∗ Corresponding
author. Email: [email protected], [email protected] (Rong An) 1
c
2
R.An
equations with Maxwell’s equations. Thus, the unknowns in (1.1)-(1.3) are the fluid velocity u, the pressure p and the magnetic field b. The MHD equations should be completed by the appropriate initial and boundary conditions. For the sake of simplicity, we consider the following initial and boundary conditions: u( x, 0) = u0 , u = 0,
b · n = 0,
b( x, 0) = b0 curl b × n = 0
in Ω,
(1.4)
on ∂Ω × [0, T ],
(1.5)
where n denotes the unit outward normal vector on ∂Ω. The initial vector functions u0 and b0 satisfy the compatibility condition div u0 = 0 and div b0 = 0. It was observed that testing (1.1) and (1.3) by u and Sb, respectively, and adding the resulting equations leads to the following energy identity: 1 S 1 d |u|2 + S|b|2 dx + |∇u|2 + |curl b|2 dx = f · udx ∀ t > 0. 2 dt Ω Re Rm Ω Ω Thus, for any prescribed initial data (u0 , b0 ) ∈ L2 (Ω) × L2 (Ω), the MHD problem (1.1)-(1.5) exists the global weak solutions. On the regularities of the weak solutions, Sermange-Temam in [14] established the existences of local strong solution with large initial data and global strong solution with small initial data. For the numerical analysis of the MHD problem, the mixed finite element approximations were first proposed and studied for the stationary MHD problem in [7], where H1 - conforming elements were used to discretize the magnetic field provided that Ω is either convex or has a C1,1 boundary ∂Ω. Inspired by the absolutely stable methods for Stokes problem in [3], a stabilized mixed finite element method for stationary MHD problem was developed by Gerbeau [4]. For the time-dependent MHD equations (1.1)-(1.3), recently, He proposed a linearized semi-implicit Euler scheme in [8], where L2 -unconditional convergence of this scheme was proved by using the negative norm technique. Other fully discrete Crank-Nicolson schemes were studied in [19, 20]. 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 corresponding magnetic field. A mixed finite element formulation based on H(curl)¨ elements (or N´ed´elec elements) for b was proposed and studied by Schotzau in [13] for the stationary MHD problem. For a review of various numerical methods for the MHD equations, we refer to Gerbeau-Bris-Leli`evre [5]. In the present work, we will consider the three-dimensional time-dependent MHD equations (1.1)-(1.3) with the initial and boundary conditions (1.4) and (1.5), and propose a linearized fully discrete scheme based on a modified Chorin’s projection scheme for Navier-Stokes equations in [15]. The projection methods were first proposed by Chorin [2] and Temam [17], and have been further developed in various directions. The main advantage of the projection methods 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. 2
3
Projection method for MHD equations
For a review of projection methods for the Navier-Stokes equations, we refer to [6]. To state our main results derived in this paper, we recall the following terminology [15]: Definition 1.1. Let X be a Banach space equipped with norm || · || X and f : [0, T ] −→ X is continuous. Let 0 = t0 < t1 < · · · < t N = T be a uniform partition of the time interval [0, T ] with time step Δt = T/N and tn = nΔt for 0 ≤ n ≤ N. We say f Δt is a weakly order α approximation of f in X if there exists a constant C independent of Δt such that Δt
N
∑ || f Δt (tn ) − f (tn )||2X ≤ C(Δt)2α ;
n =0
and we say f Δt is a strongly order α approximation of f in X if there exists a constant C independent of Δt and n such that max || f Δt (tn ) − f (tn )||2X ≤ C (Δt)2α .
0≤ n ≤ N
To our best knowledge, Prohl in [12] first proposed a projection scheme for the MHD problem. Unfortunately, the projection scheme in [12] does not allow for a discrete energy estimate. Moreover, it was proved that this scheme provided the weakly order 12 approximations of the velocity field and the magnetic field in H1 (Ω). In this paper, we propose a new linearized projection scheme, which allows for a discrete energy inequality, thus, and is unconditionally stable. We will prove that this new projection scheme provides the weakly first-order approximations of the velocity field and the magnetic field in H1 (Ω), and the strongly first-order approximations of the velocity field and the magnetic field in L2 (Ω) under the regularity assumptions 1 and 2 below. To obtain the optimal convergence order of the pressure, we need some further regularity assumption as like in [16], under which we can derive the weakly first-order approximation of the pressure in L2 (Ω). This paper is organized as follows: in the next section, we begin with some notation and lay out some regularity 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 inequality is derived in Section 3. The proofs of main results are given in Sections 4 and 5, which are split into several theorems.
2 Mathematical setting For the mathematical setting of the MHD equations (1.1)-(1.3) with the initial and boundary conditions (1.4)-(1.5), we introduce some function spaces and their associated norms. For k ∈ N + , 1 ≤ p ≤ +∞, let W k,p (Ω) denote the standard Sobolev 3
4
R.An
space. The norm in W k,p (Ω) is defined by
||u||W k,p =
⎧ ⎪ ⎪ ⎨ ( ⎪ ⎪ ⎩
∑
| β|≤k Ω
∑
|∇ β u| p dx )1/p
for 1 ≤ p < +∞,
sup |∇ β u|
for p = +∞,
| β|≤k Ω
where ∇ is the differential operator with respect to x and
∇β =
∂| β| β
β
β
∂x1 1 ∂x2 2 ∂x3 3
for the multi-index β = ( β 1 , β 2 , β 3 ) and | β| = β 1 + β 2 + β 3 with β 1 , β 2 , β 3 ≥ 0. We k,p define W0 (Ω) to be the subspace of W k,p (Ω) of functions with zero trace on ∂Ω. When p = 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,p (Ω) and L p (Ω) are used to denote the vector Sobolev spaces H k (Ω)3 , W k,p (Ω)3 and L p (Ω)3 , respectively. In particular, (·, ·) denotes the L2 (Ω) inner product. Let X be a Banach space. For some T > 0 or T = ∞, denote Bochner spaces L p (0, T; X ), 1 ≤ p < +∞ the spaces of measurable functions from the interval [0, T ] into X such that T 0
p
||u(t)|| X dt < +∞.
If p = +∞, the functions in L∞ (0, T; X ) are required to satisfy sup ||u(t)|| X < +∞.
t∈[0,T ]
The symbols C, C1 , C2 , · · · are used to denote the generic positive constants independent of the time step Δt. We define H = {u ∈ L2 (Ω), div u = 0 in Ω, u · n = 0 on ∂Ω}, V = H10 (Ω),
V0 = {u ∈ V, div u = 0 in Ω},
X = {u ∈ H1 (Ω), u · n = 0 on ∂Ω}, M=
L20 (Ω)
X0 = {u ∈ X, div u = 0 in Ω},
= { q ∈ L ( Ω ), 2
Ω
qdx = 0}.
Thanks to Poincar`e inequality, the norm in V is defined by
||u||V =
|∇u| dx 2
Ω
4
1/2 .
5
Projection method for MHD equations
Define the bilinear forms
a2 (u, v) =
1 Rm
1 a1 (u, v) = Re
Ω
∇u : ∇vdx
∀ u, v ∈ V,
(curl u · curl v + div udiv v)dx
∀ u, v ∈ X,
Ω
d(v, p) =
Ω
∀ v ∈ V, p ∈ M
pdiv vdx
and the trilinear form
1 (div u)v · wdx b(u, v, w) = (u · ∇)v · wdx + 2 Ω Ω 1 1 = (u · ∇)v · wdx − (u · ∇)w · vdx 2 Ω 2 Ω
∀ u, v, w ∈ V.
It is obvious that b(u, v, w) satisfies b(u, v, w) = −b(u, w, v)
∀ u, v, w ∈ V.
(2.1)
With the above notations, given u0 , b0 ∈ H and f ∈ L∞ (0, T; L2 (Ω)), the weak variational formulation of the MHD problem (1.1)-(1.5) is defined to find (u, p, b) ∈ L2 (0, T; V) × L2 (0, T; M) × L2 (0, T; X) such that
(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) − (u × b, curl w) = 0
∀ w ∈ X,
(2.3)
where we use the following formula: Ω
curl u · vdx =
Ω
u · curl vdx +
∂Ω
n × u · vds.
(2.4)
Corresponding to (1.1), we recall Stokes operator A [18]. Let us denote the orthogonal projection operator from L2 (Ω) onto H by PH . Then Stokes operator A is defined by Au = −PH Δu ∀ u ∈ D( A) = V0 ∩ H2 (Ω). The emphasis of this paper is to prove the temporal error estimates of the projection scheme (3.1)-(3.3) stated in Section 3. Then we make the following assumptions on the initial data and the regularities of the solution. Assumption 1: The initial data u0 , b0 and f satisfy
|| Au0 || L2 + ||b0 || H2 + sup {||f(t)|| L2 + ||ft (t)|| L2 } ≤ C. 0≤ t ≤ T
5
(2.5)
6
R.An
Assumption 2: The MHD problem (1.1)-(1.5) admits a unique local strong solution (u, p, b) on (0, T ) such that sup {|| Au(t)|| L2 + ||b(t)|| H2 + || p(t)|| H1 + ||ut (t)|| L2 + ||bt (t)|| L2 }
0≤ t ≤ T
+
T 0
||ut ||2H1 + ||bt ||2H1 + || pt ||2H1 + ||utt ||2L2 + ||btt ||2L2 dt ≤ C. (2.6)
The following known inequalities are frequently used [1, 18]:
||PH v|| H1 ≤ C ||v|| H1 ∀ v ∈ H1 (Ω), ||v|| H1 ≤ C (||curl v|| L2 + ||div v|| L2 ) ∀ v ∈ X, ||v|| Lr ≤ C ||v|| H1 (2 ≤ r ≤ 6) ∀ v ∈ H (Ω), 1
||v|| L3 ≤ ||v|| L∞ ≤
C ||v||1/2 ||v||1/2 L2 H1 1/2 1/2 C ||v|| H1 ||v|| H2
||v|| H2 ≤ C || Av|| L2
(2.7) (2.8) (2.9)
∀ v ∈ H ( Ω ),
(2.10)
∀ v ∈ H (Ω), i = 1, 2,
(2.11)
1 2
∀ v ∈ D( A ).
(2.12)
Finally, we recall a discrete version of Gronwall inequality established in [9]. Lemma 2.1. Let ak , bk , ck and γk , for integers k ≥ 0, be the nonnegative numbers such that an + τ
n
n
n
k =0
k =0
k =0
∑ bk ≤ τ ∑ γk a k + τ ∑ c k + B
for n ≥ 0.
(2.13)
Suppose that τγk < 1, for all k, and set σk = (1 − τγk )−1 . Then an + τ
n
∑
k =0
bk ≤ exp(τ
n
∑
k =0
γk σk )(τ
n
∑ ck + B)
k =0
for n ≥ 0.
(2.14)
Remark 2.1. If the first sum on the right in (2.13) extends only up to n − 1, then estimate (2.14) holds for all k > 0 with σk = 1.
3
A projection scheme and main results
Let 0 = t0 < t1 < · · · < t N = T be a uniform partition of the time interval [0, T ] with time step Δt = T/N and tn = nΔt, 0 ≤ n ≤ N, where [0, T ] is the maximal time interval such that a local strong solution exists. For 0 ≤ n ≤ N, let us denote un = u(tn ), pn = p(tn ), bn = b(tn ) and fn = f(tn ). For any sequence of functions { gn }nN=0 , we define Dt gn+1 = ( gn+1 − gn )/Δt for 0 ≤ n ≤ N − 1. We start with 0 = u0 , B0 = b0 and an arbitrary P0 ∈ M ∩ H 1 (Ω). We solve (U n , Bn ) and U0 = U (Un , Pn ) for 1 ≤ n ≤ N, respectively, by n − Un −1 U 1 n n + ∇ Pn−1 + SBn−1 × curl Bn = fn , − ΔU + (Un−1 · ∇)U Δt Re 1 n n × Bn−1 ) = 0, div Bn = 0 Dt B + curl (curl Bn ) − curl (U Rm 6
(3.1) (3.2)
7
Projection method for MHD equations
n = 0, Bn · n = 0 and curl Bn × n = 0 on ∂Ω, and with the boundary condition U n
Un − U + 2∇( Pn − Pn−1 ) = 0, Δt
div Un = 0,
(3.3)
with the boundary condition Un · n = 0 on ∂Ω. A reformulation of the problem (3.3) is 1 n div U Δt
2ΔPn = 2ΔPn−1 +
(3.4)
with the boundary condition ∇( Pn − Pn−1 ) · n = 0 on ∂Ω, and n − 2Δt∇( Pn − Pn−1 ). Un = U
(3.5)
Since the problems (3.1)-(3.3) consist of two linear equations, its solvability should follow easily if we can establish a priori energy estimates for its solutions. The following discrete energy identity just serves this purpose. n , Bn ) and (Un , Pn ) be the solution of (3.1)-(3.3). Then Theorem 3.1. For 1 ≤ n ≤ N, let (U the following energy identity holds:
||Um ||2L2 + S||Bm ||2L2 + 2(Δt)2 ||∇ Pm ||2L2 + Δt +
m
∑ (||U
n
n =1
m
2
∑ ( Re ||∇U
n =1
n
|| L2 +
2S ||curl Bn ||2L2 ) Rm
1 n − Un−1 ||2L2 + ||U − Un ||2L2 + S||Bn − Bn−1 ||2L2 ) 2
= ||u0 ||2L2 + S||b0 ||2L2 +2(Δt)2 ||∇ P0 ||2L2 + 2Δt
m
∑ (fn , U
n =1
n
) (3.6)
for every 1 ≤ m ≤ N. n
and 2SΔtBn , respectively, and adding the Proof. Testing (3.1) and (3.2) by 2ΔtU resulting equations yield n ||2 2 + ||U n − Un−1 ||2 2 − ||Un−1 ||2 2 + S(||Bn ||2 2 − ||Bn−1 ||2 2 + ||Bn − Bn−1 ||2 2 ) ||U L L L L L L 2Δt 2SΔt n n n n n 2 n − 1 ) = 2Δt(f , U ). (3.7) || L2 + ||∇U ||curl B || L2 + 2Δt(∇ P , U + Re Rm Here we use the following identity:
(u × curl v, w) = (w × u, curl v) Testing (3.5) by Un and lead to
∀ u, v, w ∈ H1 (Ω).
(3.8)
1 n n (U + U ), respectively, and adding the resulting equations 2
n ||2 2 + 1 ||Un − U n ||2 2 = Δt(∇( Pn−1 − Pn ), U n ). ||Un ||2L2 − ||U L L 2 7
(3.9)
8
R.An
The sum of (3.7) and (3.9) yields n − Un−1 ||2 2 − ||Un−1 ||2 2 + 1 ||Un − U n ||2 2 ||Un ||2L2 + ||U L L L 2 2Δt n || L2 + 2SΔt ||curl Bn ||2 2 + S(||Bn ||2L2 − ||Bn−1 ||2L2 + ||Bn − Bn−1 ||2L2 ) + ||∇U L Re Rm n n n n −1 n = 2Δt(f , U ) − Δt(∇( P + P ), U ). (3.10) n = Un + 2Δt∇( Pn − Pn−1 ) and div Un = 0, we have Observing U n ) = −2(Δt)2 (||∇ Pn ||2 2 − ||∇ Pn−1 ||2 2 ). −Δt(∇( Pn + Pn−1 ), U L L Putting the above equation into (3.10) and taking the sum for n from 1 to m ≤ N, we complete the proof of (3.6). As a consequence of (3.6), we have the following result. Corollary 3.1. The projection scheme (3.1)-(3.3) is unconditionally stable. Moreover, there exists some C > 0 such that
||Um ||2L2 + S||Bm ||2L2 + 2(Δt)2 ||∇ Pm ||2L2 + Δt m
m
1
∑ ( Re ||∇U
n =1
n
|| L2 +
S ||curl Bn ||2L2 ) Rm
− Un−1 ||2 2 + 1 ||U n − Un ||2 2 + S||Bn − Bn−1 ||2 2 ) ≤ C (3.11) + ∑ (||U L L L 2 n =1 n
for every 1 ≤ m ≤ N. Let us denote enu = un − Un ,
n, enu = un − U
enb = bn − Bn ,
P n = Pn−1 + 2( Pn − Pn−1 ),
enp = pn − Pn ,
enp = pn − P n .
The main results derived in this paper are presented in the following two theorems. Theorem 3.2. Suppose that Assumptions 1-2 are satisfied. Then there exists some positive constant C > 0 such that sup
||um − Um ||2L2 + ||bm − Bm ||2L2 ≤ C (Δt)2 ,
(3.12)
+ ||bn − Bn ||2H1 + Δt|| pn − P n ||2L2 ) ≤ C (Δt)2 .
(3.13)
1≤ m ≤ N
Δt
N
∑ (||∇(un − Un )||2L
n =1
2
8
9
Projection method for MHD equations
From (3.13), we can see that the weakly order 12 approximation of the pressure in L2 (Ω) is derived. However, this convergence order is not optimal. To improve it, we need some further regularity assumptions as like in [16]. Now, we assume that sup (||∇ut (t)|| L2 + ||curl bt (t)|| L2 + ||utt (t)|| L2 + ||btt (t)|| L2 + || pt (t)|| H1 )
0≤ t ≤ T
+
T 0
(|| Aut ||2L2 + ||bt ||2H2 + ||uttt ||2L2 + ||bttt ||2L2 + || ptt ||2H1 )dt ≤ C. (3.14)
Then we can prove the wekly first-order approximation of the pressure in L2 (Ω). Theorem 3.3. Suppose that Assumptions 1-2 and (3.14) are satisfied. If
||∇( P0 − p0 )|| L2 ≤ CΔt, then there exists some C > 0 such that Δt
4
N
∑ || pn − P n ||2L
n =1
2
≤ C (Δt)2 .
(3.15)
Proof of Theorem 3.2
For 1 ≤ n ≤ N, taking t = tn in (1.1)-(1.3), we obtain 1 un − un −1 − Δun + (un · ∇)un + ∇ pn + Sbn × curl bn = fn − Rnu , Δt Re
(4.1)
and bn − bn −1 1 − curl (curl bn ) − curl (un × bn ) = − Rnb , Δt Rm
(4.2)
where Rnu =
1 Δt
tn t n −1
(s − tn−1 )∂tt u(s)ds,
Rnb =
1 Δt
tn t n −1
(s − tn−1 )∂tt b(s)ds.
Subtracting (3.1) and (3.2) from (4.1) and (4.2), respectively, leads to 1 enu − enu−1 n − (un · ∇)un + ∇( Pn−1 − pn ) − Δ en = (Un−1 · ∇)U Δt Re u + S(Bn−1 × curl Bn − bn × curl bn ) − Rnu , (4.3) and enb − enb −1 1 n × Bn −1 ) − R n , + curl (curl enb ) = curl (un × bn ) − curl (U b Δt Rm 9
(4.4)
10
R.An
An alternative to (3.5) is enu = 2Δt∇( Pn − Pn−1 ). enu −
(4.5)
n
− (un · ∇)un and S(Bn−1 × curl Bn − bn × curl bn ) in (4.3) can The terms (Un−1 · ∇)U be rewritten as n − (un · ∇)un (Un−1 · ∇)U
(4.6)
enu − (eun−1 · ∇)un − (un−1 · ∇) enu − ((un − un−1 ) · ∇)un = (enu−1 · ∇) and S(Bn−1 × curl Bn − bn × curl bn )
(4.7)
= −S(Bn−1 × curl enb + ebn−1 × curl bn + (bn − bn−1 ) × curl bn ). n
× Bn−1 ) in (4.4) can be rewritten as The term curl (un × bn ) − curl (U n × Bn −1 ) curl (un × bn ) − curl (U
(4.8)
enu × Bn−1 ) + curl (un × (bn − bn−1 )) + curl (un × ebn−1 ). = curl (
Theorem 4.1. Suppose that Assumptions 1-2 are satisfied. Then there exists some positive constant C1 > 0 such that 2 m 2 2 m 2 ||em u || L2 + ||eb || L2 + ( Δt ) ||∇ e p || L2 + Δt
+ Proof. Testing (4.3) by
2Δt enu
m
∑ (||∇ enu ||2L
n =1
m
∑ (|| enu − enu−1 ||2L
n =1
2
2
+ ||∇enu ||2L2 + ||curl enb ||2L2 )
+ ||enu − enu ||2L2 ) ≤ C1 (Δt)2 . (4.9)
and using (u − v, 2u) = |u|2 − |v|2 + |u − v|2 , we have
2Δt enu ||2L2 ||∇ Re ≤ 2Δt|b(un − un−1 , un , enu ) + b(enu−1 , un , enu )| + 2Δt|( Rnu , enu )| enu ||2L2 − ||eun−1 ||2L2 + || enu − eun−1 ||2L2 + || n −1
n
+ 2SΔt(B × curl B − b × curl b = I1 + I2 + I3 + I4 . n
n
, enu )
+ 2Δt(∇( P
n −1
(4.10)
−p
n
), enu )
From (2.6) and the Young’s inequality, we bound I1 as enu || L2 I1 ≤ CΔt(||un − un−1 || L2 + ||enu−1 || L2 )|| Aun || L2 ||∇ Δt enu ||2L2 + CΔt(||un − un−1 ||2L2 + ||eun−1 ||2L2 ). ≤ ||∇ 4Re It is easy to show that enu || L2 ≤ I2 ≤ CΔt|| Rnu || L2 ||∇
Δt enu ||2L2 + CΔt|| Rnu ||2L2 . ||∇ 4Re
10
11
Projection method for MHD equations
¨ inequality, Poincar´e inequality and the From (2.9), (4.7), H1 (Ω) → L6 (Ω), Holder Young’s inequality, I3 is bounded by I3 = −2SΔt(Bn−1 × curl enb , enu ) − 2SΔt(enb −1 × curl bn , enu )
− 2SΔt((bn − bn−1 ) × curl bn , enu )
≤ −2SΔt(Bn−1 × curl enb , enu ) + CΔt(||enb −1 || L2 + ||bn − bn−1 || L2 )||curl bn || L3 ||∇ enu || L2 Δt enu ||2L2 − 2SΔt(Bn−1 × curl enb , ≤ ||∇ enu ) 4Re + CΔt(||enb −1 ||2L2 + ||bn − bn−1 ||2L2 ). Combining these estimates into (4.10) leads to Δt enu ||2L2 ||∇ Re ≤ CΔt(||ebn−1 ||2L2 + ||bn − bn−1 ||2L2 ) − 2SΔt(Bn−1 × curl enb , enu ) enu ||2L2 − ||enu−1 ||2L2 + || enu − enu−1 ||2L2 + ||
(4.11)
+ CΔt(||un − un−1 ||2L2 + ||eun−1 ||2L2 + || Run ||2L2 ) + 2Δt(∇( Pn−1 − pn ), enu ).
1 Testing (4.5) by enu and (enu + enu ), respectively, we obtain 2 1 enu ||2L2 + ||enu − enu ||2L2 ) = 0, (||enu ||2L2 − || 2 1 enu ||2L2 ) = Δt(∇( Pn − Pn−1 ), (||enu ||2L2 − || enu ), 2 where we use div enu = 0 in Ω and enu = 0 on ∂Ω. Taking the sum of the above identities leads to 1 enu ||2L2 + ||enu − enu ||2L2 = Δt(∇( Pn − Pn−1 ), ||enu ||2L2 − || enu ), 2 which together with (4.11) implies that 1 Δt enu − enu−1 ||2L2 + ||enu − enu ||2L2 + enu ||2L2 ||enu ||2L2 − ||enu−1 ||2L2 + || ||∇ 2 Re ≤ CΔt(||enb −1 ||2L2 + ||bn − bn−1 ||2L2 ) − 2SΔt(Bn−1 × curl enb , enu )
+ CΔt(||un − un−1 ||2L2 + ||eun−1 ||2L2 + || Rnu ||2L2 )
− Δt(∇(enp + enp−1 + pn − pn−1 ), enu ). From (4.5) again, one has
enu = enu + 2Δt∇( Pn−1 − Pn )
= enu + 2Δt∇(enp − enp−1 − ( pn − pn−1 )). 11
(4.12)
12
R.An
Then we have
− Δt(∇(enp + enp−1 + pn − pn−1 ), enu ) = −2(Δt)2 (∇(enp + enp−1 + pn − pn−1 ), ∇(enp − enp−1 − ( pn − pn−1 )) = 2(Δt)2 (||∇enp−1 ||2L2 − ||∇enp ||2L2 ) + 2(Δt)2 ||∇( pn − pn−1 )||2L2 + 4(Δt)2 (∇enp−1 , ∇( pn − pn−1 )) ≤ 2(Δt)2 (||∇enp−1 ||2L2 − ||∇enp ||2L2 ) + 2(Δt)2 ||∇( pn − pn−1 )||2L2 + 2Δt||∇( pn − pn−1 )||2L2 + 2(Δt)3 ||∇enp−1 ||2L2 . Substituting the above inequality into (4.12) we obtain
||enu ||2L2 − ||enu−1 ||2L2 + 2(Δt)2 (||∇enp ||2L2 − ||∇enp−1 ||2L2 ) 1 Δt + || enu − enu−1 ||2L2 + ||enu − enu ||2L2 + enu ||2L2 ||∇ 2 Re ≤ CΔt(||ebn−1 ||2L2 + ||bn − bn−1 ||2L2 ) − 2SΔt(Bn−1 × curl enb , enu ) + CΔt(||un − un−1 ||2L2 + ||enu−1 ||2L2 + || Rnu ||2L2 )
+ 4Δt||∇( pn − pn−1 )||2L2 + 2(Δt)3 ||∇enp−1 ||2L2 . Using (4.8) and testing (4.4) by 2SΔtenb leads to 2SΔt ||curl enb ||2L2 Rm enu × Bn−1 , curl enb ) + (un × (bn − bn−1 ), curl enb ) = 2SΔt(( S||enb ||2L2 − S||ebn−1 ||2L2 + S||enb − enb −1 ||2L2 +
+ (un × enb −1 , curl enb )) − 2SΔt( Rbn , enb )
≤ 2SΔt( enu × Bn−1 , curl enb ) + CSΔt|| Rbn || L2 ||curl enb ||2L2
+ CSΔt(||bn − bn−1 || L2 + ||enb −1 || L2 )||un || H2 ||curl enb || L2 SΔt enu × Bn−1 , curl enb ) ≤ ||curl enb ||2L2 + 2SΔt( Rm + CΔt(||bn − bn−1 ||2L2 + ||enb −1 ||2L2 + || Rnb ||2L2 ). From (3.8), the sum of the above two inequalities yields
||enu ||2L2 − ||eun−1 ||2L2 + S||enb ||2L2 − S||ebn−1 ||2L2 + 2(Δt)2 (||∇enp ||2L2 − ||∇enp−1 ||2L2 ) 1 Δt 2SΔt enu − enu−1 ||2L2 + ||enu − enu ||2L2 + enu ||2L2 + ||∇ ||curl enb ||2L2 + || 2 Re Rm ≤ CΔt(||ebn−1 ||2L2 + ||enu−1 ||2L2 + (Δt)2 ||∇enp−1 ||2L2 ) + CΔt(|| Rnu ||2L2 + || Rnb ||2L2 ) + CΔt(||un − un−1 ||2L2 + ||bn − bn−1 ||2L2 + ||∇( pn − pn−1 )||2L2 ). Taking the sum of the above inequality for n from 1 to m ≤ N and using the discrete Gronwall inequality (2.14), we conclude that there exists some positive constant C1 > 12
13
Projection method for MHD equations
0 such that 2 m 2 2 m 2 ||em u || L2 + ||eb || L2 + ( Δt ) ||∇ e p || L2 + Δt
+
m
∑ (||∇ enu ||2L
n =1
m
∑ (|| enu − enu−1 ||2L
n =1
2
2
+ ||curl enb ||2L2 )
enu ||2L2 ) ≤ C1 (Δt)2 (4.13) + ||enu −
by using N
∑ (||un − un−1 ||2L
n =1
2
+ ||bn − bn−1 ||2L2 + ||∇( pn − pn−1 )||2L2 ) ≤ CΔt, N
∑ (|| Rnu ||2L
2
n =1
+ || Rnb ||2L2 ) ≤ CΔt.
(4.14) (4.15)
Applying the projector PH to (4.5) and using (2.7), (4.13) implies the inequality (4.9) and completes the proof of Theorem 4.1. Theorem 4.1 implies that the projection scheme (3.1)-(3.3) provides the weakly and strongly first-order approximations of the velocity field and the magnetic field in H1 (Ω) and L2 (Ω), respectively. Based on the above theorem, the rate of convergence for the pressure is stated in the following theorem. Theorem 4.2. Suppose that Assumptions 1-2 are satisfied. Then there exists some positive constant C2 > 0 such that Δt
N
∑ || enp ||2L
n =1
2
≤ C2 Δt.
(4.16)
Proof. Taking the sum of (3.1) and (3.3) gives Un − Un −1 1 n n + ∇ P n + SBn−1 × curl Bn = fn . − ΔU + (Un−1 · ∇)U Δt Re Subtracting (4.17) from (4.1) yields
(4.17)
1 enu − enu−1 n − (un · ∇)un − ∇ enp − Δ en = (Un−1 · ∇)U Δt Re u + SBn−1 × curl Bn − Sbn × curl bn − Rnu . (4.18) From (4.18) and inf-sup condition [18], one has enp || L2 ≤ C sup ||
v∈V,v =0
≤
d(v, enp )
||∇v|| L2
C n enu || L2 + C || Rnu || L2 ||e − enu−1 || L2 + C ||∇ Δt u enu || L2 + ||∇eun−1 || L2 + ||∇ enu || L2 + ||un − un−1 || L2 ) + C (||∇enu−1 || L2 ||∇
+ C (||ebn−1 || H1 ||curl enb || L2 + ||curl enb || L2 + ||ebn−1 || L2 + ||bn − bn−1 || L2 ). 13
14
R.An
Taking the sum of the above inequality from 1 to N and using (3.11), (4.9), (4.14) and (4.15), we have Δt
N
∑ || enp ||2L
2
n =1
≤ CΔt + C (Δt)2
N
∑ (||enb −1 ||2H
1
n =1
+ ||∇eun−1 ||2L2 ),
which together with (3.11) and (4.9) completes the proof of Theorem 4.2.
5
Proof of Theorem 3.3
As a consequence of (4.9) and Assumption 2, we have
|| Dt Un || L2 + || Dt Bn || L2 ≤ C
for 1 ≤ n ≤ N.
(5.1)
n , Un , Bn ) can be derived in next lemma. Then, the regularities of (U Lemma 5.1. Suppose that Assumptions 1-2 are satisfied. Then there exists some τ0 > 0 such that when Δt < τ0 , there holds Δt
N
∑ (||Un ||2H
n =1
2
n ||2 2 + ||Bn ||2 2 ) ≤ C. + ||U H H
(5.2)
Proof. From (3.1) and the regularity of elliptic problem, one has n
|| H2 ≤ C (|| Dt Un || L2 + (Δt)−1 || enu − enu || L2 + ||fn || L2 ||U n
|| L2 ). + ||∇ Pn−1 || L2 + ||Bn−1 × curl Bn || L2 + ||(Un−1 · ∇)U The last three terms in the above inequality can be bounded, respectively, by
||∇ Pn−1 || L2 ≤ ||∇enp−1 || L2 + ||∇ pn−1 || L2 and
||Bn−1 × curl Bn || L2 ≤ ||enb −1 × curl Bn || L2 + ||bn−1 × curl Bn || L2 ≤ C (||ebn−1 || H1 ||Bn || H2 + ||Bn || H1 )
and n
n
n
|| L2 ≤ ||(eun−1 · ∇)U || L2 + ||(un−1 · ∇)U || L2 ||(Un−1 · ∇)U n n || H2 + ||∇U || L2 ). ≤ C (||∇eun−1 || L2 ||U Similarly, from (3.2) and the following formulation curl (u × v) = (div v)u − (u · ∇)v + (v · ∇)u − (div u)v, 14
(5.3)
15
Projection method for MHD equations
we have n · ∇)Bn−1 || L2 + ||(Bn−1 · ∇)U n || L2 + ||(div U n )Bn−1 || L2 ) ||Bn || H2 ≤ C (|| Dt Bn || L2 + ||(U n || H2 + ||∇U n || L2 ). ≤ C (|| Dt Bn || L2 + ||en−1 || H1 ||U b
Taking the sum of the above inequalities leads to N
∑ (||U
Δt
+ || P n ||2H1 + ||Bn ||2H2 )
∑ (|| Dt Un ||2L
2
∑ (||∇enp−1 ||2L
2
≤ Δt +
n =1 N
n 2 || H2
n =1 N
n =1
+ Δt
enu − enu ||2L2 ) + || Dt Bn ||2L2 + (Δt)−2 || n ||2 2 ) + ||∇ pn−1 ||2L2 + ||fn ||2L2 + ||Bn ||2H1 + ||∇U L
N
∑ (||ebn−1 ||2H ||U 1
n =1
≤ C + C3 (Δt)2
N
∑ (||U
n =1
n 2 || H2
n 2 || H2
n
||2 2 ) + ||enb −1 ||2H1 ||Bn ||2H2 + ||∇eun−1 ||2L2 ||U H
+ ||Bn ||2H2 ).
Let τ0 < 1/2C3 . Then when Δt < τ0 , we complete the proof of Lemma 4.2.
Let us denote enu − eun−1 , r np = enp − enp−1 , bn = enb − enb −1 , un = enu − eun−1 , un =
∀ 1 ≤ n ≤ N.
Then (4.9) implies
||un ||2L2 + ||bn ||2L2 + Δt
m
∑ (||∇un ||2L
n =1
2
un ||2L2 + ||curl bn ||2L2 ) ≤ C (Δt)2 . + ||∇
(5.4)
Lemma 5.2. Suppose that Assumptions 1-2 and (3.14) are satisfied. If
||∇( P0 − p0 )|| L2 ≤ CΔt, then there exists some C4 > 0 such that
||enu − enu−1 || L2 ≤ C4 (Δt)2 ,
∀ 1 ≤ n ≤ N.
(5.5)
Proof. Taking the difference of (4.4) with two consecutive indices, we have bn − bn−1 1 + curl (curl bn ) = T n − T n−1 − ( Rnb − Rbn−1 ), Δt Rm 15
(5.6)
16
R.An
where enu × Bn−1 ) + curl (un × (bn − bn−1 )) + curl (un × ebn−1 ). T n = curl ( By a simple calculation, we find that T n − T n−1 in (5.6) is equal to enu−1 × (bn−1 − bn−2 )) − curl ( eun−1 × bn−1 ) T n − T n−1 = curl ( n × Bn−1 ) + curl (
+ curl ((un − un−1 ) × (bn − bn−1 )) + curl (un−1 × (bn − 2bn−1 + bn−2 )) + curl ((un − un−1 ) × enb −1 ) + curl (un−1 × bn−1 ). (5.7)
Testing (5.6) by 2SΔtbn yields 2SΔt ||curl bn ||2L2 Rm = 2SΔt(T n − T n−1 , bn ) − 2SΔt( Rnb − Rbn−1 , bn ). S||bn ||2L2 − S||bn−1 ||2L2 + S||bn − bn−1 ||2L2 +
(5.8)
From (5.7), (3.14) and the Young’s inequality, we get 2SΔt(T n − T n−1 , bn ) un × Bn−1 , curl bn ) + 2SΔt( enu−1 × (bn−1 − bn−2 ), curl bn ) = 2SΔt( enu−1 × bn−1 , curl bn ) + 2SΔt((un − un−1 ) × (bn − bn−1 ), curl bn ) − 2SΔt(
+ 2SΔt(un−1 × (bn − 2bn−1 + bn−2 ), curl bn )
+ 2SΔt((un − un−1 ) × ebn−1 , curl bn ) + 2SΔt(un−1 × bn−1 , curl bn ) SΔt ≤ 2SΔt( ||curl bn ||2L2 un × Bn−1 , curl bn ) + 2Rm eun−1 ||2L2 + CΔt||bn−1 − bn−2 || L2 ||curl (bn−1 − bn−2 )|| L2 ||∇ n −1 2 || H2 )||bn−1 ||2L2 + CΔt|| Aun−1 ||2L2 ||bn + CΔt||un − un−1 || L2 ||∇(un − un−1 )|| L2 ||curl (bn − bn−1 )||2L2 + CΔt||un − un−1 || L2 ||∇(un − un−1 )|| L2 ||curl ebn−1 ||2L2
+ CΔt(|| Aun−1 ||2L2 + ||U
un × Bn−1 , curl bn ) + ≤ 2SΔt(
SΔt ||curl bn ||2L2 + C (Δt)5 2Rm
enu−1 ||2L2 + ||curl enb −1 ||2L2 ) + CΔt(|| Aun−1 ||2L2 + ||U + C (Δt)3 (||∇
− 2bn−1 + bn−2 ||2L2
n −1 2 || H2 )||bn−1 ||2L2 ,
where we use
||bn − 2bn−1 + bn−2 ||2L2 = ||
tn t n −1
(tn+1 − t)btt (t)dt −
4 ≤ (Δt)4 ||btt ||2L2 . 3 16
t n −1 t n −2
(tn+1 − t)btt (t)dt||2L2
17
Projection method for MHD equations
In addition, by Taylor formula with the integral residue, we have
|| Rnb − Rbn−1 ||2L2 = ||
1 n (b − 2bn−1 + bn−2 ) − (bnt − btn−1 )||2L2 Δt
≤ 2(Δt)
3
tn
t n −2
||bttt (t)||2L2 dt.
Then 2SΔt( Rnb − Rnb −1 , bn ) ≤
SΔt ||curl bn ||2L2 + C (Δt)4 2Rm
tn t n −2
||bttt (t)||2L2 dt.
Combining these estimates into (5.8) leads to S||bn ||2L2 − S||bn−1 ||2L2 + S||bn − bn−1 ||2L2 +
SΔt ||curl bn ||2L2 Rm
un × Bn−1 , curl bn ) + C (Δt)5 + C (Δt)4 ≤ 2SΔt( enu−1 ||2L2 + ||curl enb −1 ||2L2 ) + C (Δt)3 (||∇
tn
t n −2
||bttt (t)||2L2 dt
(5.9)
n−1 ||2 2 )||n−1 ||2 2 . + CΔt(|| Aun−1 ||2L2 + ||U H L b
Similarly, taking the difference of (4.18) with two consecutive indices, we have 1 un − un−1 − Δ n = −∇r np−1 − 2∇(r np − r np−1 ) Δt Re u + ∇( pn − 2pn−1 + pn−2 ) + (X n − X n−1 ) − (W n − W n−1 ) − ( Rnu − Run−1 ), (5.10) where enu − (eun−1 · ∇)un − (un−1 · ∇) enu − ((un − un−1 ) · ∇)un , X n = (enu−1 · ∇) W n = S(Bn−1 × curl enb + ebn−1 × curl bn + (bn − bn−1 ) × curl bn ). Thanks to the proof of Lemma A2 in [16], we find Δt un ||2L2 ||∇ Re + (Δt)2 (||∇r np ||2L2 − ||∇r np−1 ||2L2 + ||∇(r np − r np−1 )||2L2 )
||un ||2L2 − ||un−1 ||2L2 + ||un − un−1 ||2L2 +
≤ 2Δt(X n − X n−1 , un ) − 2Δt(W n − W n−1 , un ) + C (Δt)3 ||∇r np−1 ||2L2 + C (Δt)4
tn
t n −2
(5.11)
(||uttt (t)||2L2 + ||∇ ptt (t)||2L2 )dt.
Since X n − X n−1 = (enu−1 · ∇) n + (un−1 · ∇) eun−1 − (eun−1 · ∇)(un − un−1 ) un − ((un−1 − un−2 ) · ∇) − (un−1 · ∇)un−1 − (un−1 · ∇) eun−1
− ((un − 2un−1 + un−2 ) · ∇)un + ((un−1 − un−2 ) · ∇)(un − un−1 ), 17
18
R.An
and W n − W n−1 = SBn−1 × curl bn − Sbn−1 × curl enb −1 + S(bn−1 − bn−2 ) × curl ebn−1
+ Senb −1 × curl (bn − bn−1 ) + Sbn−1 × curl bn−1
+ S(bn − bn−1 ) × curl (bn − bn−1 ) + S(bn − 2bn−1 + bn ) × curl bn−1 , the terms 2Δt(X n − X n−1 , un ) and −2Δt(W n − W n−1 , un ) in (5.11) can be bounded as follows: Δt n−1 ||2 2 ) 2Δt(X n − X n−1 , un ||2L2 + CΔt||un−1 ||2L2 (|| Aun−1 ||2L2 + ||U un ) ≤ ||∇ H 4Re n −1 2 n n −1 2 n −1 2 + CΔt||eu || L2 || A(u − u )|| L2 + CΔt||u || L2 || Aun−1 ||2L2 eun−1 ||2L2 + CΔt||un − un−1 || L2 ||∇(un − un−1 )|| L2 ||∇
+ CΔt||un − 2un−1 + un−2 ||2L2 || Aun ||2L2
+ CΔt||un−1 − un−2 || L2 ||∇(un−1 − un−2 )|| L2 ||∇(un − un−1 )||2L2 Δt n−1 ||2 2 ) un ||2L2 + CΔt||un−1 ||2L2 (|| Aun−1 ||2L2 + ||U ≤ ||∇ H 4Re + C (Δt)4
tn
t n −1
eun−1 ||2L2 + C (Δt)5 , || Aut (t)||2L2 dt + C (Δt)3 ||∇
and
−2Δt(W n − W n−1 , un ) ≤
Δt ||∇ un ) un ||2L2 − 2SΔt(Bn−1 × curl bn , 4Re + CΔt||bn−1 ||2L2 (||bn−1 ||2H2 + ||Bn−1 ||2H2 )
+ CΔt||bn−1 − bn−2 || L2 ||curl (bn−1 − bn−2 )|| L2 ||curl enb −1 ||2L2 + CΔt||enb −1 ||2L2 ||bn − bn−1 ||2H2 + CΔt||bn−1 ||2L2 ||bn−1 ||2H2 + CΔt||bn − bn−1 || L2 ||curl (bn−1 − bn−2 )||3L2
+ CΔt||bn − 2bn−1 + bn−2 ||2L2 ||bn−1 ||2H2 Δt ≤ ||∇ un ) un ||2L2 − 2SΔt(Bn−1 × curl bn , 4Re + CΔt||bn−1 ||2L2 (||bn−1 ||2H2 + ||Bn−1 ||2H2 ) + C (Δt)5 + C (Δt)3 ||curl ebn−1 ||2L2 + C (Δt)4
tn
t n −1
||bt (t)||2H2 dt,
un , un ) = b(un−1 , un , un ) = 0 and where we use b(eun−1 ,
|| A(un − un−1 )||2L2 ≤ Δt ||bn − bn−1 )||2H2 ≤ Δt
tn
t n −1 tn t n −1
|| Aut (t)||2L2 dt, ||bt (t)||2H2 dt,
4 ||un − 2un−1 + un−2 ||2L2 ≤ (Δt)4 ||utt ||2L2 . 3 18
19
Projection method for MHD equations
Combining the above two estimates into (5.11) yields Δt ||∇ un ||2L2 2Re + (Δt)2 (||∇r np ||2L2 − ||∇r np−1 ||2L2 + ||∇(r np − r np−1 )||2L2 )
||un ||2L2 − ||un−1 ||2L2 + ||un − un−1 ||2L2 +
≤ C (Δt)3 ||∇r np−1 ||2L2 − 2SΔt(Bn−1 × curl bn , un ) + C (Δt)5 + C (Δt)
4
+ C (Δt)4
tn
t n −2 tn t n −1
(||uttt (t)||2L2 + ||∇ ptt (t)||2L2 )dt
(5.12)
(|| Aut (t)||2L2 + ||bt (t)||2H2 )dt
eun−1 ||2L2 + ||curl ebn−1 ||2L2 ) + C (Δt)3 (||∇
n −1 2 || H2 ) ||Bn−1 ||2H2 ).
+ CΔt||un−1 ||2L2 (|| Aun−1 ||2L2 + ||U
+ CΔt||bn−1 ||2L2 (||bn−1 ||2H2 +
Taking the sum of (5.9) and (5.12), we get Δt un ||2L2 ||∇ 2Re SΔt + S||bn ||2L2 − S||bn−1 ||2L2 + S||bn − bn−1 ||2L2 + ||curl bn ||2L2 Rm + (Δt)2 (||∇r np ||2L2 − ||∇r np−1 ||2L2 + ||∇(r np − r np−1 )||2L2 )
||un ||2L2 − ||un−1 ||2L2 + ||un − un−1 ||2L2 +
≤ C (Δt)3 ||∇r np−1 ||2L2 + C (Δt)5 + C (Δt)4 + C (Δt)4
tn
t n −2 tn t n −1
(||uttt (t)||2L2 + ||bttt (t)||2L2 + ||∇ ptt (t)||2L2 )dt (|| Aut (t)||2L2 + ||bt (t)||2H2 )dt
enu−1 ||2L2 + ||curl enb −1 ||2L2 ) + C (Δt)3 (||∇ + CΔt(||un−1 ||2L2 + ||bn−1 ||2L2 )(|| Aun−1 ||2L2 + ||U
+ CΔt||bn−1 ||2L2 (||bn−1 ||2H2 + ||Bn−1 ||2H2 ).
n −1 2 || H2 )
Summing up the above inequality for n from 2 to m ≤ N, using (4.9), (3.14) and the discrete Gronwall inequality, we obtain
||um ||2L2 ≤ C (Δt)4 + C ||u1 || L2 + C ||b1 ||2L2 + C (Δt)2 ||∇r1p ||2L2 = C (Δt)4 + C ||e1u || L2 + C ||e1b ||2L2 + C (Δt)2 ||∇r1p ||2L2 , where we use e0u = e0b = 0. To complete the proof of (5.5), we need to show
||e1u ||2L2 + ||e1b ||2L2 + (Δt)2 ||∇r1p ||2L2 ≤ C (Δt)4 . 19
(5.13)
20
R.An
Let n = 1 in (4.3) and (4.4). It follows from (4.6)-(4.8) and (5.3) that e1u −
Δt 1 e1u − Δt((u1 − u0 ) · ∇)u1 Δ e = Δt∇( P0 − p1 ) − Δt(u0 · ∇) Re u − SΔt(b0 × curl e1b + (b1 − b0 ) × curl b1 ) − ΔtR1u , (5.14)
and e1b +
Δt e1u × b0 ) − ΔtR1b curl (curl e1b ) = Δtcurl ( Rm + Δt((b1 − b0 ) · ∇)u1 − Δt(u1 · ∇)(b1 − b0 ). (5.15)
Testing (5.14) and (5.15) by e1u and Se1b , taking the sum of them, we find e1u ||2L2 + ||e1b ||2L2 ≤ C (Δt)2 (|| R1u ||2L2 + || R1b ||2L2 + ||∇( P0 − p0 )||2L2 + ||∇( p1 − p0 )||2L2 ) ||
+ C (Δt)2 (||∇(u1 − u0 )||2L2 + ||curl (b1 − b0 )||2L2 ) ≤ C (Δt)4 . e1u , then Since e1u = P H
e1u || L2 ≤ C (Δt)2 . ||e1u || L2 ≤ ||
On the other hand, in terms of (3.5), we have Δt||∇r1p || L2 ≤
1 e1u || L2 ) + Δt||∇( p1 − p0 )|| L2 ≤ C (Δt)2 . (||e1u || L2 + || 2
Therefore, we derive (5) and the complete the proof of (5.5).
From (5.5) and the proof of Theorem 4.2, it is easy to show Δt
N
∑ || enp ||2L
2
n =1
≤ C (Δt)2 ,
which implies that (3.15) in Theorem 3.3.
Acknowledgments The authors would like to thank the anonymous reviewers for their careful reviews and valuable comments to improve the quality of the original manuscript. This work is supported by Zhejiang Provincial Natural Science Foundation with Grant No. LY16A010017 and LY14A010020.
References [1] R.Adams, Sobolev Spaces, Academic Press, New York, 1975.
20
21
Projection method for MHD equations
[2] A.Chorin, Numerical solution of the Navier-Stokes equations, Math. Comp., 22(1968), 745-762. [3] J.Douglas Jr., J.Wang, An absolutely stabilized finite element method for the Stokes problem, Math. Comp., 52(1989), 495-508. [4] J.Gerbeau, A stabilized finite elemenet method for the incompressible magnetohydrodynamic equations, Numer. Math., 87(2000), 83-111. [5] J.Gerbeau, C.Le Bris, T.Leli`evre, Mathematical Methods for the Magnetohydrodynamics of Liquid Metals, Oxford: Oxford University Press, 2006. [6] J.Guermond, P.Minev, J.Shen, An overview of projection methods for incompressible flows, Comput. Methods Appl. Mech. Engrg., 195(2006), 6011-6045. [7] M.Gunzburger, A.Meir, J.Peterson, On the existence, uniqueness, and finite element approximation of solutions of the equations of stationary, incompresible magnetohydrodynamics, Math. Comp., 56(1991), 523-563. [8] Y.He, Unconditional convergence of the Euler semi-implicit scheme for the threedimensional incompressible MHD equations, IMA J. Numer. Anal., 35(2015),767-801. [9] J.Heywood, R.Rannacher, Finite-element approximation of the nonstationary NavierStokes problem Part IV: error analysis for second-order time discretization, SIAM J. Numer. Anal., 27(1990), 353-384. [10] W.Hughes, F.Young, The electromagnetics of fluids, Wiley, New York, 1966. [11] R.Moreau, Magneto-hydrodynamics, Kluwer Academic Publishers, 1990. [12] A.Prohl, Convergent finite element discretizations of the nonstationary incompressible magnetohydrodynamic system, ESAIM:M2AN, 42(2008), 1065-1087. ¨ [13] D.Schotzau, Mixed finite element methods for stationary incompressible magnetohydrodynamics, Numer. Math., 96(2004), 771-800. [14] M.Sermange, R.Temam, Some mathematical questions related to the MHD equations, Comm. Pure Appl. Math., 36(1983), 635-664. [15] J.Shen, On error estimates of projection methods for Navier-Stokes equations: first-order schemes, SIAM J. Numer. Anal., 29(1992):57-77. [16] J.Shen, Remarks on the pressure error estimates for the projection methods, Numer.Math., 67(1994), 513-520. [17] R.Temam, Sur l’approximation de la solution des equations de Navier-Stokes par la m´ethode des pas fractionnaires II, Arch. Rational Mech. Anal., 33(1969), 377-385. [18] R.Temam, Navier-Stokes Equations, North-Holland, Amsterdam, 1977. [19] 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. [20] Y.Zhang, Y.Hou, L.Shan, Numerical analysis of the Crank-Nicolson extrapolation time discrete scheme for magnetohydrodynamics flows, Numer. Methods Part. Diff. Equa., 31(2015), 2169-2208.
21