- Email: [email protected]
Two-level defect-correction Oseen iterative stabilized finite element methods for the stationary Navier–Stokes equations
Applied Mathematical Modelling 37 (2013) 728–741
Contents lists available at SciVerse ScienceDirect
Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm
Two-level defect-correction Oseen iterative stabilized finite element methods for the stationary Navier–Stokes equations q Pengzhan Huang a, Xinlong Feng a, Yinnian He b,a,⇑ a b
College of Mathematics and System Sciences, Xinjiang University, Urumqi 830046, PR China Center for Computational Geosciences, School of Mathematics and Statistics, Xi’an Jiaotong University, Xi’an 710049, PR China
a r t i c l e
i n f o
Article history: Received 30 May 2011 Received in revised form 18 February 2012 Accepted 29 February 2012 Available online 14 March 2012 Keywords: Defect-correction Navier–Stokes equations Local Gauss integration Two-level strategy Oseen iterative Error estimate
a b s t r a c t Two-level defect-correction Oseen iterative stabilized finite element methods for the stationary Navier–Stokes equations based on local Gauss integration are considered in this paper. The methods combine the defect-correction method and the two-level strategy with the locally stabilized method. Moreover, the stability and convergence of the presented methods are deduced. Finally, numerical tests confirm the theoretical results of the presented methods. Ó 2012 Elsevier Inc. All rights reserved.
1. Introduction In this report we consider the following stationary Navier–Stokes problem
mDu þ ðu rÞu þ rp ¼ f in X; div u ¼ 0 in X;
ð1Þ
u ¼ 0 on @ X; where X is a bounded, convex and open subset of R2 with a Lipschitz continuous boundary @ X, u ¼ ðu1 ðxÞ; u2 ðxÞÞ represents the velocity vector, p ¼ pðxÞ the pressure, f ¼ ðf1 ðxÞ; f2 ðxÞÞ the prescribed body force, and m > 0 the viscosity, which is inversely proportional to the Reynolds number. The incompressible Navier–Stokes problem is one of the main system studied in pipe flow, flow around airfoils, weather, blood flow and convective heat transfer inside industrial furnaces. Lots of works are devoted to this system, and the finite element methods and finite volume methods have been two classes of the most successful methods. However, solving the Navier–Stokes equations numerically is still a challenge because of theirs large computational scale, especially for large Reynolds number. Besides, it is well-known that the discretization of the Navier–Stokes equations by finite element methods and finite volume methods may generally bring out two shortcomings: the violation of the discrete inf-sup condition and spurious oscillations due to the domination of nonlinear convection term. q This work is in parts supported by the NSF of China (Nos. 61163027, 10901131, 10971166, 10961024), and the Natural Science Foundation of Xinjiang Province (No. 2010211B04). ⇑ Corresponding author at: Center for Computational Geosciences, School of Mathematics and Statistics, Xi’an Jiaotong University, Xi’an 710049, PR China. E-mail addresses: [email protected] (P. Huang), [email protected] (X. Feng), [email protected] (Y. He).
0307-904X/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.apm.2012.02.051
P. Huang et al. / Applied Mathematical Modelling 37 (2013) 728–741
729
To save the CPU time of a numerical method, a common choice is two-level method. The basic idea of two-level discretization strategy is to capture the ‘‘large eddies’’, ‘‘low modes’’, or ‘‘global solution envelope’’ by computing an initial approximation on a very coarse mesh. The fine structures are captured by solving one linear system on a fine mesh. Some details of two-level can be found in the works of Xu [1,2], Layton et al. [3–5], Zhang and He [6], Ervin et al. [7], He et al. [8,9], Li [10] and Shang [11]. The defect-correction method is an iterative improvement technique which is well established for solving nonlinear steady-state problems (see [12] for a survey of the technique). The basic idea of the defect-correction method is to solve stabilized artificial viscosity nonlinear equations in the defect step and correct the residual by a linearized problem in the correction step for a few steps (e.g., see [13,14]). This method has been successively studied in many references. Layton [15] initially investigated the defect-correction method for the incompressible Navier–Stokes equations with high Reynolds number. In [16], Axelsson and Layton applied the defect-correction method for convection–diffusion problems. Subsequently, Layton et al. also provided some further studies based on the defect-correction method for Navier–Stokes equations (see [14] and the references therein). In the recent, Kaya et al. [13] considered the synthesis of a subgrid stabilization method with defect-correction method for the stationary Navier–Stokes equations. Liu and Hou [17] studied a two-level defect-correction method for the steady-state Navier–Stokes equations based on the spectral method. Besides, Wang [18] presented a new defect-correction method for the Navier–Stokes equations. Employing finite element methods in solving the Navier–Stokes equations, the inf-sup condition has played an important role because it ensures a stability and accuracy of the underlying numerical schemes. A pair of finite element spaces that are used to approximate the velocity and the pressure unknowns are said to be stable if they satisfy the inf-sup condition. Intuitively speaking, the inf-sup condition is something that enforces a certain correlation between two finite element spaces so that they both have the required properties when employed for the Navier–Stokes equations. However, due to computational convenience and efficiency in practice, some mixed finite element pairs which do not satisfy the inf-sup condition are also popular. Thus, much attention has been paid to the study of the stabilized method for the Navier–Stokes problem. Recent studies have focused on stabilization of the lowest equal-order finite element pair P1 —P 1 (linear functions) or Q 1 —Q 1 (bilinear functions) using the projection of the pressure onto the piecewise constant space [19,20]. This stabilization technique is free of stabilization parameters and does not require any calculation of high-order derivatives or edge-based data structures. Therefore, this method is gaining more and more popularity in computational fluid dynamics. The method we study in this paper is to combine the defect-correction method with the two-level discretization for solving the two-dimensional steady Navier–Stokes problem based on the locally stabilized method with Gaussian quadrature rule, which includes three algorithms: m defect steps by Oseen iteration and once correction by Oseen iteration (mOseen-1-Oseen); m defect steps by Oseen iteration and once correction by Simple iteration (m-Oseen-1-Simple); m defect steps by Oseen iteration and once correction by Newton iteration (m-Oseen-1-Newton). The remainder of this paper is organized as follows. In Section 2, we introduce the notations, an abstract functional setting of the Navier–Stokes problem and some well-known results used throughout this paper. A stabilized finite element strategy is recalled in Section 3. Two-level defect-correction Oseen iterative algorithms are given in Section 4. Then in Section 5, numerical experiments are shown to verify the theoretical results completely. Finally, we end with a short conclusion in Section 6. 2. The functional setting of the Navier–Stokes equations As in [8,21], for the mathematical setting of problem (1), we introduce the following Hilbert spaces:
X ¼ H10 ðXÞ2 ;
Y ¼ L2 ðXÞ2 ;
M ¼ L20 ðXÞ ¼
q 2 L2 ðXÞ :
Z
q dx ¼ 0 :
X
The spaces L2 ðXÞm ; m ¼ 1; 2, are equipped with the L2 -scalar product ð; Þ and L2 -norm k k0 . The space X is endowed with the usual scalar product ðru; rv Þ and the norm kruk0 . Standard definitions are used for the Sobolev spaces W m;p ðXÞ, with the norm k km;p , m; p P 0. We will write Hm ðXÞ for W m;2 ðXÞ and k km for k km;2 . Next, let the closed subset V of X be given by
V ¼ fv 2 X : div v ¼ 0g; and denote by H the closed subset of Y, i.e.,
H ¼ fv 2 Y : div v ¼ 0;
v nj@X ¼ 0g:
We denote the Stokes operator by A ¼ PD, where P is the L2 -orthogonal projection of Y onto H and set DðAÞ ¼ H2 ðXÞ2 \ V. Moreover, we need a further assumption on X provided in [22,21]. Next, we make a regularity assumption on the Stokes problem. (A). Assume that X is smooth so that the unique solution ðv ; qÞ 2 X M of the steady Stokes problem
Dv þ rq ¼ g in X; div v ¼ 0 in X;
v ¼ 0 on @ X;
730
P. Huang et al. / Applied Mathematical Modelling 37 (2013) 728–741
for any prescribed g 2 Y exists and satisfies
kv k2 þ kqk1 6 c0 kgk0 ;
ð2Þ
where c0 > 0 is a constant depending on X. Subsequently, c (with or without a subscript) will denote a positive constant which may stand for different values at its different occurrences. We remark that the validity of assumption (A) is known [23] if @ X is of C 2 or if X is a two-dimensional convex polygon. From the assumption (A), it is well known [24] that
kv k2 6 c1 kAv k0 ;
krv k0 6 c2 kAv k0 ; 8v 2 DðAÞ;
ð3Þ
kv k0 6 c2 krv k0 ;
8v 2 X:
ð4Þ
We define the continuous bilinear forms að; Þ, dð; Þ and a generalized bilinear form Bðð; Þ; ð; ÞÞ on X X; X M and ðX MÞ ðX MÞ, respectively, by
aðu; v Þ ¼ mðru; rv Þ; dðv ; qÞ ¼ ðq; divv Þ;
8u; v 2 X; 8v 2 X; 8q 2 M;
and
Bððu; pÞ; ðv ; qÞÞ ¼ aðu; v Þ dðv ; pÞ þ dðu; qÞ;
8ðu; pÞ; ðv ; qÞ 2 X M;
and a trilinear form bð; ; Þ on X X X by
1 1 1 bðu; v ; wÞ ¼ ððu rÞv ; wÞ þ ððdivuÞv ; wÞ ¼ ððu rÞv ; wÞ ððu rÞw; v Þ; 2 2 2
8u; v ; w 2 X:
With the above notations, the variational formulation of problem (1) reads as follows: Find ðu; pÞ 2 X M such that for all ðv ; qÞ 2 X M,
Bððu; pÞ; ðv ; qÞÞ þ bðu; u; v Þ ¼ ðf ; v Þ:
ð5Þ
It is easy to verify that bð; ; Þ and Bðð; Þ; ð; ÞÞ satisfy the following important properties [8,21]:
bðu; v ; wÞ ¼ bðu; w; v Þ; jbðu; v ; wÞj þ jbðv ; u; wÞj þ jbðw; u; v Þj 6 c3 kruk0 kAv k0 kwk0 ;
ð6Þ
jbðu; v ; wÞj 6 Nkruk0 krv k0 krwk0 ;
ð7Þ
and
jBððu; pÞ; ðv ; qÞÞj 6 cðmkruk0 þ kpk0 Þðkrv k0 þ kqk0 Þ; sup ðv ;qÞ2XM
jBððu; pÞ; ðv ; qÞÞj P b0 ðmkruk0 þ kpk0 Þ; krv k0 þ kqk0
where N ¼ sup
u;v ;w2X
jbðu;v ;wÞj , kruk0 krv k0 krwk0
8ðu; pÞ; ðv ; qÞ 2 X M;
8ðu; pÞ 2 X M;
and b0 is a positive constant depends only on X and f.
The following existence and uniqueness of solution of (5) are classical results [21]. Theorem 2.1. Given f 2 X 0 , there exists at least a solution pair ðu; pÞ 2 X M which satisfies (5) and
kruk0 6 m1 kf k1 ;
kf k1 ¼ sup
ðf ; v Þ
v 2X krv k0
:
ð8Þ
And if m and f satisfy the following uniqueness condition:
c ¼ m2 Nkf k1 < 1;
ð9Þ
then the solution pair ðu; pÞ of problem (5) is unique. Furthermore, we have the following regularity result [21,9]. Theorem 2.2. Assume that (A) holds and f 2 Y. Then the solution pair ðu; pÞ 2 DðAÞ ðH1 ðXÞ \ MÞ of problem (5) satisfies the following regularity
mkAuk0 þ kpk1 6 ckf k0 :
ð10Þ
3. A stabilized finite element based on local Gauss integration From now on, H is a real positive parameter tending to 0. The finite element subspace X H M H of X M is characterized by K H , a partitioning of X into triangles K with the mesh size H, assumed to be uniformly regular in the usual sense.
P. Huang et al. / Applied Mathematical Modelling 37 (2013) 728–741
731
Moreover, the fine mesh partition K h can be thought of as generated from K H by a mesh refinement process. Also, we introduce finite element space X h M h based on K h . Furthermore, let K l be a finite element partition with mesh size l. Here, l ¼ h or H and H h. Then we define
n Þ2 \ X : uj 2 P1 ðKÞ2 ; X l ¼ u 2 C 0 ðX K
o
8K 2 K l ;
n Þ \ M : qj 2 P1 ðKÞ; M l ¼ q 2 C 0 ðX K
o
8K 2 K l ;
where P1 ðKÞ represents the space of linear functions on K. Let the discrete analogue of the space V be given by
V l ¼ fv l 2 X l : dðv l ; ql Þ ¼ 0; 8ql 2 Ml g: Next, let P l : Y ! V l denote the L2 -orthogonal projection defined by
ðP l u; v l Þ ¼ ðu; v l Þ;
8 u 2 Y;
v l 2 V l:
Moreover, a discrete analogue Al ¼ P l Dl of the Stokes operator A is defined through the condition that ðDl ul ; v l Þ ¼ ðrul ; rv l Þ for all ul ; v l 2 X l . Note that the lowest equal-order pair X l M l does not satisfy the discrete inf–sup condition
sup
dðv l ; ql Þ
v l 2X l krv l k0
P b1 kql k0 ;
8ql 2 Ml ;
where the constant b1 > 0 is independent of l. In order to fulfill this condition, a stabilized generalized bilinear term is used [10,25]:
Bl ððul ; pl Þ; ðv ; qÞÞ ¼ Bððul ; pl Þ; ðv ; qÞÞ þ Gðpl ; qÞ; where Gðpl ; qÞ can be defined by
Gðpl ; qÞ ¼ a
P
K2K l
Z
K;2
pl q dx
Z K;1
pl q dx ;
8pl ; q 2 Ml ;
R where a > 0 is a stabilization parameter and K;i gðxÞdx indicates a local Gauss integral over K that is exact for polynomials of degree i; i ¼ 1; 2. In particular, the trial function pl 2 M l must be projected to piecewise constant space W l defined below when i ¼ 1 for any q 2 Ml . Let Pl be a L2 -projection operator, which is defined by
ðp; ql Þ ¼ ðPl p; ql Þ;
8p 2 L2 ðXÞ; ql 2 W l :
ð11Þ
Here W l L2 ðXÞ denotes the piecewise constant space associated with the triangulation K l . The following properties of the projection operator Pl can be proved [20,25]:
kPl pk0 6 c5 kpk0 ; 8p 2 L2 ðXÞ; kp Pl pk0 6 c6 lkpk1 ;
ð12Þ
8p 2 H1 ðXÞ:
ð13Þ
As a result of (11), the bilinear form Gð; Þ can be expressed as
Gðpl ; qÞ ¼ aðpl Pl pl ; q Pl qÞ;
pl ; q 2 M l :
ð14Þ
Now, we consider the following stabilized problem: find ðul ; pl Þ 2 X l M l such that
Bl ððul ; pl Þ; ðv ; qÞÞ þ bðul ; ul ; v Þ þ
rl aðul ; v Þ ¼ ðf ; v Þ 8ðv ; qÞ 2 X l M l ; m
ð15Þ
where r > 0 is an artificial viscosity parameter which is used as a stabilizing factor. The following theorem establishes the weak coercivity of the bilinear form Bl ðð; Þ; ð; ÞÞ for the finite element pair X l Ml . Theorem 3.1 ([10,25]). For all ðul ; pl Þ; ðv ; qÞ 2 X l M l , the bilinear form Bl ðð; Þ; ð; ÞÞ satisfies the continuity property
jBl ððul ; pl Þ; ðv ; qÞÞj 6 c7 ðmkrul k0 þ kpl k0 Þðkrv k0 þ kqk0 Þ;
ð16Þ
and the weak coercivity property
sup ðv ;qÞ2X l M l
jBl ððul ; pl Þ; ðv ; qÞÞj P b2 ðmkrul k0 þ kpl k0 Þ; krv k0 þ kqk0
where b2 > 0 is independent of l.
ð17Þ
732
P. Huang et al. / Applied Mathematical Modelling 37 (2013) 728–741
As in [26,25,18], we have the following results. Theorem 3.2. Let the exact solution ðu; pÞ be in DðAÞ ðH1 ðXÞ \ MÞ. Then, ðul ; pl Þ satisfies the following stability and error estimate:
ðrl þ mÞkrul k0 6 kf k1 ;
ðrl þ mÞkAl ul k0 6 ckf k0 ;
ð18Þ
mku ul k0 þ lðmkrðu ul Þk0 þ kp pl k0 Þ 6 c8 l2 kf k0 :
ð19Þ
4. Two-level defect-correction Oseen iterative algorithms The two-level defect-correction Oseen iterative algorithms based on local Gauss integration we study are as follows. Algorithm 4.1 (m-Oseen-1-Oseen) m Step I. Solve a defect Navier–Stokes problem on coarse mesh, i.e., find ðum H ; pH Þ 2 X H M H such that for all ðv ; qÞ 2 X H M H ,
BH ððunH ; pnH Þ; ðv ; qÞÞ þ
rH n n aðuH ; v Þ þ bðun1 H ; uH ; v Þ ¼ ðf ; v Þ; m
ð20Þ
where n ¼ 1; 2; . . . ; m. Here, ðu0H ; p0H Þ is determined by
BH ððu0H ; p0H Þ; ðv ; qÞÞ þ
rH 0 aðuH ; v Þ ¼ ðf ; v Þ; 8ðv ; qÞ 2 X H M H : m
ð21Þ
Step II. Solve a correction Navier–Stokes problem on fine mesh, i.e., find ðumh ; pmh Þ 2 X h M h such that for all ðv ; qÞ 2 X h M h ,
Bh ððumh ; pmh Þ; ðv ; qÞÞ þ
rh rh m aðumh ; v Þ þ bðum aðuH ; v Þ: H ; umh ; v Þ ¼ ðf ; v Þ þ m m
ð22Þ
Theorem 4.1. Under the assumptions of Theorem 2.1, ðum H Þ defined by scheme (20), satisfies 1 krum H k0 6 ðrH þ mÞ kf k1 ;
kAH um H k0
ð23Þ
1
6 cðrH þ mÞ kf k0 :
ð24Þ
Proof. It follows from (21) that (20) holds for m ¼ 0. Assuming that (20) holds for m ¼ J, we want to prove that it holds for m ¼ J þ 1. Taking ðv ; qÞ ¼ ðuHJþ1 ; pHJþ1 Þ 2 X H M H in (20) with n ¼ J þ 1 and using (14), we arrive at 1 kruJþ1 H k0 6 ðrH þ mÞ kf k1 :
Also, taking
v ¼ AH uHJþ1
ðrH þ m
ð25Þ
and q ¼ 0 in (20) with n ¼ J þ 1 and using (9), we get
ÞkAH uJþ1 H k0
1=2
1=2
1=2
Jþ1 6 kf k0 þ c2 NkruJH k0 kruJþ1 H k0 kAH uH k0
6
m 2
kAH uJþ1 H k0 þ kf k0 þ c
N2
m4
kf k21 kf k0 6
m 2
6
m 2
J 2 Jþ1 1 2 kAH uJþ1 H k0 þ kf k0 þ c m N kruH k0 kruH k0
kAH uJþ1 H k0 þ ckf k0 :
Combining this estimate with (25) shows that (23) and (24) are valid for m ¼ J þ 1. Thus, we have completed the proof of the theorem. h m Next, we will establish the bound of the error ðum H uH ; pH pH Þ for all m P 1.
m Lemma 4.1. Under the assumptions of Theorem 2.1, ðum H ; pH Þ defined by scheme (20) satisfies the following bounds:
ðrH þ mÞkrðuH um H Þk0 6 1 kpH pm H k0 6 3b2
for all m P 1.
N
N
!m kf k1 2
ðrH þ mÞ !m
ðrH þ mÞ2
kf k1
kf k1 ;
kf k1 ;
ð26Þ ð27Þ
733
P. Huang et al. / Applied Mathematical Modelling 37 (2013) 728–741
Proof. Setting ðen ; gn Þ ¼ ðuH unH ; pH pnH Þ, subtracting (20) from (15) with
BH ððe ; g Þ; ðv ; qÞÞ þ bðe n
n
n1
; uH ; v Þ þ
n bðun1 h ;e ;
l ¼ H, we have
v Þ þ m rHaðe ; v Þ ¼ 0; 8n P 1: 1
n
ð28Þ
Taking ðv ; qÞ ¼ ðen ; gn Þ in (28), we arrive at
ðrH þ mÞkren k0 6 Nkren1 k0 kruH k0 6 NðrH þ mÞ1 kren1 k0 kf k1 ; which yields
!m
N
m
kre k0 6
ðrH þ mÞ2
kre0 k0 ;
kf k1
ð29Þ
for all m P 1. Next, we deduce from (7), (18), (21), (15) and the uniqueness condition that
ðrH þ mÞkre0 k0 6 NkruH k20 6
N ðrH þ mÞ2
kf k21 6 kf k1 :
Combining this inequality with (29) yields (26). Furthermore, using (17), (18), (23), (28) with n ¼ m, and Theorem 2.1, we get
kgm k0 6 b1 Nkrum1 k0 krem k0 þ Nkrem1 k0 kruH k0 þ rHkrem k0 2 H 6 b1 2
N ðrH þ mÞ
1 m kf k1 ðrH þ mÞkrem k0 þ ðrH þ mÞkrem1 k0 þ b1 2 rHkre k0 6 3b2 2
N ðrH þ mÞ
!m kf k1 2
kf k1 :
Thus, (27) follows. Combining Lemma 4.1 and Theorem 3.2, it is valid for the defect step of Algorithm 4.1 that. m Theorem 4.2. Let ðu; pÞ 2 DðAÞ ðH1 ðXÞ \ MÞ and ðum H ; pH Þ 2 X H M H be the solutions of the problem (5) and (20), respectively. Then, under the assumption of Lemma 4.1, the following estimate holds:
mkrðu
um H Þk0
þ kp
pm H k0
6 cH þ c
!m
N ðrH þ mÞ2
kf k1
kf k1 ;
for all m P 1. Theorem 4.3. Under the assumptions of Theorem 2.1, umh defined by scheme (22), satisfies
1 rh þ kf k1 ; rh þ m ðrH þ mÞðrh þ mÞ 1 rh þ kf k0 : kAh umh k0 6 c rh þ m ðrH þ mÞðrh þ mÞ
krumh k0 6
ð30Þ
Proof. We derive from (3), (7), (22) and Theorem 4.1 that
krumh k0 6 ðrh þ mÞ1 ð1 þ rhðrH þ mÞ1 Þkf k1 6 1=2
rH þ rh þ m kf k1 ; ðrH þ mÞðrh þ mÞ
1=2
1=2
þ rhkAH um ðrh þ mÞkAh umh k0 6 kf k0 þ c2 Nkrum H k0 krumh k0 kAh umh k0 H k0 rh rH þ rh þ m m þ kAh umh k0 ; 6 ckf k0 þ1þ rH þ m rH þ m 2 which is (30) after some simple calculation. h Theorem 4.4. Let ðu; pÞ 2 DðAÞ ðH1 ðXÞ \ MÞ and ðumh ; pmh Þ 2 X h M h be the solutions of the problem (5) and (22), respectively. Then, under the assumptions of Theorem 2.1, the following estimate holds: 2
mkrðu umh Þk0 þ kp pmh k0 6 cðh þ H Þ þ c
N ðrH þ mÞ2
!m kf k1
kf k1 :
ð31Þ
734
P. Huang et al. / Applied Mathematical Modelling 37 (2013) 728–741
Proof. Setting ðem ; gm Þ ¼ ðuh umh ; ph pmh Þ, subtracting (22) from (15) with
l ¼ h, we have
m m 1 m m Bh ððem ; gm Þ; ðv ; qÞÞ þ bðuh um H ; uh ; v Þ þ bðuH ; e ; v Þ þ m rhðaðe ; v Þ þ aðuH ; v ÞÞ ¼ 0:
ð32Þ
Taking ðv ; qÞ ¼ ðe ; g Þ in (32) and using (6), we arrive at m
m
m ðrh þ mÞkrem k0 6 c3 kAh uh k0 kuh um H k0 þ rhkruH k0
6 which yields
rh kf k þ c kA u k ðku uk0 þ ku uH k0 þ kuH um H k0 Þ; rH þ m 1 3 h h 0 h 2
m
!m
N
2
ðrh þ mÞkre k0 6 cðrh þ h þ H Þ þ c
ðrH þ mÞ2
kf k1
kf k1 :
ð33Þ
Furthermore, using (17), (18), (23), (32) and Theorem 4.3, we get
m m m m kgm k0 6 b1 Nkrum 2 H k0 kre k0 þ c 3 kuh uH k0 kAh uh k0 þ rhðkre k0 þ kruH k0 Þ !m N 2 2 kf k1 kf k1 : 6 cðrh þ h þ H Þ þ c ðrH þ mÞ2
ð34Þ
Combining (33) and (34) with Theorem 3.2, we finish the proof of the theorem. h Remark 4.1. Since it is valid that
1
m
kf k1
rH þ rh þ m r2 Hh kf k1 ¼ kf k > 0: ðrH þ mÞðrh þ mÞ mðrH þ mÞðrh þ mÞ 1
This result implies that the stability of Algorithm 4.1 is better than that of the general Galerkin finite element method. Remark 4.2. Setting
r ¼ 0 in (20)–(22), we get the general Oseen two-level scheme. By Theorems 4.2 and 4.4, we obtain
mkrðu umH Þk0 þ kp pmH k0 6 cH þ c
N
m
m
kf k1 2
kf k1 ;
and
mkrðu umh Þk0 þ kp pmh k0 6 cðh þ H2 Þ þ c
N
kf k1 2
m
m kf k1 :
If we choose H ¼ h, then we obtain the Oseen defect-correction method. Remark 4.3. From Theorem 3.2, for the usual stabilized finite element solution ðuh ; ph Þ, which involves solving one large Navier–Stokes problem on a fine mesh with mesh size h, we have the following error estimate:
mkrðu uh Þk0 þ kp ph k0 6 ch: Furthermore, if we choose H such that h ¼ OðH2 Þ for Algorithm 4.1, then we get the convergence rate of same order as the usual stabilized finite element method from Theorem 4.4. However, our method is more efficient than the one-level scheme. Algorithm 4.2 (m-Oseen-1-Simple) m Step I. Solve a defect Navier–Stokes problem on coarse mesh, i.e., find ðum H ; pH Þ 2 X H M H by (20). Step II. Solve a correction Navier–Stokes problem on fine mesh, i.e., find ðumh ; pmh Þ 2 X h M h such that for all ðv ; qÞ 2 X h M h ,
Bh ððumh ; pmh Þ; ðv ; qÞÞ þ
rh rh m m aðumh ; v Þ þ bðum aðuH ; v Þ: H ; uH ; v Þ ¼ ðf ; v Þ þ m m
ð35Þ
Next, we will study the stability and error estimate of Algorithm 4.2.
Theorem 4.5. Under the assumptions of Theorem 2.1, ðumh ; pmh Þ defined by scheme (35) satisfies
krumh k0 6
1 1 kf k1 ; þ rh þ m rH þ m
ð36Þ
P. Huang et al. / Applied Mathematical Modelling 37 (2013) 728–741
735
and the following error estimate:
N
2
mkrðu umh Þk0 þ kp pmh k0 6 cðh þ H Þ þ c
ðrH þ mÞ2
!m kf k1
kf k1 :
ð37Þ
Proof. Taking ðv ; qÞ ¼ ðumh ; pmh Þ in (35), using (7) and Theorem 4.1, we get
krumh k0 6 ðrh þ mÞ1 ð1 þ rhðrH þ mÞ1 þ NðrH þ mÞ2 kf k1 Þkf k1 6
r H þ r h þ 2m kf k1 ; ðrH þ mÞðrh þ mÞ
which is (36). Next, setting ðem ; gm Þ ¼ ðuh umh ; ph pmh Þ, we derive from (15) and (35) that m Bh ððem ; gm Þ; ðv ; qÞÞ þ m1 rhðaðem ; v Þ þ aðum H ; v ÞÞ þ bðuh uH ; uh ; v Þ þ bðuH ; uh uH ; v Þ m m þ bðuH um H ; uh ; v Þ þ bðuH ; uH uH ; v Þ ¼ 0;
8ðv ; qÞ 2 X h Mh :
ð38Þ
Setting ðv ; qÞ ¼ ðem ; gm Þ in (38), using (6) and (7) and Theorems 3.1 and 4.1, we obtain m m m ðrh þ mÞkrem k0 6 c3 kuh uH k0 ðkAh uh k0 þ kAH um H k0 Þ þ rhkruH k0 þ Nðkruh k0 þ kruH k0 ÞkrðuH uH Þk0 !m N 2 6 cðrh þ h þ H2 Þ þ c kf k1 kf k1 : ðrH þ mÞ2
ð39Þ
Moreover, it follows from (17) that
m m m m kgm k0 6 b1 c3 kuh uH k0 ðkAh uh k0 þ kAH um 2 H k0 Þ þ rhðkre k0 þ kruH k0 Þ þ Nðkruh k0 þ kruH k0 ÞkrðuH uH Þk0 !m N 2 2 kf k1 kf k1 : 6 cðrh þ h þ H Þ þ c ðrH þ mÞ2
ð40Þ
Hence, by combining (39) and (40) with Theorem 3.2, we have completed the proof of (37). The proof ends. h Algorithm 4.3 (m-Oseen-1-Newton) m Step I. Solve a defect Navier–Stokes problem on coarse mesh, i.e., find ðum H ; pH Þ 2 X H M H by (20). Step II. Solve a correction Navier–Stokes problem on fine mesh, i.e., find ðumh ; pmh Þ 2 X h M h such that for all ðv ; qÞ 2 X h M h ,
Bh ððumh ; pmh Þ; ðv ; qÞÞ þ
rh rh m m m aðumh ; v Þ þ bðumh ; um aðuH ; v Þ þ bðum H ; v Þ þ bðuH ; umh ; v Þ ¼ ðf ; v Þ þ H ; uH ; v Þ: m m
ð41Þ
Next, we will study the stability and error estimate of Algorithm 4.3. Theorem 4.6. Under the assumptions of Theorem 2.1, ðumh ; pmh Þ defined by scheme (41) satisfies
krumh k0 6
! 1 1 þ kf k1 ; Hm rH þ rrhþhmm rh þ rrHþ m
ð42Þ
and the following error estimate:
N
2
mkrðu umh Þk0 þ kp pmh k0 6 cðh þ H Þ þ c
ðrH þ mÞ2
!2m kf k1
kf k1 :
Proof. Taking ðv ; qÞ ¼ ðumh ; pmh Þ in (41) and using (7), we get m m 2 ðrh þ mÞkrumh k0 6 rhkrum H k0 þ Nkrumh k0 kruH k0 þ kf k1 þ NkruH k0 :
Using uniqueness condition and Theorem 4.1, we have
!
r2 Hh þ rHm þ rhm rh Nkf k1 kf k1 : krumh k0 6 þ1þ rH þ m rH þ m ðrH þ mÞ2 By some simple calculation, we have
ð43Þ
736
P. Huang et al. / Applied Mathematical Modelling 37 (2013) 728–741
krumh k0 6
rH þ rh þ 2m kf k ; r2 Hh þ rhm þ rHm 1
which is (42). Next, setting ðem ; gm Þ ¼ ðuh umh ; ph pmh Þ, we derive from (15) and (41) that m m m m Bh ððem ; gm Þ; ðv ; qÞÞ þ m1 rhaðem ; v Þ þ m1 rhaðum H ; v Þ þ bðe ; uh ; v Þ þ bðumh ; e ; v Þ þ bðumh uH ; umh uH ; v Þ ¼ 0:
ð44Þ
Setting ðv ; qÞ ¼ ðe ; g Þ in (44), using (6), (7), (9) and Theorems 3.1 and 4.1, we obtain m
m
m ðrh þ mÞð1 cÞkrem k20 6 Bh ððem ; gm Þ; ðem ; gm ÞÞ þ m1 rhaðem ; em Þ þ bðem ; um H;e Þ m m m 6 jbðuh uH ; uh uH ; em Þj þ jbðuh uH ; uH um H ; e Þj þ jbðuH uH ; uh uH ; e Þj m m 1 m m þ jbðuH um H ; uH uH ; e Þj þ jm rhaðuH ; e Þj 2 m 6 Nkrem k0 krðuh uH Þk20 þ krðuH um H Þk0 þ 2krðuH uH Þk0 krðuh uH Þk0 m þ rhkrum H k0 kre k0 :
ð45Þ
Moreover, it follows from (19), (23), (26) and Young inequality that
ðrh þ mÞkre k0 6 cðH þ rhÞ þ c
!2m
N
2
m
ðrH þ mÞ2
kf k1
kf k1 :
ð46Þ
Finally, one finds from (17) and (45) that
kgm k0 6 b1 2
rhðkrem k0 þ krumH k0 Þ þ Nkrem k0 krumH k0 þ Nðkrðuh uH Þk20 þ krðuH umH Þk20
þ 2krðuH um H Þk0 krðuh uH Þk0 Þ 6 cðH2 þ rhÞ þ c
N ðrH þ mÞ
!2m kf k1 2
kf k1 :
Hence, by combining this estimate and (46) with Theorem 3.2, we have completed the proof of (43). The proof ends.
h
Remark 4.4. From (30), (36) and (42), we find that the stability of Algorithm 4.1 is the best among these algorithms. h
Table 1 Comparisons of the one-level with the two-level defect-correction Oseen iterative methods. Methods
H
h
CPU-time
krðuuapp Þk0 kruk0
kppapp k0 kpk0
uH1 -Rate
pL2 -Rate
One-level Algorithm 4.1 Algorithm 4.2 Algorithm 4.3
– 1/4 1/4 1/4
1/16 1/16 1/16 1/16
0.454 0.235 0.204 0.219
0.249855 0.249881 0.249887 0.249826
0.0359607 0.0359257 0.0358945 0.0360782
– – – –
– – – –
One-level Algorithm 4.1 Algorithm 4.2 Algorithm 4.3
– 1/5 1/5 1/5
1/25 1/25 1/25 1/25
1.094 0.563 0.469 0.501
0.142012 0.142025 0.142018 0.142009
0.0165533 0.0165378 0.0165141 0.0166518
1.2659 1.2660 1.2661 1.2657
1.7384 1.7384 1.7396 1.7325
One-level Algorithm 4.1 Algorithm 4.2 Algorithm 4.3
– 1/6 1/6 1/6
1/36 1/36 1/36 1/36
2.329 1.141 0.954 1.100
0.0917065 0.0917145 0.0917115 0.0917041
0.00881857 0.00881067 0.00879852 0.00890549
1.1993 1.1993 1.1993 1.1993
1.7270 1.7269 1.7267 1.7163
One-level Algorithm 4.1 Algorithm 4.2 Algorithm 4.3
– 1/7 1/7 1/7
1/49 1/49 1/49 1/49
4.375 2.032 1.688 1.766
0.0643529 0.0643586 0.0643553 0.0643531
0.00520091 0.00519653 0.00519072 0.00528934
1.1489 1.1489 1.1490 1.1488
1.7127 1.7125 1.7117 1.6898
One-level Algorithm 4.1 Algorithm 4.2 Algorithm 4.3
– 1/8 1/8 1/8
1/64 1/64 1/64 1/64
7.469 3.422 2.875 2.985
0.0478061 0.0478105 0.0478077 0.0478079
0.00330493 0.00330241 0.00329995 0.00340521
1.1130 1.1129 1.1130 1.1128
1.6978 1.6975 1.6961 1.6490
One-level Algorithm 4.1 Algorithm 4.2 Algorithm 4.3
– 1/9 1/9 1/9
1/81 1/81 1/81 1/81
12.469 5.609 4.704 4.875
0.0370025 0.0370059 0.0370042 0.0370054
0.00222314 0.00222168 0.00222093 0.00234419
1.0875 1.0875 1.0874 1.0873
1.6832 1.6827 1.6810 1.5850
737
P. Huang et al. / Applied Mathematical Modelling 37 (2013) 728–741
Remark 4.5. If we choose H such that h ¼ OðH2 Þ for Algorithms 4.2 and 4.3, then we get the convergence rate of same order as the usual stabilized finite element method from Theorems 4.5 and 4.6. However, the given methods are more efficient than the one-level schemes. Remark 4.6. To the author’s knowledge, Algorithm 4.2 is the most efficient in these three methods for the small Reynolds number. Because the trilinear form is the most simple. Meanwhile, the coefficient matrices of Algorithms 4.2 and 4.3 are symmetric, while the one of Algorithm 4.1 is not.
5. Numerical experiments In this section we present numerical experiments to check the numerical theory developed in the previous sections and illustrate the efficiency of the two-level defect-correction Oseen iterative methods based on local Gauss integration. In the given experiments, the pressure and velocity are approximated by the lowest equal-order finite element pairs defined with respect to the same uniform triangulation. Meanwhile, the algorithms are implemented using public domain finite element software [27]. Besides, we take r ¼ 0:01 as in [14] and a ¼ 1. We consider the exact solution problem. Let X be the unit square in R2 . The exact solution for the velocity u ¼ ðu1 ; u2 Þ and the pressure p is given as follows:
pðx; yÞ ¼ 10ð2x 1Þð2y 1Þ; u1 ðx; yÞ ¼ 10x2 ðx 1Þ2 yðy 1Þð2y 1Þ; u2 ðx; yÞ ¼ 10xðx 1Þð2x 1Þy2 ðy 1Þ2 ; and the right-hand side f ¼ ðf1 ðx; yÞ; f2 ðx; yÞÞ is determined by original problem (1). Our goal in this test is to validate the properties and merits of the given two-level methods as compared with the onelevel method. Hence, we compare the solution generated at a fixed value of h for the one-level method with the results obtained using two-level methods where the finest grid has the same mesh spacing h. In all two-level computations, h is given, so H is chosen such that h ¼ H2 . Meanwhile, we have used the equal-order finite element P 1 P 1 with a fixed tolerance as 1
0.8
0.8
0.6
0.6 Y
Y
1
0.4
0.4
0.2
0.2
0
0
0.2
0.4
X
0.6
0.8
0
1
0.8
0.8
0.6
0.6
0.2
0.4
0
0.2
0.4
X
0.6
0.8
1
0.6
0.8
1
Y
1
Y
1
0
0.4
0.4
0.2
0.2
0
0
0.2
0.4
Fig. 1. Velocity streamlines for
X
0.6
0.8
1
0
X
m ¼ 1=3200: (a) standard Galerkin method, (b) Algorithm 4.2, (c) Algorithm 4.3, and (d) Algorithm 4.1.
738
P. Huang et al. / Applied Mathematical Modelling 37 (2013) 728–741
0.8
0.8
0.6
0.6 Y
1
Y
1
0.4
0.4
0.2
0.2
0
0
0.2
0.4
X
0.6
0.8
Fig. 2. Velocity streamlines for
0
1
0
0.2
0.4
X
0.6
0.8
1
m ¼ 1=5000: (a) Algorithm 4.2 and (b) Algorithm 4.1.
1
0.03
0.8
0.025 0.02 Y
Y
0.6
0.015
0.4 0.01 0.2
0
0.005
0
0.2
0.4
0.6
0.8
1
0
0.01
0.02
X Fig. 3. Velocity streamlines for
m ¼ 1=10000 by Algorithm 4.1: (a) whole and (b) detail in the left bottom corner.
1 0.9
0.03
X
0.6 Ghia et al. Algorithm 4.1
Ghia et al. Algorithm 4.1
0.4
0.8 0.2 v−velocity
0.7
y
0.6 0.5 0.4 0.3
0 −0.2 −0.4
0.2 −0.6
0.1 0 −0.5
0
u−velocity
0.5
1
−0.8
0
Fig. 4. Horizontal (a) and vertical velocity (b) for
0.2
0.4
x
0.6
0.8
1
m ¼ 1=3200.
1.0E5. We consider the case of the viscidity m ¼ 1, and pick six values of h, i.e., 1/16, 1/25, 1/36, 1/49, 1/64 and 1/81. The CPU-time, the L2 -error of the pressure, H1 -error of the velocity for the one-level method and the two-level methods for different values of h are tabulated in Table 1. From Table 1, we can see that the one-level method and two-level methods work well and keep the convergence rates just like the theoretical analysis. As expected, the two-level methods spend less time than one-level method under nearly the
739
P. Huang et al. / Applied Mathematical Modelling 37 (2013) 728–741
1
0.6 Ghia et al. Erturk et al. Algorithm 4.1
0.9 0.8
0.2 v−velocity
0.7
y
0.6 0.5 0.4 0.3
0 −0.2 −0.4
0.2
−0.6
0.1 0
Ghia et al. Erturk et al. Algorithm 4.1
0.4
−0.5
0
u−velocity
0.5
−0.8
1
0
Fig. 5. Horizontal (a) and vertical velocity (b) for
(a)
0.8
Ghia et al. Erturk et al. Algorithm 4.1
0.6
0.8
1
m ¼ 1=5000.
Ghia et al. Erturk et al. Algorithm 4.1
0.2 v−velocity
0.6 y
x
0.4
0.7
0.5 0.4 0.3
0 −0.2 −0.4
0.2
−0.6
0.1 0 −0.5
0.4
(b) 0.6
1 0.9
0.2
0
u−velocity
0.5
1
−0.8
0
Fig. 6. Horizontal (a) and vertical velocity (b) for
0.2
0.4
x
0.6
0.8
1
m ¼ 1=10000.
same relative error. Meanwhile, we can observe that Algorithms 4.2 and 4.3 are more efficient than Algorithm 4.1, which validate the analysis in Remark 4.6. However, from Remark 4.4, we know that the stability of Algorithm 4.1 is better than that of Algorithms 4.2 and 4.3. So, to show stability of Algorithm 4.1, we test a popular benchmark problem, the lid driven cavity, which has been analyzed in [28,29]. This problem is chosen because some benchmark data is available for comparison. Let the computation be carried out in the region X ¼ fðx; yÞj0 < x; y < 1g. We impose the normal component of the velocity to be zero on @ X and the tangential component to be zero except along y ¼ 1 where it is set to one. The exact solution of this problem is unknown. So, we take Ghia et al.’s algorithm [28] computed on a very fine mesh ðh ¼ 1=129Þ as the ‘‘exact’’ solution for the purpose of comparison. Meanwhile, we also give the results of the standard Galerkin method to show the good stability and efficiency of the two-level defect-correction Oseen iterative methods for the smaller viscidity problem. The computed results are shown, based on the standard Galerkin method, Algorithms 4.2 and 4.3, and Algorithm 4.1 for a set of different smaller viscidity coefficients m = 1/3200, 1/5000, 1/10000 with different mesh sizes h ¼ 1=64; 1=81; 1=121, respectively. In Figs. 1–3, we present the velocity streamlines of the cavity flow. The standard Galerkin method and Algorithm 4.3 are divergent despite m ¼ 1=5000. Algorithms 4.1 and 4.2 get almost the same results at this case. When the viscidity coefficient decreases to 1/10000, Algorithm 4.2 is divergent too. As we know, the position of the main vortex moves towards the center of the cavity when the viscidity decreases, and additional second vortex may appear in the right bottom corner of the cavity and a third vortex appears at the left bottom corner. When viscidity is 1/10000, six resolved vortices are captured. Above results for Algorithm 4.1 are consistent with those in [29] very well, which suggest its efficiency for solving the Navier–Stokes equations at high Reynolds numbers. In Figs. 4–6, we draw the x component of velocity along the vertical centerline and y component of velocity along the horizontal centerline based on Algorithm 4.1. Compared with the results obtained by Ghia et al. and Erturk et al., we can find that they are in excellent agreement. Moreover, the pressure contours of the cavity flow at different viscidity coefficients
740
P. Huang et al. / Applied Mathematical Modelling 37 (2013) 728–741
(a) 1
Y
0.6
0.4
0.2
0
0.2
0.4
X
0.6
0.8
0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.157104 0.120856 0.1 0.0802086 0.0546398 0.034689 0.0108654 0 -0.0316355 -0.058813 -0.0709235 -0.0784258
0.8
0.6 Y
0.8
0
(b) 1
0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.156319 0.131462 0.1 0.0824395 0.0751833 0.0681661 0.0533265 0.00204348 -0.029392 -0.0501363 -0.0622554 -0.0736008 -0.0814365
0.4
0.2
0
1
0
0.2
(c) 1
X
0.6
0.8
1
0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0698888 0.05293 0.051765 0.0411888 0.0403703 0.0378173 0.0339114 0.0292205 0.0180187 0 -0.0334256 -0.0533253 -0.0656891 -0.0760043
0.8
Y
0.6
0.4
0.2
0
0.4
0
Fig. 7. Pressure level lines for
0.2
0.4
X
0.6
0.8
1
m ¼ 1=3200 (a), m ¼ 1=5000 (b) and m ¼ 1=10000 (c) based on Algorithm 4.1.
are given in Fig. 7. From this figure, we can see that there is a singular point at the up right corner. Therefore, Algorithm 4.1 can capture this classical model well. 6. Conclusions In this work we have presented the two-level defect-correction Oseen iterative stabilized finite element methods in solving the steady Navier–Stokes equations based on local Gauss integration. The main feature of our method is to combine the defect-correction method and the two-level discretization. These methods are m-Oseen-1-Oseen, m-Oseen-1-Simple, and mOseen-1-Newton. The theoretical analysis and numerical tests show m-Oseen-1-Oseen is a favorite algorithm among these three algorithms for the considered problem. This idea can be extended to the time-dependent problems and other fluids dynamical equations. Acknowledgements The authors thank the editor and reviewers for their valuable comments and suggestions which helped us to improve the results of this paper. References [1] [2] [3] [4] [5] [6]
J. Xu, A novel two-grid method for semilinear elliptic equations, SIAM J. Sci. Comput. 15 (1994) 231–237. J. Xu, Two-grid discretization techniques for linear and nonlinear PDEs, SIAM J. Numer. Anal. 33 (1996) 1759–1778. W. Layton, A two level discretization method for the Navier–Stokes equations, Comput. Math. Appl. 26 (1993) 33–38. W. Layton, W. Lenferink, Two-level Picard and modified Picard methods for the Navier–Stokes equations, Appl. Math. Comput. 69 (1995) 263–274. W. Layton, L. Tobiska, A two-level method with backtracking for the Navier–Stokes equations, SIAM J. Numer. Anal. 35 (1998) 2035–2054. Y. Zhang, Y.N. He, A two-level finite element method for the stationary Navier–Stokes equations based on a stabilized local projection, Numer. Methods Part. Diff. Eq. 27 (2011) 460–477. [7] V. Ervin, W. Layton, J. Maubach, A posteriori error estimators for a two-level finite element method for the Navier–Stokes equations, Numer. Methods Part. Diff. Eq. 12 (1996) 333–346. [8] Y.N. He, K.T. Li, Two-level stabilized finite element methods for the steady Navier–Stokes problem, Computing 74 (2005) 337–351. [9] Y.N. He, A.W. Wang, A simplified two-level method for the steady Navier–Stokes equations, Comput. Methods Appl. Mech. Eng. 197 (2008) 1568–1576.
P. Huang et al. / Applied Mathematical Modelling 37 (2013) 728–741
741
[10] J. Li, Investigations on two kinds of two-level stabilized finite element methods for the stationary Navier–Stokes equations, Appl. Math. Comput. 182 (2006) 1470–1481. [11] Y.Q. Shang, A parallel two-level linearization method for incompressible flow problems, Appl. Math. Lett. 24 (2011) 364–369. [12] H.J. Stetter, The defect correction principle and discretization methods, Numer. Math. 29 (1978) 425–443. [13] S. Kaya, W. Layton, B. Riviere, Subgrid stabilized defect correction methods for the Navier–Stokes equations, SIAM J. Numer. Anal. 44 (2006) 1639– 1654. [14] W. Layton, H. Lee, J. Peterson, A defect-correction method for the incompressible Navier–Stokes equations, Appl. Math. Comput. 129 (2002) 1–19. [15] W. Layton, Solution algorithm for incompressible viscous flows at high Reynolds number, Vestnik Mosk. Gos. Univ. Ser. 15 (1996) 25–35. [16] O. Axelsson, W. Layton, Defect correction methods for convection dominated, convection diffusion equations, RAIRO J. Numer. Anal. 24 (1990) 423– 455. [17] Q.F. Liu, Y.R. Hou, A two-level defect-correction method for Navier–Stokes equations, Bull. Aust. Math. Soc. 81 (2010) 442–454. [18] K. Wang, A new defect correction method for the Navier–Stokes equations at high Reynolds numbers, Appl. Math. Comput. 216 (2010) 3252–3264. [19] P. Bochev, C.R. Dohrmann, M.D. Gunzburger, Stabilization of low-order mixed finite elements for the Stokes equations, SIAM J. Numer. Anal. 44 (2006) 82–101. [20] J. Li, Y.N. He, Z.X. Chen, A new stabilized finite element method for the transient Navier–Stokes equations, Comput. Methods Appl. Mech. Eng. 197 (2007) 22–35. [21] Y.N. He, J. Li, Convergence of three iterative methods based on the finite element discretization for the stationary Navier–Stokes equations, Comput. Methods Appl. Mech. Eng. 198 (2009) 1351–1359. [22] J.G. Heywood, R. Rannacher, Finite-element approximations of the nonstationary Navier–Stokes problem. Part I: Regularity of solutions and secondorder spatial discretization, SIAM J. Numer. Anal. 19 (1982) 275–311. [23] R.B. Kellogg, J.E. Osborn, A regularity result for the Stokes problem in a convex polygon, J. Funct. Anal. 21 (1976) 397–431. [24] S. Larsson, The long-time behavior of finite-element approximations of solutions to semilinear parabolic problems, SIAM J. Numer. Anal. 26 (1989) 348–365. [25] J. Li, Y.N. He, A stabilized finite element method based on two local Gauss integrations for the Stokes equations, J. Comput. Appl. Math. 214 (2008) 58– 65. [26] Y.N. He, A.W. Wang, L.Q. Mei, Stabilized finite-element method for the stationary Navier–Stokes equations, J. Eng. Math. 51 (2005) 367–380. [27] FreeFem++, version 2.23.
Contents lists available at SciVerse ScienceDirect
Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm
Two-level defect-correction Oseen iterative stabilized finite element methods for the stationary Navier–Stokes equations q Pengzhan Huang a, Xinlong Feng a, Yinnian He b,a,⇑ a b
College of Mathematics and System Sciences, Xinjiang University, Urumqi 830046, PR China Center for Computational Geosciences, School of Mathematics and Statistics, Xi’an Jiaotong University, Xi’an 710049, PR China
a r t i c l e
i n f o
Article history: Received 30 May 2011 Received in revised form 18 February 2012 Accepted 29 February 2012 Available online 14 March 2012 Keywords: Defect-correction Navier–Stokes equations Local Gauss integration Two-level strategy Oseen iterative Error estimate
a b s t r a c t Two-level defect-correction Oseen iterative stabilized finite element methods for the stationary Navier–Stokes equations based on local Gauss integration are considered in this paper. The methods combine the defect-correction method and the two-level strategy with the locally stabilized method. Moreover, the stability and convergence of the presented methods are deduced. Finally, numerical tests confirm the theoretical results of the presented methods. Ó 2012 Elsevier Inc. All rights reserved.
1. Introduction In this report we consider the following stationary Navier–Stokes problem
mDu þ ðu rÞu þ rp ¼ f in X; div u ¼ 0 in X;
ð1Þ
u ¼ 0 on @ X; where X is a bounded, convex and open subset of R2 with a Lipschitz continuous boundary @ X, u ¼ ðu1 ðxÞ; u2 ðxÞÞ represents the velocity vector, p ¼ pðxÞ the pressure, f ¼ ðf1 ðxÞ; f2 ðxÞÞ the prescribed body force, and m > 0 the viscosity, which is inversely proportional to the Reynolds number. The incompressible Navier–Stokes problem is one of the main system studied in pipe flow, flow around airfoils, weather, blood flow and convective heat transfer inside industrial furnaces. Lots of works are devoted to this system, and the finite element methods and finite volume methods have been two classes of the most successful methods. However, solving the Navier–Stokes equations numerically is still a challenge because of theirs large computational scale, especially for large Reynolds number. Besides, it is well-known that the discretization of the Navier–Stokes equations by finite element methods and finite volume methods may generally bring out two shortcomings: the violation of the discrete inf-sup condition and spurious oscillations due to the domination of nonlinear convection term. q This work is in parts supported by the NSF of China (Nos. 61163027, 10901131, 10971166, 10961024), and the Natural Science Foundation of Xinjiang Province (No. 2010211B04). ⇑ Corresponding author at: Center for Computational Geosciences, School of Mathematics and Statistics, Xi’an Jiaotong University, Xi’an 710049, PR China. E-mail addresses: [email protected] (P. Huang), [email protected] (X. Feng), [email protected] (Y. He).
0307-904X/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.apm.2012.02.051
P. Huang et al. / Applied Mathematical Modelling 37 (2013) 728–741
729
To save the CPU time of a numerical method, a common choice is two-level method. The basic idea of two-level discretization strategy is to capture the ‘‘large eddies’’, ‘‘low modes’’, or ‘‘global solution envelope’’ by computing an initial approximation on a very coarse mesh. The fine structures are captured by solving one linear system on a fine mesh. Some details of two-level can be found in the works of Xu [1,2], Layton et al. [3–5], Zhang and He [6], Ervin et al. [7], He et al. [8,9], Li [10] and Shang [11]. The defect-correction method is an iterative improvement technique which is well established for solving nonlinear steady-state problems (see [12] for a survey of the technique). The basic idea of the defect-correction method is to solve stabilized artificial viscosity nonlinear equations in the defect step and correct the residual by a linearized problem in the correction step for a few steps (e.g., see [13,14]). This method has been successively studied in many references. Layton [15] initially investigated the defect-correction method for the incompressible Navier–Stokes equations with high Reynolds number. In [16], Axelsson and Layton applied the defect-correction method for convection–diffusion problems. Subsequently, Layton et al. also provided some further studies based on the defect-correction method for Navier–Stokes equations (see [14] and the references therein). In the recent, Kaya et al. [13] considered the synthesis of a subgrid stabilization method with defect-correction method for the stationary Navier–Stokes equations. Liu and Hou [17] studied a two-level defect-correction method for the steady-state Navier–Stokes equations based on the spectral method. Besides, Wang [18] presented a new defect-correction method for the Navier–Stokes equations. Employing finite element methods in solving the Navier–Stokes equations, the inf-sup condition has played an important role because it ensures a stability and accuracy of the underlying numerical schemes. A pair of finite element spaces that are used to approximate the velocity and the pressure unknowns are said to be stable if they satisfy the inf-sup condition. Intuitively speaking, the inf-sup condition is something that enforces a certain correlation between two finite element spaces so that they both have the required properties when employed for the Navier–Stokes equations. However, due to computational convenience and efficiency in practice, some mixed finite element pairs which do not satisfy the inf-sup condition are also popular. Thus, much attention has been paid to the study of the stabilized method for the Navier–Stokes problem. Recent studies have focused on stabilization of the lowest equal-order finite element pair P1 —P 1 (linear functions) or Q 1 —Q 1 (bilinear functions) using the projection of the pressure onto the piecewise constant space [19,20]. This stabilization technique is free of stabilization parameters and does not require any calculation of high-order derivatives or edge-based data structures. Therefore, this method is gaining more and more popularity in computational fluid dynamics. The method we study in this paper is to combine the defect-correction method with the two-level discretization for solving the two-dimensional steady Navier–Stokes problem based on the locally stabilized method with Gaussian quadrature rule, which includes three algorithms: m defect steps by Oseen iteration and once correction by Oseen iteration (mOseen-1-Oseen); m defect steps by Oseen iteration and once correction by Simple iteration (m-Oseen-1-Simple); m defect steps by Oseen iteration and once correction by Newton iteration (m-Oseen-1-Newton). The remainder of this paper is organized as follows. In Section 2, we introduce the notations, an abstract functional setting of the Navier–Stokes problem and some well-known results used throughout this paper. A stabilized finite element strategy is recalled in Section 3. Two-level defect-correction Oseen iterative algorithms are given in Section 4. Then in Section 5, numerical experiments are shown to verify the theoretical results completely. Finally, we end with a short conclusion in Section 6. 2. The functional setting of the Navier–Stokes equations As in [8,21], for the mathematical setting of problem (1), we introduce the following Hilbert spaces:
X ¼ H10 ðXÞ2 ;
Y ¼ L2 ðXÞ2 ;
M ¼ L20 ðXÞ ¼
q 2 L2 ðXÞ :
Z
q dx ¼ 0 :
X
The spaces L2 ðXÞm ; m ¼ 1; 2, are equipped with the L2 -scalar product ð; Þ and L2 -norm k k0 . The space X is endowed with the usual scalar product ðru; rv Þ and the norm kruk0 . Standard definitions are used for the Sobolev spaces W m;p ðXÞ, with the norm k km;p , m; p P 0. We will write Hm ðXÞ for W m;2 ðXÞ and k km for k km;2 . Next, let the closed subset V of X be given by
V ¼ fv 2 X : div v ¼ 0g; and denote by H the closed subset of Y, i.e.,
H ¼ fv 2 Y : div v ¼ 0;
v nj@X ¼ 0g:
We denote the Stokes operator by A ¼ PD, where P is the L2 -orthogonal projection of Y onto H and set DðAÞ ¼ H2 ðXÞ2 \ V. Moreover, we need a further assumption on X provided in [22,21]. Next, we make a regularity assumption on the Stokes problem. (A). Assume that X is smooth so that the unique solution ðv ; qÞ 2 X M of the steady Stokes problem
Dv þ rq ¼ g in X; div v ¼ 0 in X;
v ¼ 0 on @ X;
730
P. Huang et al. / Applied Mathematical Modelling 37 (2013) 728–741
for any prescribed g 2 Y exists and satisfies
kv k2 þ kqk1 6 c0 kgk0 ;
ð2Þ
where c0 > 0 is a constant depending on X. Subsequently, c (with or without a subscript) will denote a positive constant which may stand for different values at its different occurrences. We remark that the validity of assumption (A) is known [23] if @ X is of C 2 or if X is a two-dimensional convex polygon. From the assumption (A), it is well known [24] that
kv k2 6 c1 kAv k0 ;
krv k0 6 c2 kAv k0 ; 8v 2 DðAÞ;
ð3Þ
kv k0 6 c2 krv k0 ;
8v 2 X:
ð4Þ
We define the continuous bilinear forms að; Þ, dð; Þ and a generalized bilinear form Bðð; Þ; ð; ÞÞ on X X; X M and ðX MÞ ðX MÞ, respectively, by
aðu; v Þ ¼ mðru; rv Þ; dðv ; qÞ ¼ ðq; divv Þ;
8u; v 2 X; 8v 2 X; 8q 2 M;
and
Bððu; pÞ; ðv ; qÞÞ ¼ aðu; v Þ dðv ; pÞ þ dðu; qÞ;
8ðu; pÞ; ðv ; qÞ 2 X M;
and a trilinear form bð; ; Þ on X X X by
1 1 1 bðu; v ; wÞ ¼ ððu rÞv ; wÞ þ ððdivuÞv ; wÞ ¼ ððu rÞv ; wÞ ððu rÞw; v Þ; 2 2 2
8u; v ; w 2 X:
With the above notations, the variational formulation of problem (1) reads as follows: Find ðu; pÞ 2 X M such that for all ðv ; qÞ 2 X M,
Bððu; pÞ; ðv ; qÞÞ þ bðu; u; v Þ ¼ ðf ; v Þ:
ð5Þ
It is easy to verify that bð; ; Þ and Bðð; Þ; ð; ÞÞ satisfy the following important properties [8,21]:
bðu; v ; wÞ ¼ bðu; w; v Þ; jbðu; v ; wÞj þ jbðv ; u; wÞj þ jbðw; u; v Þj 6 c3 kruk0 kAv k0 kwk0 ;
ð6Þ
jbðu; v ; wÞj 6 Nkruk0 krv k0 krwk0 ;
ð7Þ
and
jBððu; pÞ; ðv ; qÞÞj 6 cðmkruk0 þ kpk0 Þðkrv k0 þ kqk0 Þ; sup ðv ;qÞ2XM
jBððu; pÞ; ðv ; qÞÞj P b0 ðmkruk0 þ kpk0 Þ; krv k0 þ kqk0
where N ¼ sup
u;v ;w2X
jbðu;v ;wÞj , kruk0 krv k0 krwk0
8ðu; pÞ; ðv ; qÞ 2 X M;
8ðu; pÞ 2 X M;
and b0 is a positive constant depends only on X and f.
The following existence and uniqueness of solution of (5) are classical results [21]. Theorem 2.1. Given f 2 X 0 , there exists at least a solution pair ðu; pÞ 2 X M which satisfies (5) and
kruk0 6 m1 kf k1 ;
kf k1 ¼ sup
ðf ; v Þ
v 2X krv k0
:
ð8Þ
And if m and f satisfy the following uniqueness condition:
c ¼ m2 Nkf k1 < 1;
ð9Þ
then the solution pair ðu; pÞ of problem (5) is unique. Furthermore, we have the following regularity result [21,9]. Theorem 2.2. Assume that (A) holds and f 2 Y. Then the solution pair ðu; pÞ 2 DðAÞ ðH1 ðXÞ \ MÞ of problem (5) satisfies the following regularity
mkAuk0 þ kpk1 6 ckf k0 :
ð10Þ
3. A stabilized finite element based on local Gauss integration From now on, H is a real positive parameter tending to 0. The finite element subspace X H M H of X M is characterized by K H , a partitioning of X into triangles K with the mesh size H, assumed to be uniformly regular in the usual sense.
P. Huang et al. / Applied Mathematical Modelling 37 (2013) 728–741
731
Moreover, the fine mesh partition K h can be thought of as generated from K H by a mesh refinement process. Also, we introduce finite element space X h M h based on K h . Furthermore, let K l be a finite element partition with mesh size l. Here, l ¼ h or H and H h. Then we define
n Þ2 \ X : uj 2 P1 ðKÞ2 ; X l ¼ u 2 C 0 ðX K
o
8K 2 K l ;
n Þ \ M : qj 2 P1 ðKÞ; M l ¼ q 2 C 0 ðX K
o
8K 2 K l ;
where P1 ðKÞ represents the space of linear functions on K. Let the discrete analogue of the space V be given by
V l ¼ fv l 2 X l : dðv l ; ql Þ ¼ 0; 8ql 2 Ml g: Next, let P l : Y ! V l denote the L2 -orthogonal projection defined by
ðP l u; v l Þ ¼ ðu; v l Þ;
8 u 2 Y;
v l 2 V l:
Moreover, a discrete analogue Al ¼ P l Dl of the Stokes operator A is defined through the condition that ðDl ul ; v l Þ ¼ ðrul ; rv l Þ for all ul ; v l 2 X l . Note that the lowest equal-order pair X l M l does not satisfy the discrete inf–sup condition
sup
dðv l ; ql Þ
v l 2X l krv l k0
P b1 kql k0 ;
8ql 2 Ml ;
where the constant b1 > 0 is independent of l. In order to fulfill this condition, a stabilized generalized bilinear term is used [10,25]:
Bl ððul ; pl Þ; ðv ; qÞÞ ¼ Bððul ; pl Þ; ðv ; qÞÞ þ Gðpl ; qÞ; where Gðpl ; qÞ can be defined by
Gðpl ; qÞ ¼ a
P
K2K l
Z
K;2
pl q dx
Z K;1
pl q dx ;
8pl ; q 2 Ml ;
R where a > 0 is a stabilization parameter and K;i gðxÞdx indicates a local Gauss integral over K that is exact for polynomials of degree i; i ¼ 1; 2. In particular, the trial function pl 2 M l must be projected to piecewise constant space W l defined below when i ¼ 1 for any q 2 Ml . Let Pl be a L2 -projection operator, which is defined by
ðp; ql Þ ¼ ðPl p; ql Þ;
8p 2 L2 ðXÞ; ql 2 W l :
ð11Þ
Here W l L2 ðXÞ denotes the piecewise constant space associated with the triangulation K l . The following properties of the projection operator Pl can be proved [20,25]:
kPl pk0 6 c5 kpk0 ; 8p 2 L2 ðXÞ; kp Pl pk0 6 c6 lkpk1 ;
ð12Þ
8p 2 H1 ðXÞ:
ð13Þ
As a result of (11), the bilinear form Gð; Þ can be expressed as
Gðpl ; qÞ ¼ aðpl Pl pl ; q Pl qÞ;
pl ; q 2 M l :
ð14Þ
Now, we consider the following stabilized problem: find ðul ; pl Þ 2 X l M l such that
Bl ððul ; pl Þ; ðv ; qÞÞ þ bðul ; ul ; v Þ þ
rl aðul ; v Þ ¼ ðf ; v Þ 8ðv ; qÞ 2 X l M l ; m
ð15Þ
where r > 0 is an artificial viscosity parameter which is used as a stabilizing factor. The following theorem establishes the weak coercivity of the bilinear form Bl ðð; Þ; ð; ÞÞ for the finite element pair X l Ml . Theorem 3.1 ([10,25]). For all ðul ; pl Þ; ðv ; qÞ 2 X l M l , the bilinear form Bl ðð; Þ; ð; ÞÞ satisfies the continuity property
jBl ððul ; pl Þ; ðv ; qÞÞj 6 c7 ðmkrul k0 þ kpl k0 Þðkrv k0 þ kqk0 Þ;
ð16Þ
and the weak coercivity property
sup ðv ;qÞ2X l M l
jBl ððul ; pl Þ; ðv ; qÞÞj P b2 ðmkrul k0 þ kpl k0 Þ; krv k0 þ kqk0
where b2 > 0 is independent of l.
ð17Þ
732
P. Huang et al. / Applied Mathematical Modelling 37 (2013) 728–741
As in [26,25,18], we have the following results. Theorem 3.2. Let the exact solution ðu; pÞ be in DðAÞ ðH1 ðXÞ \ MÞ. Then, ðul ; pl Þ satisfies the following stability and error estimate:
ðrl þ mÞkrul k0 6 kf k1 ;
ðrl þ mÞkAl ul k0 6 ckf k0 ;
ð18Þ
mku ul k0 þ lðmkrðu ul Þk0 þ kp pl k0 Þ 6 c8 l2 kf k0 :
ð19Þ
4. Two-level defect-correction Oseen iterative algorithms The two-level defect-correction Oseen iterative algorithms based on local Gauss integration we study are as follows. Algorithm 4.1 (m-Oseen-1-Oseen) m Step I. Solve a defect Navier–Stokes problem on coarse mesh, i.e., find ðum H ; pH Þ 2 X H M H such that for all ðv ; qÞ 2 X H M H ,
BH ððunH ; pnH Þ; ðv ; qÞÞ þ
rH n n aðuH ; v Þ þ bðun1 H ; uH ; v Þ ¼ ðf ; v Þ; m
ð20Þ
where n ¼ 1; 2; . . . ; m. Here, ðu0H ; p0H Þ is determined by
BH ððu0H ; p0H Þ; ðv ; qÞÞ þ
rH 0 aðuH ; v Þ ¼ ðf ; v Þ; 8ðv ; qÞ 2 X H M H : m
ð21Þ
Step II. Solve a correction Navier–Stokes problem on fine mesh, i.e., find ðumh ; pmh Þ 2 X h M h such that for all ðv ; qÞ 2 X h M h ,
Bh ððumh ; pmh Þ; ðv ; qÞÞ þ
rh rh m aðumh ; v Þ þ bðum aðuH ; v Þ: H ; umh ; v Þ ¼ ðf ; v Þ þ m m
ð22Þ
Theorem 4.1. Under the assumptions of Theorem 2.1, ðum H Þ defined by scheme (20), satisfies 1 krum H k0 6 ðrH þ mÞ kf k1 ;
kAH um H k0
ð23Þ
1
6 cðrH þ mÞ kf k0 :
ð24Þ
Proof. It follows from (21) that (20) holds for m ¼ 0. Assuming that (20) holds for m ¼ J, we want to prove that it holds for m ¼ J þ 1. Taking ðv ; qÞ ¼ ðuHJþ1 ; pHJþ1 Þ 2 X H M H in (20) with n ¼ J þ 1 and using (14), we arrive at 1 kruJþ1 H k0 6 ðrH þ mÞ kf k1 :
Also, taking
v ¼ AH uHJþ1
ðrH þ m
ð25Þ
and q ¼ 0 in (20) with n ¼ J þ 1 and using (9), we get
ÞkAH uJþ1 H k0
1=2
1=2
1=2
Jþ1 6 kf k0 þ c2 NkruJH k0 kruJþ1 H k0 kAH uH k0
6
m 2
kAH uJþ1 H k0 þ kf k0 þ c
N2
m4
kf k21 kf k0 6
m 2
6
m 2
J 2 Jþ1 1 2 kAH uJþ1 H k0 þ kf k0 þ c m N kruH k0 kruH k0
kAH uJþ1 H k0 þ ckf k0 :
Combining this estimate with (25) shows that (23) and (24) are valid for m ¼ J þ 1. Thus, we have completed the proof of the theorem. h m Next, we will establish the bound of the error ðum H uH ; pH pH Þ for all m P 1.
m Lemma 4.1. Under the assumptions of Theorem 2.1, ðum H ; pH Þ defined by scheme (20) satisfies the following bounds:
ðrH þ mÞkrðuH um H Þk0 6 1 kpH pm H k0 6 3b2
for all m P 1.
N
N
!m kf k1 2
ðrH þ mÞ !m
ðrH þ mÞ2
kf k1
kf k1 ;
kf k1 ;
ð26Þ ð27Þ
733
P. Huang et al. / Applied Mathematical Modelling 37 (2013) 728–741
Proof. Setting ðen ; gn Þ ¼ ðuH unH ; pH pnH Þ, subtracting (20) from (15) with
BH ððe ; g Þ; ðv ; qÞÞ þ bðe n
n
n1
; uH ; v Þ þ
n bðun1 h ;e ;
l ¼ H, we have
v Þ þ m rHaðe ; v Þ ¼ 0; 8n P 1: 1
n
ð28Þ
Taking ðv ; qÞ ¼ ðen ; gn Þ in (28), we arrive at
ðrH þ mÞkren k0 6 Nkren1 k0 kruH k0 6 NðrH þ mÞ1 kren1 k0 kf k1 ; which yields
!m
N
m
kre k0 6
ðrH þ mÞ2
kre0 k0 ;
kf k1
ð29Þ
for all m P 1. Next, we deduce from (7), (18), (21), (15) and the uniqueness condition that
ðrH þ mÞkre0 k0 6 NkruH k20 6
N ðrH þ mÞ2
kf k21 6 kf k1 :
Combining this inequality with (29) yields (26). Furthermore, using (17), (18), (23), (28) with n ¼ m, and Theorem 2.1, we get
kgm k0 6 b1 Nkrum1 k0 krem k0 þ Nkrem1 k0 kruH k0 þ rHkrem k0 2 H 6 b1 2
N ðrH þ mÞ
1 m kf k1 ðrH þ mÞkrem k0 þ ðrH þ mÞkrem1 k0 þ b1 2 rHkre k0 6 3b2 2
N ðrH þ mÞ
!m kf k1 2
kf k1 :
Thus, (27) follows. Combining Lemma 4.1 and Theorem 3.2, it is valid for the defect step of Algorithm 4.1 that. m Theorem 4.2. Let ðu; pÞ 2 DðAÞ ðH1 ðXÞ \ MÞ and ðum H ; pH Þ 2 X H M H be the solutions of the problem (5) and (20), respectively. Then, under the assumption of Lemma 4.1, the following estimate holds:
mkrðu
um H Þk0
þ kp
pm H k0
6 cH þ c
!m
N ðrH þ mÞ2
kf k1
kf k1 ;
for all m P 1. Theorem 4.3. Under the assumptions of Theorem 2.1, umh defined by scheme (22), satisfies
1 rh þ kf k1 ; rh þ m ðrH þ mÞðrh þ mÞ 1 rh þ kf k0 : kAh umh k0 6 c rh þ m ðrH þ mÞðrh þ mÞ
krumh k0 6
ð30Þ
Proof. We derive from (3), (7), (22) and Theorem 4.1 that
krumh k0 6 ðrh þ mÞ1 ð1 þ rhðrH þ mÞ1 Þkf k1 6 1=2
rH þ rh þ m kf k1 ; ðrH þ mÞðrh þ mÞ
1=2
1=2
þ rhkAH um ðrh þ mÞkAh umh k0 6 kf k0 þ c2 Nkrum H k0 krumh k0 kAh umh k0 H k0 rh rH þ rh þ m m þ kAh umh k0 ; 6 ckf k0 þ1þ rH þ m rH þ m 2 which is (30) after some simple calculation. h Theorem 4.4. Let ðu; pÞ 2 DðAÞ ðH1 ðXÞ \ MÞ and ðumh ; pmh Þ 2 X h M h be the solutions of the problem (5) and (22), respectively. Then, under the assumptions of Theorem 2.1, the following estimate holds: 2
mkrðu umh Þk0 þ kp pmh k0 6 cðh þ H Þ þ c
N ðrH þ mÞ2
!m kf k1
kf k1 :
ð31Þ
734
P. Huang et al. / Applied Mathematical Modelling 37 (2013) 728–741
Proof. Setting ðem ; gm Þ ¼ ðuh umh ; ph pmh Þ, subtracting (22) from (15) with
l ¼ h, we have
m m 1 m m Bh ððem ; gm Þ; ðv ; qÞÞ þ bðuh um H ; uh ; v Þ þ bðuH ; e ; v Þ þ m rhðaðe ; v Þ þ aðuH ; v ÞÞ ¼ 0:
ð32Þ
Taking ðv ; qÞ ¼ ðe ; g Þ in (32) and using (6), we arrive at m
m
m ðrh þ mÞkrem k0 6 c3 kAh uh k0 kuh um H k0 þ rhkruH k0
6 which yields
rh kf k þ c kA u k ðku uk0 þ ku uH k0 þ kuH um H k0 Þ; rH þ m 1 3 h h 0 h 2
m
!m
N
2
ðrh þ mÞkre k0 6 cðrh þ h þ H Þ þ c
ðrH þ mÞ2
kf k1
kf k1 :
ð33Þ
Furthermore, using (17), (18), (23), (32) and Theorem 4.3, we get
m m m m kgm k0 6 b1 Nkrum 2 H k0 kre k0 þ c 3 kuh uH k0 kAh uh k0 þ rhðkre k0 þ kruH k0 Þ !m N 2 2 kf k1 kf k1 : 6 cðrh þ h þ H Þ þ c ðrH þ mÞ2
ð34Þ
Combining (33) and (34) with Theorem 3.2, we finish the proof of the theorem. h Remark 4.1. Since it is valid that
1
m
kf k1
rH þ rh þ m r2 Hh kf k1 ¼ kf k > 0: ðrH þ mÞðrh þ mÞ mðrH þ mÞðrh þ mÞ 1
This result implies that the stability of Algorithm 4.1 is better than that of the general Galerkin finite element method. Remark 4.2. Setting
r ¼ 0 in (20)–(22), we get the general Oseen two-level scheme. By Theorems 4.2 and 4.4, we obtain
mkrðu umH Þk0 þ kp pmH k0 6 cH þ c
N
m
m
kf k1 2
kf k1 ;
and
mkrðu umh Þk0 þ kp pmh k0 6 cðh þ H2 Þ þ c
N
kf k1 2
m
m kf k1 :
If we choose H ¼ h, then we obtain the Oseen defect-correction method. Remark 4.3. From Theorem 3.2, for the usual stabilized finite element solution ðuh ; ph Þ, which involves solving one large Navier–Stokes problem on a fine mesh with mesh size h, we have the following error estimate:
mkrðu uh Þk0 þ kp ph k0 6 ch: Furthermore, if we choose H such that h ¼ OðH2 Þ for Algorithm 4.1, then we get the convergence rate of same order as the usual stabilized finite element method from Theorem 4.4. However, our method is more efficient than the one-level scheme. Algorithm 4.2 (m-Oseen-1-Simple) m Step I. Solve a defect Navier–Stokes problem on coarse mesh, i.e., find ðum H ; pH Þ 2 X H M H by (20). Step II. Solve a correction Navier–Stokes problem on fine mesh, i.e., find ðumh ; pmh Þ 2 X h M h such that for all ðv ; qÞ 2 X h M h ,
Bh ððumh ; pmh Þ; ðv ; qÞÞ þ
rh rh m m aðumh ; v Þ þ bðum aðuH ; v Þ: H ; uH ; v Þ ¼ ðf ; v Þ þ m m
ð35Þ
Next, we will study the stability and error estimate of Algorithm 4.2.
Theorem 4.5. Under the assumptions of Theorem 2.1, ðumh ; pmh Þ defined by scheme (35) satisfies
krumh k0 6
1 1 kf k1 ; þ rh þ m rH þ m
ð36Þ
P. Huang et al. / Applied Mathematical Modelling 37 (2013) 728–741
735
and the following error estimate:
N
2
mkrðu umh Þk0 þ kp pmh k0 6 cðh þ H Þ þ c
ðrH þ mÞ2
!m kf k1
kf k1 :
ð37Þ
Proof. Taking ðv ; qÞ ¼ ðumh ; pmh Þ in (35), using (7) and Theorem 4.1, we get
krumh k0 6 ðrh þ mÞ1 ð1 þ rhðrH þ mÞ1 þ NðrH þ mÞ2 kf k1 Þkf k1 6
r H þ r h þ 2m kf k1 ; ðrH þ mÞðrh þ mÞ
which is (36). Next, setting ðem ; gm Þ ¼ ðuh umh ; ph pmh Þ, we derive from (15) and (35) that m Bh ððem ; gm Þ; ðv ; qÞÞ þ m1 rhðaðem ; v Þ þ aðum H ; v ÞÞ þ bðuh uH ; uh ; v Þ þ bðuH ; uh uH ; v Þ m m þ bðuH um H ; uh ; v Þ þ bðuH ; uH uH ; v Þ ¼ 0;
8ðv ; qÞ 2 X h Mh :
ð38Þ
Setting ðv ; qÞ ¼ ðem ; gm Þ in (38), using (6) and (7) and Theorems 3.1 and 4.1, we obtain m m m ðrh þ mÞkrem k0 6 c3 kuh uH k0 ðkAh uh k0 þ kAH um H k0 Þ þ rhkruH k0 þ Nðkruh k0 þ kruH k0 ÞkrðuH uH Þk0 !m N 2 6 cðrh þ h þ H2 Þ þ c kf k1 kf k1 : ðrH þ mÞ2
ð39Þ
Moreover, it follows from (17) that
m m m m kgm k0 6 b1 c3 kuh uH k0 ðkAh uh k0 þ kAH um 2 H k0 Þ þ rhðkre k0 þ kruH k0 Þ þ Nðkruh k0 þ kruH k0 ÞkrðuH uH Þk0 !m N 2 2 kf k1 kf k1 : 6 cðrh þ h þ H Þ þ c ðrH þ mÞ2
ð40Þ
Hence, by combining (39) and (40) with Theorem 3.2, we have completed the proof of (37). The proof ends. h Algorithm 4.3 (m-Oseen-1-Newton) m Step I. Solve a defect Navier–Stokes problem on coarse mesh, i.e., find ðum H ; pH Þ 2 X H M H by (20). Step II. Solve a correction Navier–Stokes problem on fine mesh, i.e., find ðumh ; pmh Þ 2 X h M h such that for all ðv ; qÞ 2 X h M h ,
Bh ððumh ; pmh Þ; ðv ; qÞÞ þ
rh rh m m m aðumh ; v Þ þ bðumh ; um aðuH ; v Þ þ bðum H ; v Þ þ bðuH ; umh ; v Þ ¼ ðf ; v Þ þ H ; uH ; v Þ: m m
ð41Þ
Next, we will study the stability and error estimate of Algorithm 4.3. Theorem 4.6. Under the assumptions of Theorem 2.1, ðumh ; pmh Þ defined by scheme (41) satisfies
krumh k0 6
! 1 1 þ kf k1 ; Hm rH þ rrhþhmm rh þ rrHþ m
ð42Þ
and the following error estimate:
N
2
mkrðu umh Þk0 þ kp pmh k0 6 cðh þ H Þ þ c
ðrH þ mÞ2
!2m kf k1
kf k1 :
Proof. Taking ðv ; qÞ ¼ ðumh ; pmh Þ in (41) and using (7), we get m m 2 ðrh þ mÞkrumh k0 6 rhkrum H k0 þ Nkrumh k0 kruH k0 þ kf k1 þ NkruH k0 :
Using uniqueness condition and Theorem 4.1, we have
!
r2 Hh þ rHm þ rhm rh Nkf k1 kf k1 : krumh k0 6 þ1þ rH þ m rH þ m ðrH þ mÞ2 By some simple calculation, we have
ð43Þ
736
P. Huang et al. / Applied Mathematical Modelling 37 (2013) 728–741
krumh k0 6
rH þ rh þ 2m kf k ; r2 Hh þ rhm þ rHm 1
which is (42). Next, setting ðem ; gm Þ ¼ ðuh umh ; ph pmh Þ, we derive from (15) and (41) that m m m m Bh ððem ; gm Þ; ðv ; qÞÞ þ m1 rhaðem ; v Þ þ m1 rhaðum H ; v Þ þ bðe ; uh ; v Þ þ bðumh ; e ; v Þ þ bðumh uH ; umh uH ; v Þ ¼ 0:
ð44Þ
Setting ðv ; qÞ ¼ ðe ; g Þ in (44), using (6), (7), (9) and Theorems 3.1 and 4.1, we obtain m
m
m ðrh þ mÞð1 cÞkrem k20 6 Bh ððem ; gm Þ; ðem ; gm ÞÞ þ m1 rhaðem ; em Þ þ bðem ; um H;e Þ m m m 6 jbðuh uH ; uh uH ; em Þj þ jbðuh uH ; uH um H ; e Þj þ jbðuH uH ; uh uH ; e Þj m m 1 m m þ jbðuH um H ; uH uH ; e Þj þ jm rhaðuH ; e Þj 2 m 6 Nkrem k0 krðuh uH Þk20 þ krðuH um H Þk0 þ 2krðuH uH Þk0 krðuh uH Þk0 m þ rhkrum H k0 kre k0 :
ð45Þ
Moreover, it follows from (19), (23), (26) and Young inequality that
ðrh þ mÞkre k0 6 cðH þ rhÞ þ c
!2m
N
2
m
ðrH þ mÞ2
kf k1
kf k1 :
ð46Þ
Finally, one finds from (17) and (45) that
kgm k0 6 b1 2
rhðkrem k0 þ krumH k0 Þ þ Nkrem k0 krumH k0 þ Nðkrðuh uH Þk20 þ krðuH umH Þk20
þ 2krðuH um H Þk0 krðuh uH Þk0 Þ 6 cðH2 þ rhÞ þ c
N ðrH þ mÞ
!2m kf k1 2
kf k1 :
Hence, by combining this estimate and (46) with Theorem 3.2, we have completed the proof of (43). The proof ends.
h
Remark 4.4. From (30), (36) and (42), we find that the stability of Algorithm 4.1 is the best among these algorithms. h
Table 1 Comparisons of the one-level with the two-level defect-correction Oseen iterative methods. Methods
H
h
CPU-time
krðuuapp Þk0 kruk0
kppapp k0 kpk0
uH1 -Rate
pL2 -Rate
One-level Algorithm 4.1 Algorithm 4.2 Algorithm 4.3
– 1/4 1/4 1/4
1/16 1/16 1/16 1/16
0.454 0.235 0.204 0.219
0.249855 0.249881 0.249887 0.249826
0.0359607 0.0359257 0.0358945 0.0360782
– – – –
– – – –
One-level Algorithm 4.1 Algorithm 4.2 Algorithm 4.3
– 1/5 1/5 1/5
1/25 1/25 1/25 1/25
1.094 0.563 0.469 0.501
0.142012 0.142025 0.142018 0.142009
0.0165533 0.0165378 0.0165141 0.0166518
1.2659 1.2660 1.2661 1.2657
1.7384 1.7384 1.7396 1.7325
One-level Algorithm 4.1 Algorithm 4.2 Algorithm 4.3
– 1/6 1/6 1/6
1/36 1/36 1/36 1/36
2.329 1.141 0.954 1.100
0.0917065 0.0917145 0.0917115 0.0917041
0.00881857 0.00881067 0.00879852 0.00890549
1.1993 1.1993 1.1993 1.1993
1.7270 1.7269 1.7267 1.7163
One-level Algorithm 4.1 Algorithm 4.2 Algorithm 4.3
– 1/7 1/7 1/7
1/49 1/49 1/49 1/49
4.375 2.032 1.688 1.766
0.0643529 0.0643586 0.0643553 0.0643531
0.00520091 0.00519653 0.00519072 0.00528934
1.1489 1.1489 1.1490 1.1488
1.7127 1.7125 1.7117 1.6898
One-level Algorithm 4.1 Algorithm 4.2 Algorithm 4.3
– 1/8 1/8 1/8
1/64 1/64 1/64 1/64
7.469 3.422 2.875 2.985
0.0478061 0.0478105 0.0478077 0.0478079
0.00330493 0.00330241 0.00329995 0.00340521
1.1130 1.1129 1.1130 1.1128
1.6978 1.6975 1.6961 1.6490
One-level Algorithm 4.1 Algorithm 4.2 Algorithm 4.3
– 1/9 1/9 1/9
1/81 1/81 1/81 1/81
12.469 5.609 4.704 4.875
0.0370025 0.0370059 0.0370042 0.0370054
0.00222314 0.00222168 0.00222093 0.00234419
1.0875 1.0875 1.0874 1.0873
1.6832 1.6827 1.6810 1.5850
737
P. Huang et al. / Applied Mathematical Modelling 37 (2013) 728–741
Remark 4.5. If we choose H such that h ¼ OðH2 Þ for Algorithms 4.2 and 4.3, then we get the convergence rate of same order as the usual stabilized finite element method from Theorems 4.5 and 4.6. However, the given methods are more efficient than the one-level schemes. Remark 4.6. To the author’s knowledge, Algorithm 4.2 is the most efficient in these three methods for the small Reynolds number. Because the trilinear form is the most simple. Meanwhile, the coefficient matrices of Algorithms 4.2 and 4.3 are symmetric, while the one of Algorithm 4.1 is not.
5. Numerical experiments In this section we present numerical experiments to check the numerical theory developed in the previous sections and illustrate the efficiency of the two-level defect-correction Oseen iterative methods based on local Gauss integration. In the given experiments, the pressure and velocity are approximated by the lowest equal-order finite element pairs defined with respect to the same uniform triangulation. Meanwhile, the algorithms are implemented using public domain finite element software [27]. Besides, we take r ¼ 0:01 as in [14] and a ¼ 1. We consider the exact solution problem. Let X be the unit square in R2 . The exact solution for the velocity u ¼ ðu1 ; u2 Þ and the pressure p is given as follows:
pðx; yÞ ¼ 10ð2x 1Þð2y 1Þ; u1 ðx; yÞ ¼ 10x2 ðx 1Þ2 yðy 1Þð2y 1Þ; u2 ðx; yÞ ¼ 10xðx 1Þð2x 1Þy2 ðy 1Þ2 ; and the right-hand side f ¼ ðf1 ðx; yÞ; f2 ðx; yÞÞ is determined by original problem (1). Our goal in this test is to validate the properties and merits of the given two-level methods as compared with the onelevel method. Hence, we compare the solution generated at a fixed value of h for the one-level method with the results obtained using two-level methods where the finest grid has the same mesh spacing h. In all two-level computations, h is given, so H is chosen such that h ¼ H2 . Meanwhile, we have used the equal-order finite element P 1 P 1 with a fixed tolerance as 1
0.8
0.8
0.6
0.6 Y
Y
1
0.4
0.4
0.2
0.2
0
0
0.2
0.4
X
0.6
0.8
0
1
0.8
0.8
0.6
0.6
0.2
0.4
0
0.2
0.4
X
0.6
0.8
1
0.6
0.8
1
Y
1
Y
1
0
0.4
0.4
0.2
0.2
0
0
0.2
0.4
Fig. 1. Velocity streamlines for
X
0.6
0.8
1
0
X
m ¼ 1=3200: (a) standard Galerkin method, (b) Algorithm 4.2, (c) Algorithm 4.3, and (d) Algorithm 4.1.
738
P. Huang et al. / Applied Mathematical Modelling 37 (2013) 728–741
0.8
0.8
0.6
0.6 Y
1
Y
1
0.4
0.4
0.2
0.2
0
0
0.2
0.4
X
0.6
0.8
Fig. 2. Velocity streamlines for
0
1
0
0.2
0.4
X
0.6
0.8
1
m ¼ 1=5000: (a) Algorithm 4.2 and (b) Algorithm 4.1.
1
0.03
0.8
0.025 0.02 Y
Y
0.6
0.015
0.4 0.01 0.2
0
0.005
0
0.2
0.4
0.6
0.8
1
0
0.01
0.02
X Fig. 3. Velocity streamlines for
m ¼ 1=10000 by Algorithm 4.1: (a) whole and (b) detail in the left bottom corner.
1 0.9
0.03
X
0.6 Ghia et al. Algorithm 4.1
Ghia et al. Algorithm 4.1
0.4
0.8 0.2 v−velocity
0.7
y
0.6 0.5 0.4 0.3
0 −0.2 −0.4
0.2 −0.6
0.1 0 −0.5
0
u−velocity
0.5
1
−0.8
0
Fig. 4. Horizontal (a) and vertical velocity (b) for
0.2
0.4
x
0.6
0.8
1
m ¼ 1=3200.
1.0E5. We consider the case of the viscidity m ¼ 1, and pick six values of h, i.e., 1/16, 1/25, 1/36, 1/49, 1/64 and 1/81. The CPU-time, the L2 -error of the pressure, H1 -error of the velocity for the one-level method and the two-level methods for different values of h are tabulated in Table 1. From Table 1, we can see that the one-level method and two-level methods work well and keep the convergence rates just like the theoretical analysis. As expected, the two-level methods spend less time than one-level method under nearly the
739
P. Huang et al. / Applied Mathematical Modelling 37 (2013) 728–741
1
0.6 Ghia et al. Erturk et al. Algorithm 4.1
0.9 0.8
0.2 v−velocity
0.7
y
0.6 0.5 0.4 0.3
0 −0.2 −0.4
0.2
−0.6
0.1 0
Ghia et al. Erturk et al. Algorithm 4.1
0.4
−0.5
0
u−velocity
0.5
−0.8
1
0
Fig. 5. Horizontal (a) and vertical velocity (b) for
(a)
0.8
Ghia et al. Erturk et al. Algorithm 4.1
0.6
0.8
1
m ¼ 1=5000.
Ghia et al. Erturk et al. Algorithm 4.1
0.2 v−velocity
0.6 y
x
0.4
0.7
0.5 0.4 0.3
0 −0.2 −0.4
0.2
−0.6
0.1 0 −0.5
0.4
(b) 0.6
1 0.9
0.2
0
u−velocity
0.5
1
−0.8
0
Fig. 6. Horizontal (a) and vertical velocity (b) for
0.2
0.4
x
0.6
0.8
1
m ¼ 1=10000.
same relative error. Meanwhile, we can observe that Algorithms 4.2 and 4.3 are more efficient than Algorithm 4.1, which validate the analysis in Remark 4.6. However, from Remark 4.4, we know that the stability of Algorithm 4.1 is better than that of Algorithms 4.2 and 4.3. So, to show stability of Algorithm 4.1, we test a popular benchmark problem, the lid driven cavity, which has been analyzed in [28,29]. This problem is chosen because some benchmark data is available for comparison. Let the computation be carried out in the region X ¼ fðx; yÞj0 < x; y < 1g. We impose the normal component of the velocity to be zero on @ X and the tangential component to be zero except along y ¼ 1 where it is set to one. The exact solution of this problem is unknown. So, we take Ghia et al.’s algorithm [28] computed on a very fine mesh ðh ¼ 1=129Þ as the ‘‘exact’’ solution for the purpose of comparison. Meanwhile, we also give the results of the standard Galerkin method to show the good stability and efficiency of the two-level defect-correction Oseen iterative methods for the smaller viscidity problem. The computed results are shown, based on the standard Galerkin method, Algorithms 4.2 and 4.3, and Algorithm 4.1 for a set of different smaller viscidity coefficients m = 1/3200, 1/5000, 1/10000 with different mesh sizes h ¼ 1=64; 1=81; 1=121, respectively. In Figs. 1–3, we present the velocity streamlines of the cavity flow. The standard Galerkin method and Algorithm 4.3 are divergent despite m ¼ 1=5000. Algorithms 4.1 and 4.2 get almost the same results at this case. When the viscidity coefficient decreases to 1/10000, Algorithm 4.2 is divergent too. As we know, the position of the main vortex moves towards the center of the cavity when the viscidity decreases, and additional second vortex may appear in the right bottom corner of the cavity and a third vortex appears at the left bottom corner. When viscidity is 1/10000, six resolved vortices are captured. Above results for Algorithm 4.1 are consistent with those in [29] very well, which suggest its efficiency for solving the Navier–Stokes equations at high Reynolds numbers. In Figs. 4–6, we draw the x component of velocity along the vertical centerline and y component of velocity along the horizontal centerline based on Algorithm 4.1. Compared with the results obtained by Ghia et al. and Erturk et al., we can find that they are in excellent agreement. Moreover, the pressure contours of the cavity flow at different viscidity coefficients
740
P. Huang et al. / Applied Mathematical Modelling 37 (2013) 728–741
(a) 1
Y
0.6
0.4
0.2
0
0.2
0.4
X
0.6
0.8
0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.157104 0.120856 0.1 0.0802086 0.0546398 0.034689 0.0108654 0 -0.0316355 -0.058813 -0.0709235 -0.0784258
0.8
0.6 Y
0.8
0
(b) 1
0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.156319 0.131462 0.1 0.0824395 0.0751833 0.0681661 0.0533265 0.00204348 -0.029392 -0.0501363 -0.0622554 -0.0736008 -0.0814365
0.4
0.2
0
1
0
0.2
(c) 1
X
0.6
0.8
1
0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0698888 0.05293 0.051765 0.0411888 0.0403703 0.0378173 0.0339114 0.0292205 0.0180187 0 -0.0334256 -0.0533253 -0.0656891 -0.0760043
0.8
Y
0.6
0.4
0.2
0
0.4
0
Fig. 7. Pressure level lines for
0.2
0.4
X
0.6
0.8
1
m ¼ 1=3200 (a), m ¼ 1=5000 (b) and m ¼ 1=10000 (c) based on Algorithm 4.1.
are given in Fig. 7. From this figure, we can see that there is a singular point at the up right corner. Therefore, Algorithm 4.1 can capture this classical model well. 6. Conclusions In this work we have presented the two-level defect-correction Oseen iterative stabilized finite element methods in solving the steady Navier–Stokes equations based on local Gauss integration. The main feature of our method is to combine the defect-correction method and the two-level discretization. These methods are m-Oseen-1-Oseen, m-Oseen-1-Simple, and mOseen-1-Newton. The theoretical analysis and numerical tests show m-Oseen-1-Oseen is a favorite algorithm among these three algorithms for the considered problem. This idea can be extended to the time-dependent problems and other fluids dynamical equations. Acknowledgements The authors thank the editor and reviewers for their valuable comments and suggestions which helped us to improve the results of this paper. References [1] [2] [3] [4] [5] [6]
J. Xu, A novel two-grid method for semilinear elliptic equations, SIAM J. Sci. Comput. 15 (1994) 231–237. J. Xu, Two-grid discretization techniques for linear and nonlinear PDEs, SIAM J. Numer. Anal. 33 (1996) 1759–1778. W. Layton, A two level discretization method for the Navier–Stokes equations, Comput. Math. Appl. 26 (1993) 33–38. W. Layton, W. Lenferink, Two-level Picard and modified Picard methods for the Navier–Stokes equations, Appl. Math. Comput. 69 (1995) 263–274. W. Layton, L. Tobiska, A two-level method with backtracking for the Navier–Stokes equations, SIAM J. Numer. Anal. 35 (1998) 2035–2054. Y. Zhang, Y.N. He, A two-level finite element method for the stationary Navier–Stokes equations based on a stabilized local projection, Numer. Methods Part. Diff. Eq. 27 (2011) 460–477. [7] V. Ervin, W. Layton, J. Maubach, A posteriori error estimators for a two-level finite element method for the Navier–Stokes equations, Numer. Methods Part. Diff. Eq. 12 (1996) 333–346. [8] Y.N. He, K.T. Li, Two-level stabilized finite element methods for the steady Navier–Stokes problem, Computing 74 (2005) 337–351. [9] Y.N. He, A.W. Wang, A simplified two-level method for the steady Navier–Stokes equations, Comput. Methods Appl. Mech. Eng. 197 (2008) 1568–1576.
P. Huang et al. / Applied Mathematical Modelling 37 (2013) 728–741
741
[10] J. Li, Investigations on two kinds of two-level stabilized finite element methods for the stationary Navier–Stokes equations, Appl. Math. Comput. 182 (2006) 1470–1481. [11] Y.Q. Shang, A parallel two-level linearization method for incompressible flow problems, Appl. Math. Lett. 24 (2011) 364–369. [12] H.J. Stetter, The defect correction principle and discretization methods, Numer. Math. 29 (1978) 425–443. [13] S. Kaya, W. Layton, B. Riviere, Subgrid stabilized defect correction methods for the Navier–Stokes equations, SIAM J. Numer. Anal. 44 (2006) 1639– 1654. [14] W. Layton, H. Lee, J. Peterson, A defect-correction method for the incompressible Navier–Stokes equations, Appl. Math. Comput. 129 (2002) 1–19. [15] W. Layton, Solution algorithm for incompressible viscous flows at high Reynolds number, Vestnik Mosk. Gos. Univ. Ser. 15 (1996) 25–35. [16] O. Axelsson, W. Layton, Defect correction methods for convection dominated, convection diffusion equations, RAIRO J. Numer. Anal. 24 (1990) 423– 455. [17] Q.F. Liu, Y.R. Hou, A two-level defect-correction method for Navier–Stokes equations, Bull. Aust. Math. Soc. 81 (2010) 442–454. [18] K. Wang, A new defect correction method for the Navier–Stokes equations at high Reynolds numbers, Appl. Math. Comput. 216 (2010) 3252–3264. [19] P. Bochev, C.R. Dohrmann, M.D. Gunzburger, Stabilization of low-order mixed finite elements for the Stokes equations, SIAM J. Numer. Anal. 44 (2006) 82–101. [20] J. Li, Y.N. He, Z.X. Chen, A new stabilized finite element method for the transient Navier–Stokes equations, Comput. Methods Appl. Mech. Eng. 197 (2007) 22–35. [21] Y.N. He, J. Li, Convergence of three iterative methods based on the finite element discretization for the stationary Navier–Stokes equations, Comput. Methods Appl. Mech. Eng. 198 (2009) 1351–1359. [22] J.G. Heywood, R. Rannacher, Finite-element approximations of the nonstationary Navier–Stokes problem. Part I: Regularity of solutions and secondorder spatial discretization, SIAM J. Numer. Anal. 19 (1982) 275–311. [23] R.B. Kellogg, J.E. Osborn, A regularity result for the Stokes problem in a convex polygon, J. Funct. Anal. 21 (1976) 397–431. [24] S. Larsson, The long-time behavior of finite-element approximations of solutions to semilinear parabolic problems, SIAM J. Numer. Anal. 26 (1989) 348–365. [25] J. Li, Y.N. He, A stabilized finite element method based on two local Gauss integrations for the Stokes equations, J. Comput. Appl. Math. 214 (2008) 58– 65. [26] Y.N. He, A.W. Wang, L.Q. Mei, Stabilized finite-element method for the stationary Navier–Stokes equations, J. Eng. Math. 51 (2005) 367–380. [27] FreeFem++, version 2.23.