Local projection stabilized method on unsteady Navier–Stokes equations with high Reynolds number using equal order interpolation

Applied Mathematics and Computation 243 (2014) 465–481

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Local projection stabilized method on unsteady Navier–Stokes equations with high Reynolds number using equal order interpolation Gang Chen a, Minfu Feng a,⇑, Hong Zhou b a b

School of Mathematics, Sichuan University, Chengdu, China Exploration and Development Research Institute, PetroChina Southwest Oil & Gasfield Company, China

a r t i c l e

i n f o

Keywords: Unsteady Navier–Stokes equations High Reynolds number Local projection stabilized Pressure stability condition Crank–Nicolson method

a b s t r a c t In this paper we propose and analyze a stabilized method for unsteady Navier–Stokes equations with high Reynolds number, using local projection stabilized method to control spurious oscillations in the velocities due to dominant convection, or in the pressure due to the velocity–pressure coupling. Using equal-order conforming elements in space and Crank–Nicolson difference in time, we derive a fully discrete formulation. We prove stability and convergence of the approximate solution. The error estimates hold irrespective of the Reynolds number, provided the exact solution is smooth. This result is comparable with the streamline diffusion and continuous interior penalty methods. Ó 2014 Elsevier Inc. All rights reserved.

1. Introduction The stable and accurate mixed finite element methods (FEMs) for the Navier–Stokes equations (NSEs) may suffer from violating inf–sup stability condition and oscillating approximate solutions caused by high Reynolds number. The streamline diffusion (SD) method has been a popular method to tackle these two issues in the past two decades, due to its good stability and high accuracy. It was first proposed by Brooks and Hughes [1]. Johnson et al. [2] analyzed this method and extended it to time-dependent problems using time–space elements. Johnson et al. [3] also proposed an SD method on NSEs based on a streamfunction-vorticity formulation with divergence-free discrete velocities. Hansbo and Szepessy developed a velocity– pressure SD method using time–space elements for the incompressible NSEs [4]. More related work of stabilized methods for Stokes and NSEs can be found in [5–12]. When using SD method to solve time-dependent problems, one has to solve the discretization problem on d þ 1-dimensional space–time domain. Che Sun and his coworkers proposed and developed the finite difference streamline diffusion (FDSD) method [6,7], by using finite difference discrete in time, which only need to solve the discretization problem on d-dimensional space domain. This method not only reduces the computational work, but also keeps the good features of SD method. The FDSD method was applied to unsteady NSEs in [8,9]. However, the SD/FDSD methods have some undesirable features: they introduce addition nonphysical coupling terms between velocity and pressure; they produce inaccuracy numerical solutions near the boundary; they have to calculate second derivative when using high order elements. To overcome those disadvantages, alternative stabilized methods have been developed recently: the variational multiscale (VMS) methods [13–17], the orthogonal subscales methods [18], the continuous interior penalty (CIP) methods [19] ⇑ Corresponding author. E-mail address: [email protected] (M. Feng). http://dx.doi.org/10.1016/j.amc.2014.05.086 0096-3003/Ó 2014 Elsevier Inc. All rights reserved.

466

G. Chen et al. / Applied Mathematics and Computation 243 (2014) 465–481

and the local projection stabilize (LPS) methods [20–24]. These stabilized methods not only have good accuracy, but also avoid the undesirable features of SD/FDSD methods. Pressure projection [25,26] is also a favored method for increasing pressure stability. It was used to solve unsteady NSEs using equal-order elements successfully [27,28]. The subgrid scale eddy viscosity model is a numerical stabilization of convection dominated and underresolved flow. This approach adds an artificial viscosity only on the fine scales. Layton generalized this concept for the stationary convection diffusion problem [13], with error estimates almost comparable with the SD method. The subgrid scale eddy viscosity method was applied on unsteady NSEs [14,15] to derive error estimates dependent on a reduced Reynolds number (that means partially dependent on Reynolds number). This idea was used to derive a fully-discrete scheme using inf–sup unstable elements for unsteady NSEs [16,17]. Both orthogonal subscales and LPS methods involve some kinds of orthogonal property, which can be viewed as special cases of subgrid scale eddy viscosity method. The orthogonal subscales and LPS methods not only can be used to stabilize the velocity, but also can be used to stabilize pressure. In [18], orthogonal subscales method on Oseen problems was analysed, where the projection space is assumed to be continues. On the other hand, LPS method assumes the projection space to be discontinues, which was easier to be implemented in parallel computing. LPS method was first introduced to stabilize the Stokes problems [20]. Since then, several researches were done to develop the LPS method. LPS method on stationary Oseen equation was analysed in [21,22,24,23]. LPS method on stationary NSEs was also analyzed in [24], but the error estimates are dependent on Reynolds number. It is natural to wonder how to use LPS method to solve unsteady NSEs and to obtain comparable numerical results with SD/FESD [4,12,8] and CIP [19] methods. This motivates our research. The key question for the (stabilized) FEMs for NSEs with high Reynolds number is: how to derive the error estimates held irrespective of the Reynolds number. The SD [4,12,8] and CIP [19] methods deal with this point by adding nonlinear and jump stabilized terms, respectively. In this paper, we propose and analyze a Crank–Nicolson scheme to solve unsteady NSEs with high Reynolds number, using LPS method to stabilize both velocities and pressure. Unlike the SD/FDSD methods [4,12,8] and semi-discrete CIP method [19], there are no nonlinear or jump stabilized terms introduced in our scheme. This makes our method much easier to be implemented. For the initial data we use Ritz-projection instead of L2 -projection to avoid the inaccuracy pressure close to initial moment. The almost absolute stability and error estimates held irrespective of the Reynolds number are proved. With suitable choice of parameters, our method’s error estimates are quite comparable with the SD [4] and CIP [19] methods’ error estimates. We implement two numerical experiments to confirm and illustrate our theoretical analysis. An outline of the paper is as follows. In Section 2, we present necessary notations. In Section 3 we propose and analyze the stability of our method. In Section 4 we give error estimates for our scheme. In Section 5 we give some numerical experiments. In Section 6 we conclude the whole paper. Throughout this paper, we use C to denote a positive constant independent of Dt; h and m, not necessarily the same at each occurrence. 2. Basic notations Let X 2 Rd ðd ¼ 2; 3Þ be a bounded domain with polygonal or polyhedral boundary C ¼ @ X. Let W m;p ðXÞ; W m;p 0 ðXÞ denote the m-order Sobolev space on X; k  k and j  j denote the norm and semi-norm on these spaces. When p ¼ 2; Hm 0 ðXÞ ¼ m m;p W m;p ðXÞ and k  km ¼ k  km;p ; j  jm ¼ j  jm;p . We denote the inner product of Hm ðXÞ by ð; Þm and 0 ðXÞ; H ðXÞ ¼ W ð; Þ ¼ ð; Þ0 . Let X denote a Banach space, the mapping /ðx; tÞ : ½0; T ! X, and

k/kL2 ð0;T;XÞ ¼

Z

T 0

k/k2X ðtÞdt

1=2 ;

k/kL1 ð0;T;XÞ ¼ sup k/kX ðtÞ:

ð2:1Þ

06t6T

Vector analogs of the Sobolev spaces along with vector-valued functions are denoted by upper and lower case bold face font, respectively, e.g., H 10 ðXÞ; L2 ðXÞ and u. Let I ¼ ½0; T, where T is a fixed positive constant. The flow of an incompressible fluid is governed by the incompressible Navier–Stokes equations

8 ut þ u  r u  m D u þ r p ¼ f > > > < r  u ¼ 0 in X  I; > u ¼ 0 on C  I; > > : uðx; 0Þ ¼ u0 ðxÞ in X;

in X  I; ð2:2Þ

where u ¼ uðx; tÞ 2 Rd denotes the velocities, p ¼ pðx; tÞ 2 R denotes the pressure and f ¼ f ðx; tÞ 2 Rd denotes the body forces,

m ¼ Re1 denotes the viscosity coefficient, Re denotes the Reynolds number. Defining V ¼ H 10 ðXÞ; Q ¼ L20 ðXÞ :¼ L2 ðXÞ=R and Bðu; p; v ; qÞ ¼ mðru; rv Þ  ðr  v ; pÞ þ ðr  u; qÞ; bðw; u; v Þ ¼

1 1 ðw  ru; v Þ  ðw  rv ; uÞ: 2 2

ð2:3Þ ð2:4Þ

467

G. Chen et al. / Applied Mathematics and Computation 243 (2014) 465–481

Problem (2.2) can clearly be written in the form: find ðu; pÞ 2 V  Q , for all ðv ; qÞ 2 V  Q , such that

ðut ; v Þ þ Bðu; p; v ; qÞ þ bðu; u; v Þ ¼ ðf ; v Þ:

ð2:5Þ

Let Dt ¼ T=N be time step, where N is a positive integer, t n ¼ nDt; tnþ1=2 ¼ ðn þ 1=2ÞDt. We make use of the following notation

f n :¼ f ðx; t n Þ; f nþ1=2 :¼

un :¼ uðx; tn Þ; pn :¼ pðx; t n Þ;

f nþ1 þ f n ; 2

unþ1=2 :¼

ð2:6Þ

unþ1 þ un ; 2

ð2:7Þ

nþ1 ¼ pðx; t nþ1=2 Þ: p

ð2:8Þ

Then from (2.2), it is immediate to see that,

nþ1 ; v ; qÞ ¼ ðdt unþ1  ut ðx; tnþ1=2 Þ; v Þ þ mðrðunþ1=2  uðx; t nþ1=2 ÞÞ; rv Þ þ ðr  ðunþ1=2 ðdt unþ1 ; v Þ þ Bðunþ1=2 ; p  uðx; t nþ1=2 ÞÞ; qÞ þ ðbðunþ1=2 ; unþ1=2 ; v Þ  bðuðx; t nþ1=2 Þ; uðx; t nþ1=2 Þ; v ÞÞ þ ðf ðx; tnþ1=2 Þ; v Þ :¼ Rnþ1 ðv Þ þ ðf ðx; tnþ1=2 Þ; v Þ; nþ1

where dt u

¼ ðu

nþ1

ð2:9Þ

n

 u Þ=Dt. The following lemma is trivial.

Lemma 2.1. Suppose wtt ; wttt 2 L2 ð0; T; L2 ðXÞÞ, we have the estimates

kwt ðx; tnþ1=2 Þ  dt wnþ1 k20 6 C Dt3 kwttt k2L2 ðtn ;tnþ1 ;L2 ðXÞÞ ; kwðx; t nþ1=2 Þ 

wnþ1=2 k20

6 C Dt

3

ð2:10Þ

kwtt k2L2 ðtn ;tnþ1 ;L2 ðXÞÞ :

ð2:11Þ

Using Lemma 2.1 we can give an estimate for Rnþ1 ðv Þ. Lemma 2.2. We assume u is sufficiency smooth, then there holds

jRnþ1 ðv Þj 6 C Dt 3=2 kv k0 ð1 þ kunþ1=2 k1;1 Þ  ðkuttt kL2 ðtn ;tnþ1 ;L2 ðXÞÞ þ kutt kL2 ðtn ;tnþ1 ;H2 ðXÞÞ Þ:

ð2:12Þ

Let T h ¼ fKg be a quasi-uniform simplices, quadrilaterals or hexahedra mesh. For all K 2 T h , let hK be the diameter of K and h ¼ maxK2T h hK . Let Pkþ1 ðT h Þ represent the k þ 1-order continuous piecewise polynomial on decomposition T h ; P dc k ðT h Þ represent k-order discontinuous piecewise polynomial on decomposition T h , where k P 0 is an integer, when k ¼ 0, we can write b 1 b P dc 0 ðT h Þ as P 0 ðT h Þ ; let P k ðKÞ; k P 0 be the k-order polynomial on K; let P r ðT h Þ ¼ fv h 2 H ðXÞ : v h jK 2 P r ðKÞ; 8K 2 T h g, where j

K K P br ðKÞ ¼ P r ðKÞ þ Pdþ1 i¼1 ki P j1 ðKÞ; j ¼ 1; 2; ki ; i ¼ 1; 2; . . . ; d þ 1 are the d þ 1 barycentric coordinate functions of the element K.

Let Y h 2 H1 ðXÞ be a finite element space of continuous, piecewise polynomial function defined over T h satisfying: there is an interpolation operator ih : H1 ðXÞ ! Y h such that ih : H10 ðXÞ ! Y h \ H10 ðXÞ and rþ1

kw  ih wk0;K þ hjw  ih wj1;K 6 ChK kukrþ1;K :

ð2:13Þ 2

2

The finite element space are defined as V h ¼ Y h \ V and Q h ¼ Y h \ Q . Let Ph : L ðXÞ ! Y h be the L -projection, i.e. for any

v 2 L2 ðXÞ, find Ph v 2 Y h such that ðPh v ; wh Þ ¼ ðv ; wh Þ;

8wh 2 Y h :

ð2:14Þ

From the results in [19,29] we have the following result Lemma 2.3. When u 2 Hrþ1 ðXÞ, there holds rþ1

ku  Ph uk0 þ hju  Ph uj1 6 Ch

kukrþ1 ; rþ1

when u 2 W 1;1 ðXÞ: In addition, we also have the following approximate properties for Ph in [19], ku  Ph uk0;1 þ hju  Ph uj1;1 6 Ch

Lemma 2.4. For all u 2 W 1;1 ðXÞ and

kuk1;1 ;

ð2:15Þ ð2:16Þ

v h 2 Y h there holds

kPh ðuvh Þ  uvh k0 6 Chkuk1;1 kv h k0 ;

ð2:17Þ

kPh ðuvh Þ  Ph ðuÞv h k0 6 Chkuk1;1 kv h k0 ;

ð2:18Þ

kPh uk0;1 6 C P kuk0;1 ;

ð2:19Þ

kPh uk1;1 6 C P kuk1;1 :

ð2:20Þ

468

G. Chen et al. / Applied Mathematics and Computation 243 (2014) 465–481

3. Local projection stabilized method A Galerkin finite element discretization of (2.5) is unstable in the case of small viscosity or inf–sup unstable pairs. We propose an approach which uses LPS method to solve those two issues in this part. 3.1. LPS operators By a macro element M we denote the union of one or more neighboring cells K 2 T h . The diameter of a macro element M is denoted by hM . We assume that the decomposition of X into macro elements M 2 Mh is non-overlapping and also shape regular, in particular

hM  hK ;

8K  M; 8M 2 Mh 

ð3:1Þ

Note that the case ðMh ¼ T h Þ is accepted. Let Dh denote a discontinuous finite element space defined on Mh . We define Dh ðMÞ ¼ fv h jM : v h 2 Dh g and v h ðMÞ ¼ fv h jM : v h 2 v h ; g. We define mapping ph restricted to M as ph jM is the L2 -projection from L2 ðMÞ onto Dh ðMÞ; 8M 2 Mh . Define jh : L2 ðXÞ ! L2 ðXÞ as jh ¼ id  ph , where id is the identity. We assume Dh and jh satisfying the following property: Assumption 3.1. Dh and

kjh qk0;M 6

jh satisfy the following properties:

l ChM kqkl;M ;

8q 2 Hl ðMÞ; 8M 2 Mh ;

ð3:2Þ

where 0 6 l 6 r and l is a fixed integer. Assumption 3.2. Let the local inf–sup condition

9b > 0; 8h > 0; 8M 2 Mh :

inf

wh 2Dh ðMÞ;wh –0

ðv h ; wh ÞM

sup

v h 2v h ðMÞ;v h –0 kv h k0;M kwh k0;M

Pb>0

ð3:3Þ

R P be satisfied, where ðv h ; wh ÞM ¼ M2Mh M v h wh dx. The cases fulfillments of Assumptions 3.1 and 3.2 will be presented in Section 5. By (2.13), Assumption 3.2 and the technology in [22,30] we get the following lemma: Lemma 3.1. There exist two interpolation operators J h : Q ! Q h and I h : V ! V h satisfying the following orthogonally and approximation properties:

ðp  J h p; qh Þ ¼ 0; 8qh 2 Dh ; 8p 2 Q ; ðu  I h u; v h Þ ¼ 0; 8v h 2 Dh ; 8u 2 V; kp  J h pk0 þ hjp  J h pj1 6 Ch

rþ1

ku  I h uk0 þ hju  I h uj1 6 Ch

kpkrþ1 ;

rþ1

ð3:4Þ ð3:5Þ

8p 2 Q \ Hrþ1 ðXÞ; 8u 2 V \ H

kukrþ1 ;

rþ1

8u 2 V \ W

ku  I h uk0;1 þ hju  I h uj1;1 6 Chkuk1;1 ;

ðXÞ;

1;1

ðXÞ:

ð3:6Þ ð3:7Þ ð3:8Þ

3.2. Fully discrete stabilized scheme using LPS method In this subsection, we consider an fully discrete scheme using LPS method for (2.2). We use FEM in space and Crank– Nicolson difference in time. For given initial velocities u0h , the discretization read: nþ1 Find ðunþ1 ;p h h Þ 2 V h  Q h , for all ðv h ; qh Þ 2 V h  Q h such that nþ1=2

ðdt unþ1 h ; v h Þ þ bðuh

nþ1=2

; uh

nþ1=2

; v h Þ þ Bh ðuh

nþ1 ;p h ; v h ; qh Þ ¼ ðf ðx; t nþ1=2 Þ; v h Þ;

ð3:9Þ

where n ¼ 0; 1; 2; . . . ; N  1 and

unþ1 þ unh h ; 2 Bh ðu; p; v ; qÞ ¼ Bðu; p; v ; qÞ þ Sh ðu; p; v ; qÞ;

unþ1=2 ¼ h

Sh ðu; p; v ; qÞ ¼

X

a1;M ðjh ru; jh rv ÞM þ

M2Mh

ð3:10Þ ð3:11Þ X

a2;M ðjh rp; jh rqÞM

ð3:12Þ

M2Mh

nþ1 and unþ1 is the approximate solution of unþ1 ; p is the approximate solution of pðx; tnþ1=2 Þ. The symbols a1;M ; a2;M are posih h tive parameters satisfying

a1;M ¼ C 1 hmM ; a2;M ¼ C 2 hmM ; 8M 2 Mh ;

ð3:13Þ

G. Chen et al. / Applied Mathematics and Computation 243 (2014) 465–481

469

where C 1 ; C 2 are fixed positive constants to be chosen in calculation. m ¼ 1 or m ¼ 2. From (3.1) and (3.13), it’s obviously to see that

a1;M  hm ; a2;M  hm ; 8M 2 Mh :

ð3:14Þ

nþ1 nþ1 ¼ pðx; t nþ1=2 Þ. If pnþ1 Remark 3.1. p obtained by (3.9) is the approximation of p , the approximation of pnþ1 ¼ pðx; t nþ1 Þ, h h is required, we process the following step. nþ1 Let pnh þ pnþ1 ¼ 2p ; n ¼ 0; 1; . . . ; N  1, for any given pjh ; 0 6 j 6 N, we can get pih ; i ¼ 0; 1; 2; . . . ; N; i – j by solving h h nþ1 nþ1 n h ; n ¼ 0; 1; . . . ; N  1. When j ¼ 1, we give two different methods to get p1h : ph þ ph ¼ 2p (1) We can use the backwards stabilized method to get p1h ,

ðdt u1h ; v h Þ þ bðu1h ; u1h ; v h Þ þ Bh ðu1h ; p1h ; v h ; qh Þ ¼ ðf 1 ; v h Þ:

ð3:15Þ

(2) We can simply use (3.9) with n ¼ 0 and the time step Dt 0 ¼ 2Dt to get an approximate value of pðx; DtÞ (i.e. p1h ). nþ2 nþ1 (3) We can also get approximate value pnþ1 of pnþ1 by letting pnþ1 ¼ 12 ðp þp h h h h Þ; n ¼ 0; 1; . . . ; N  2 and 3 N 1  N1 Nh  pN1 pNh ¼ 2p p p ¼  . This skill is used in our numerical experiment in Section 5.2. h 2 h 2 h To avoid pressure oscillations in the transient solution close to initial time, we should choose the initial data as Ritz-projection onto a space of discretely divergence-free functions. The traditional approach of the initial data using L2 -projection is seeking u0h 2 V h , such that

ðu0h ; v h Þ ¼ ðu0 ; v h Þ;

8v h 2 V h :

ð3:16Þ

Numerical experiments in Section 5.3 shows this approach produces inaccuracy pressure close to initial moment. The similar subject was studied by Burman and Fernández [31] on transient Stokes equations, and he proposed to use Ritz-projection instead of L2 -projection to avoid possible inaccuracy solutions. In this paper, we use the idea of Burman and Fernández [31], proposing the following Ritz-projection to approach the initial data. The numerical experiments in Section 5.3 show the benefit of using Ritz-projection over L2 -projection. The Ritz-projection to approach the initial data is defined as Find ðu0h ; k0h Þ 2 V h  Q h , for all ðv h ; qh Þ 2 V h  Q h such that

ðu0h ; v h Þ þ Bh ðu0h ; k0h ; v h ; qh Þ ¼ ðu0 ; v h Þ þ Bðu0 ; 0; v h ; qh Þ:

ð3:17Þ

This problem is well-posed due to Theorem 3.3 and it’s easy to see that

ku0h k20 þ kðu0h ; k0h Þk2 6 ku0 k20 þ mkru0 k20 ; ku0h  u0 k20 þ kðu0h  u0 ; k0h Þk2 6 Cðh

2rþ2

ð3:18Þ 2r

þ mh þ h

2lþm

Þku0 k2kþ1 ;

ð3:19Þ

where 1=2

1=2

kðu; pÞk2 ¼ mkruk20 þ ka1;M jh ruk20 þ ka2;M jh rpk20 :

ð3:20Þ

3.3. Stability This subsection presents the analysis of stabilities. At first, we present the stability of velocities. nþ1 Theorem 3.2 (Velocities’ stability). Let ðunþ1 ;p Þ 2 V h  Q h denote the solution of (3.9), then for any arbitrary positive h h constant C 0 , when Dt 6 C 0 , we have the following result

max kuhiþ1 k20 þ 2Dt

06i6N1

  N1  X T iþ1 2  kf k2L1 ð0;T;L2 ðXÞÞ þ ku0 k20 þ mkru0 k20 : kðuiþ1=2 ; Þk 6 C exp p h h C 0  Dt i¼0

ð3:21Þ

nþ1 Proof. Testing (3.9) with ðv h ; qh Þ ¼ ðunþ1=2 ;p Þ gives h h

1 C0 1 nþ1=2 nþ1=2 nþ1=2 2 2 2 n 2 nþ1 ðkunþ1 ;p Þ 6 kf ðx; t nþ1=2 Þk20 þ ku k0 h k0  kuh k0 Þ þ kðuh h Þk ¼ ðf ðx; t nþ1=2 Þ; uh 2Dt 2C 0 h 2 C0 1 2 n 2 6 kf ðx; t nþ1=2 Þk20 þ ðkunþ1 h k0 þ kuh k0 Þ: 4C 0 2

ð3:22Þ

470

G. Chen et al. / Applied Mathematics and Computation 243 (2014) 465–481

Multiplying (3.22) with 2Dt and adding it from 0 to n gives 2 kunþ1 h k0 þ 2Dt

  n n n X X Dt Dt X Dt 2 i 2 hiþ1 Þk2 6 C 0 Dt kf ðx; t iþ1=2 Þk20 þ ku0h k20 kðuhiþ1=2 ; p kunþ1 k þ ku k þ 1 þ 2C 0 h 0 C 0 i¼1 h 0 2C 0 i¼0 i¼0 6 C 0 Dt

n nþ1 X Dt X 3 3 kf ðx; t iþ1=2 Þk20 þ kuih k20 þ ku0 k20 þ mkru0 k20 ; C 2 2 0 i¼0 i¼1

ð3:23Þ

when Dt < C 0 , by discrete Gronwall’s inequality we have 2 kunþ1 h k0 þ 2Dt

  n  X T iþ1=2 hiþ1 Þk2 6 C exp kf k2L1 ð0;T;L2 ðXÞÞ þ ku0 k20 þ mkru0 k20 : kðuh ; p C  D t 0 i¼0

ð3:24Þ



Now, we discuss the stability of pressure and the existence of unique solution. Theorem 3.3 (Pressure’s stability). For any ph 2 Q h , there exists a positive constant b independent of h; Dt and

m such that

ðr  v h ; ph Þ sup þ Chkjh rph k0 P bkph k0 : jv h j1 v h 2V h ;v –0 Proof. By the continuous inf–sup condition, for all ph 2 Q h  Q , there exists

ð3:25Þ

v 2V

such that

r  v ¼ ph ; jv j1 6 Ckph k0 :

ð3:26Þ

Using Green’s formula and (3.7), we can obtain

kph k20 ¼ ðr  v ; ph Þ ¼ ðr  ðv  I h v Þ; ph Þ þ ðr  I h v ; ph Þ ¼ ðv  I h v ; jh rph Þ þ ðr  I h v ; ph Þ 6 Chjv j1 kjh rph k0 þ ðr  I h v ; ph Þ;

ð3:27Þ

which implies (3.25). h Remark 3.4. In scheme (3.9), for given h; Dt, by Brouwer’s fixed point Theorem, the existence of unþ1 followed immediately. h nþ1 is followed by Babuska’s theorem and Theorem 3.3. Then for given h; Dt and unþ1 h , the existence of ph 4. Error estimates In this section, we assume that the solution ðu; pÞ of (2.2) satisfies the following regularity condition:

  max kukL1 ð0;T;W 1;1 ðXÞÞ ; kukL1 ð0;T;Hrþ1 ðXÞÞ ; kpkL1 ð0;T;Hrþ1 ðXÞÞ ; kutt kL2 ð0;T;H2 ðXÞÞ ; kuttt kL2 ð0;T;L2 ðXÞÞ ; kut kL2 ð0;T;Hrþ1 ðXÞÞ 6 C:

ð4:1Þ

We use the following definitions:

ðenu ; enp Þ :¼ ðun  unh ; pn  pnh Þ ¼ ðun  I h un ; pn  J h pn Þ þ ðI h un  unh ; J h pn  pnh Þ :¼ ðgn ; bn Þ þ ðnn ; cn Þ; ðeunþ1=2 ; enþ1 p Þ :¼

ð4:2Þ

!  n   n  eu þ enþ1 g þ gnþ1 nþ1 nn þ nnþ1 nþ1 nþ1 nþ1 u nþ1  p nþ1    : ¼ þ ;p ;  J ; J  p p p p h h h h 2 2 2

nþ1 Þ þ ðnnþ1=2 ; c nþ1 Þ: ¼ ðgnþ1=2 ; b

ð4:3Þ

2

Let p0 ðuÞ be the L -projection of u onto the space of piecewise constant functions with respect to the decomposition Mh , then one has the approximation properties

kp0 ðuÞk0;1 6 C p kuk0;1 ;

ð4:4Þ

ku  p0 ðuÞk0;1 6 Chkuk1;1 :

ð4:5Þ

nþ1 Theorem 4.1 (Velocities’ error estimates). Let ðu; pÞ 2 V  Q denote the solution of (2.2), ðunþ1 ;p Þ 2 V h  Q h denote the h h 1=2 1=2 solution of (3.9), when C 3 Dt 6 12 ; h 6 1; q1 ¼ C 1 C 2 C p C P 6 2  C24 and C 4 > 1 is an arbitrary fixed positive constant, there holds 2 max keiþ1 u k0 þ q0 Dt

06i6N1

N1 X 2r 2rþ2m 2rþm 2lþm 2 kðeiþ1=2 ; eiþ1 þh þh þ Dt 4 Þ; u p Þk 6 Cðmh þ h i¼0

ð4:6Þ

471

G. Chen et al. / Applied Mathematics and Computation 243 (2014) 465–481

where C 3 ¼ CkukL1 ð0;T;W 1;1 ðXÞÞ þ 12 and q0 ¼ 1  q21 P C14 . Proof. Subtracting (3.19) from (2.9) gives

nþ1 ; v h ; q Þ þ ðbðunþ1=2 ; unþ1=2 ; v h Þ nþ1 ; v h ; qh Þ ¼ ðdt gnþ1 ; v h Þ  Bh ðgnþ1=2 ; b ðdt nnþ1 ; v h Þ þ Bh ðnnþ1=2 ; c h h h nþ1 ; v h ; qh Þ þ Rnþ1 ðv h Þ:  bðunþ1=2 ; unþ1=2 ; v h ÞÞ þ Sh ðunþ1=2 ; p nþ1=2

Testing (4.7) with ðv h ; qh Þ ¼ ðn

ð4:7Þ

nþ1 Þ gives ;c

1 nþ1 ; nnþ1=2 ; c nþ1 Þk2 ¼ ðdt gnþ1 ; nnþ1=2 Þ  Bh ðgnþ1=2 ; b nþ1 Þ ðknnþ1 k20  knn k20 Þ þ kðnnþ1=2 ; c 2Dt ; uhnþ1=2 ; nnþ1=2 Þ  bðunþ1=2 ; unþ1=2 ; nnþ1=2 ÞÞ þ ðbðunþ1=2 h nþ1 Þ þ Rnþ1 ðnnþ1=2 Þ nþ1 ; nnþ1=2 ; c þ Sh ðunþ1=2 ; p :¼

5 X

Ri :

ð4:8Þ

i¼1

Using Cauchy–Schwarz inequality and Young’s inequality we can get

jR1 j 6 kdt gnþ1 k0 knnþ1=2 k0 6 Ch

2rþ2

1 Dt1 kut k2L2 ðtn ;tnþ1 ;Hrþ1 ðXÞÞ þ knnþ1=2 k20 ; 6

ð4:9Þ

Since q0 P C14 , there holds 0 < q1 0 6 C 4 . Combing (3.4)–(3.7), we obtain,

nþ1 Þj þ nþ1 Þj þ jðjh r  nnþ1=2 ; b jR2 j 6 mjðrgnþ1=2 ; rnnþ1=2 Þj þ jðgnþ1=2 ; jh rc þ

X

a1;M ðjh rgnþ1=2 ; jh rnnþ1=2 ÞM

M2Mh

X

a2;M ðjh rbnþ1 ; jh rcnþ1 ÞM

M2Mh 2r

2rþ2m

6 Cðmh þ h

þh

2rþm

q0

Þþ

8

nþ1 Þk2 : kðnnþ1=2 ; c

ð4:10Þ

Using triangle inequality, Green’s formula, r  u ¼ 0 and the fact bðw; v ; v Þ ¼ 0, we can bound R3 as

jR3 j 6 jbðnnþ1=2 ; nnþ1=2 ; gnþ1=2 Þj þ jbðgnþ1=2 ; gnþ1=2 ; nnþ1=2 Þj þ jbðunþ1=2 ; gnþ1=2 ; nnþ1=2 Þj þ jðnnþ1=2  runþ1=2 ; nnþ1=2 Þj  1 þ jðgnþ1=2  runþ1=2 ; nnþ1=2 Þj þ r  uhnþ1=2 ; unþ1=2 nnþ1=2 2 6 X R3;i : :¼

ð4:11Þ

i¼1

Using Green’s formula, Hölder’s inequality, inverse inequality, (3.8) and (3.7) we obtain

1 jR3;1 j 6 jðnnþ1=2  rnnþ1=2 ; gnþ1=2 Þj þ jðr  nnþ1=2 ; nnþ1=2 gnþ1=2 Þj 2 6 Cknnþ1=2 k0 krnnþ1=2 k0 kgnþ1=2 k0;1 6 CkukL1 ð0;T;W 1;1 ðXÞÞ knnþ1=2 k20 ;

ð4:12Þ

nþ1=2

jR3;2 j 6 Ckgnþ1=2 k0 krgnþ1=2 k0;1 kn 2rþ2

6 Ch Using

Green’s

k0

1 þ knnþ1=2 k20 : 6

formula,

ð4:13Þ

r  u ¼ 0,

triangle

inequality,

Hölder’s

inequality,

(3.8),

(4.4)

and

(4.5)

and

jh ðp0 ðunþ1=2 Þ  rnnþ1=2 Þ ¼ p0 ðunþ1=2 Þ  jh rnnþ1=2 we have jR3;3 j ¼ jðunþ1=2  rnnþ1=2 ; gnþ1=2 Þj 6 jððunþ1=2  p0 ðunþ1=2 ÞÞ  rnnþ1=2 ; gnþ1=2 Þj þ jðjh ðp0 ðunþ1=2 Þ  rnnþ1=2 Þ; gnþ1=2 Þj 6 Ch

rþ2

krnnþ1=2 k0 þ Ch

rþ1

kjh rnnþ1=2 k0 6 Cðh

2rþ2

þh

2rþ2m

1 q nþ1 Þk2 : Þ þ knnþ1=2 k20 þ 0 kðnnþ1=2 ; c 6 8

ð4:14Þ

472

G. Chen et al. / Applied Mathematics and Computation 243 (2014) 465–481

Using Hölder’s inequality, inverse inequality and (3.7), one gets

jR3;4 j 6 Cknnþ1=2 k20 krunþ1=2 k0;1 6 CkukL1 ð0;T;W 1;1 ðXÞÞ knnþ1=2 k20 ;

ð4:15Þ

jR3;5 j 6 Ckgnþ1=2 k0 krunþ1=2 k0;1 knnþ1=2 k0

R

rþ1

knnþ1=2 k0 1 2rþ2 6 Ch þ knnþ1=2 k20 : 6 6 Ch

ð4:16Þ

Ph ðunþ1=2 nnþ1=2 Þdx. Testing (3.9) with ðv h ; ph Þ ¼ ð0; Ph ðunþ1=2 nnþ1=2 Þ  C 0 Þ gives X nþ1=2 nþ1=2 ðr  uh ; Ph ðunþ1=2 nnþ1=2 ÞÞ þ a2;M ðjh rpnþ1 n ÞÞM ¼ 0: h ; jh rPh ðu

Let C 0 ¼

X

ð4:17Þ

M2Mh

Using (4.36) we estimate R3;6 as

jR3;6 j 6 jðr  nnþ1=2 ; unþ1=2 nnþ1=2  Ph ðunþ1=2 nnþ1=2 ÞÞj þ jðr  gnþ1=2 ; unþ1=2 nnþ1=2  Ph ðunþ1=2 nnþ1=2 ÞÞj X nþ1=2 nþ1=2 a2;M jðjh rpnþ1 n ÞÞM j þ h ; jh rPh ðu M2Mh 3 X R3;6;i : :¼

ð4:18Þ

i¼1

Using (2.17) and inverse inequality (3.8) we get

jR3;6;1 j 6 Chkunþ1=2 k1;1 krnnþ1=2 k0 knnþ1=2 k0 6 CkukL1 ð0;T;W 1;1 ðXÞÞ knnþ1=2 k20 ; jR3;6;2 j 6 Ch 6 Ch

rþ1

kn

ð4:19Þ

nþ1=2

k0 1 nþ1=2 2 þ kn k0 : 6

2rþ2

ð4:20Þ

Since

kjh rPh ðunþ1=2 nnþ1=2 ÞÞk0;M 6 kjh rðPh ðunþ1=2 nnþ1=2 Þ  Ph ðunþ1=2 Þnnþ1=2 Þk0;M þ kjh rðPh ðunþ1=2  p0 ðPh ðunþ1=2 ÞÞnnþ1=2 Þk0;M þ kjh rp0 ðPh ðunþ1=2 Þnnþ1=2 Þk0;M 6 CkukL1 ð0;T;W 1;1 ðXÞÞ knnþ1=2 k0;M þ C p C P kunþ1=2 k0;1 kjh rnnþ1=2 k0;M

ð4:21Þ

and l

nþ1 k0;M þ kjh rjh p nþ1 k0;M þ Ch kp nþ1 nþ1 k0;M 6 kjh rc nþ1 klþ1;M kjh rp h k0;M 6 kjh rc and q1 ¼

C 1=2 C 1=2 1 2 CpCP

we get

  a2;M kjh rcnþ1 k0;M þ Chl kpnþ1 klþ1;M

X

jR3;6;3 j 6

ð4:22Þ

M2Mh

   CkukL1 ð0;T;W 1;1 ðXÞÞ knnþ1=2 k0;M þ C p C P kunþ1=2 k0;1 kjh rnnþ1=2 k0;M 6 Ch

2lþm

þ

q0

m

q

m

nþ1 Þk2 þ Ch kuk2L1 ð0;T;W 1;1 ðXÞÞ knnþ1=2 k20 kðnnþ1=2 ; c 8 m nþ1 k0;M kjh rnnþ1=2 k0;M C 2 C p C P hM kjh rc

þ

X M2Mh 2lþm

nþ1 Þk2 þ Ch kuk2L1 ð0;T;W 1;1 ðXÞÞ knnþ1=2 k20 þ 0 kðnnþ1=2 ; c 8 X 1=2 1=2 m nþ1 k0;M kjh rnnþ1=2 k0;M þ q1 C 1 C 2 hM kjh rc

¼ Ch

M2Mh

6 Ch

2lþm

þ

q0 8

m

nþ1 Þk2 þ Ch kuk2L1 ð0;T;W 1;1 ðXÞÞ knnþ1=2 k20 þ kðnnþ1=2 ; c

q1 2

nþ1 Þk2 : kðnnþ1=2 ; c

ð4:23Þ

We continue to estimate R4 and R5 as 2rþm

q0

nþ1 Þk2 ; kðnnþ1=2 ; c   1 jR5 j 6 knnþ1=2 k20 þ C Dt 3 kuttt k2L2 ðtn ;tnþ1 ;L2 ðXÞÞ þ kutt k2L2 ðtn ;tnþ1 ;H2 ðXÞÞ : 6

jR4 j 6 Ch

þ

8

ð4:24Þ ð4:25Þ

G. Chen et al. / Applied Mathematics and Computation 243 (2014) 465–481

473

From (4.8) to (4.36), we now have

1 q nþ1 Þk2 6 Cðmh2r þ h2rþ2m þ h2rþm þ h2lþm Þ ðknnþ1 k20  knn k20 Þ þ 0 kðnnþ1=2 ; c 2Dt 2    1  nþ1 2 kn k0 þ knn k20 þ CkukL1 ð0;T;W 1;1 ðXÞÞ þ 2 þ C Dt 3 ðkuttt k2L2 ðtn ;tnþ1 ;L2 ðXÞÞ þ kutt k2L2 ðtn ;tnþ1 ;H2 ðXÞÞ Þ 2rþ2

þ Ch

Dt 1 kut k2L2 ðtn ;tnþ1 ;Hrþ1 ðXÞÞ :

ð4:26Þ

Multiplying (4.26) with 2Dt and adding it from 0 to n gives

knnþ1 k20 þ q0 Dt

n n X X iþ1 Þk2 6 Cðmh2r þ h2rþ2m þ h2rþm þ h2lþm þ Dt 4 Þ þ C 3 Dt kniþ1 k20 : kðniþ1=2 ; c i¼0

ð4:27Þ

i¼0

By discrete Gronwall’s inequality, when C 3 Dt 6 12 we have

knnþ1 k20 þ q0 Dt

n X iþ1 Þk2 6 Cðmh2r þ h2rþ2m þ h2rþm þ h2lþm þ Dt 4 Þ: kðniþ1=2 ; c

ð4:28Þ

i¼0

By triangle inequality we get our final conclusion. m

m

Remark 4.1. The choices of Dt; a1;M ¼ C 1 hM and a2;M ¼ C 2 hM can be done as follows. 2 (1) Since C p ! 1 and C P ! 1 when h ! 0, and q1 ¼ C 1=2 C 1=2 1 2 C p C P 6 2  C 4 , then we have C 2 < 4C 1 . So C 2 6 C 1 is good enough. (2) When Dt ! 0, C 3 Dt 6 12 can be fulfilled.

We define 2r

EðuÞ2 ¼ ðmh þ h

2rþ2m

þh

2rþm

þh

2lþm

þ Dt 4 Þ:

ð4:29Þ

Then we get

keunþ1=2 k20;1 6 2kgnþ1=2 k20;1 þ 2knnþ1=2 k20;1 6 Ch

2rþ2

d

d

þ Ch EðuÞ2 6 Ch EðuÞ2 :

ð4:30Þ

Theorem 4.2 (Pressure’s error estimates). Under the condition of Theorem 4.1, there holds

Dt

N1 N1 X X d 2 2 2 2 keiþ1 kdt eiþ1 p k 0 6 C Dt u k0 þ Cð1 þ h EðuÞ ÞEðuÞ : i¼0

ð4:31Þ

i¼0

Proof. By testing (4.7) with qh ¼ 0 and rearranging it we get nþ1 nþ1 nþ1=2 ðr  v h ; enþ1 ; unþ1=2 ; v h Þ  bðunþ1=2 ; uhnþ1=2 ; v h ÞÞ p Þ ¼ ðdt eu ; v h Þ  mðreu ; rv h Þ  ðbðu h X X a1;M ðjh renþ1=2 ; j h rv h ÞM þ a1;M ðjh runþ1=2 ; jh rv h ÞM þ Rnþ1 ðv h Þ  u M2Mh

M2Mh

6 X Ri : :¼

ð4:32Þ

i¼1

Bounding Ri term by term gets nþ1=2 jR1 þ R2 j 6 Ckdt enþ1 k0 jv h j1 u k0 jv h j1 þ mkreu

jR3 j 6 jðu

nþ1=2

 rv

ð4:33Þ

nþ1=2 Þj h ; eu

1  rv h ; enþ1=2 Þj þ jðr  unþ1=2 ; enþ1=2 v h Þj u u h 2 nþ1=2 nþ1=2 ;u ; v h Þjj þ jbðeu

þ

jðeunþ1=2

k0;1 Þkeunþ1=2 k0 jv h j1 6 ðkunþ1=2 k1;1 þ kenþ1=2 u 1 þ jðr  unþ1=2 ; enþ1=2 v h Þ :¼ R3;1 þ R3;2 : u h 2

ð4:34Þ

474

G. Chen et al. / Applied Mathematics and Computation 243 (2014) 465–481

We estimate R3;1 as d=2

jR3;1 j 6 Cð1 þ h

EðuÞÞEðuÞjv j1

ð4:35Þ

Similar to (4.36), we can get

ðr  uh ; Ph ðnnþ1=2 v h ÞÞ þ

X

nþ1=2 a2;M ðjh rpnþ1 v h ÞÞM ¼ 0: h ; jh rPh ðn

ð4:36Þ

M2Mh

Then we can estimate R3;2 as

1 1 jðr  eunþ1=2 ; gnþ1=2 v h Þ þ jðr  unþ1=2 ; nnþ1=2 v h Þ h 2 2 1 1 k0 kunþ1=2 k1;1 kv h k0 þ jðr  unþ1=2 ; nnþ1=2 v h  Ph ðnnþ1=2 v h ÞÞ þ jðr  uhnþ1=2 ; Ph ðnnþ1=2 v h ÞÞj 6 Chkrenþ1=2 u h 2 2 1 ¼ Chkreunþ1=2 k0 kunþ1=2 k1;1 kv h k0 þ jðr  eunþ1=2 ; nnþ1=2 v h  Ph ðnnþ1=2 v h ÞÞ 2 1 X nþ1 a2;M ðjh rph ; jh rPh ðnnþ1=2 v h ÞÞM j þ j 2 M2M

R3;2 6

h

6 Chkr

enþ1=2 k0 j u

6 Cð1 þ h

d=2

v h j1 þ Chkrenþ1=2 k0 knnþ1=2 k0 kv h k1;1 þ u

EðuÞÞEðuÞjv h j1 þ

X

1 X nþ1=2 j a2;M ðjh rpnþ1 v h ÞÞM j h ; jh rPh ðn 2 M2M h

nþ1

ða2;M kjh rc

mþl

k0;M þ h

Þðknk0;M þ kjh rnnþ1=2 kÞkv h k1;1

M2Mh

6 Cð1 þ h

d=2

EðuÞÞEðuÞjv h j1 þ Ch

d=2

nþ1 Þk0 EðuÞ þ kðnnþ1=2 ; c nþ1 Þk20 Þjv h j1 : ðkðnnþ1=2 ; c

ð4:37Þ

We continue to bound R4 and R5 1=2

1=2

nþ1 Þk0 jv h j1 : jR4 þ R5 j 6 Cka1;M jh renþ1=2 k0 jv h j1 þ Cka1;M jh runþ1=2 k0 jv h j1 6 kðnnþ1=2 ; c u

ð4:38Þ

And R6 is already bounded in Lemma 2.2. So there holds

Dt

N1 X iþ1 Þj2 jðr  v h ; c i¼0

jv

2 h j1

! 2

iþ1 k20 þ Ch kjh rc

6 C Dt

N1 X d 2 2 2 kdt eiþ1 u k0 þ Cð1 þ h EðuÞ ÞEðuÞ :

ð4:39Þ

i¼0

We get our final conclusion by Theorem 3.3 and triangle inequality.

h

2 Remark 4.2. A direct bound of kdt enþ1 u k0 is

2 kdt enþ1 u k0

6 C Dt

2

N X Dt keiu k20

!

2r

2rþ2m

6 C Dt 2 ðmh þ h

þh

2rþm

þh

2lþm

þ Dt4 Þ:

ð4:40Þ

i¼0

2rþ1 4

Remark 4.3. The choice m ¼ 1; l ¼ r; Dt  h 2 max keiþ1 u k0 6 Ch

2rþ1

06i6N1

Dt

;

N1 X 2rþ1 2 2 keiþ1 : p k0 6 Cð1 þ Dt Þh

;

m 6 h admits ð4:41Þ ð4:42Þ

i¼0

In this case, our results are comparable with the results obtained by SD [4] and CIP [19] methods. The choice r m ¼ 2; l ¼ r  1; Dt  h2 ; m 6 h admits 2r

2 max keiþ1 u k0 6 Ch ;

06i6N1

Dt

N 1 X 2rd 2r 2 2 keiþ1 þh Þh : p k0 6 Cð1 þ Dt

ð4:43Þ

ð4:44Þ

i¼0

5. Numerical experiments In this section, we calculate two standard benchmarks for NSEs: the Green–Taylor vortex and flow around a cylinder. In all experiments, the algorithms are implemented using public domain finite element software Freefem++ [32]. We give the

475

G. Chen et al. / Applied Mathematics and Computation 243 (2014) 465–481 Table 5.1 Available cases and convergence of max16i6N kui  uih k20 . Name

Yh

Dh

r

l

m

a1;M

a2;M

Convergence rate

Case I

P b1 ðT h Þ

P 0 ðT h Þ

1

1

1

0:1hM

0:1hM

h þ Dt4

Case II

P b2

2

1

2

hM

hM

h þ Dt4

1

ð

2

2

3 4

Case III

ðT h Þ P 1 ðT h Þ

P0 T h Þ P 0 ðT 2h Þ

1

1

1

0:1hM

0:1hM

h þ Dt4

Case IV

P 2 ðT h Þ

2

2

1

0:1hM

0:1hM

h þ Dt4

Case V

P b2

P dc 1 ðT 2h Þ P dc 1 ðT h Þ

2

2

1

0:1hM

0:1hM

h þ Dt4

2

ðT h Þ

3 5 5

cases satisfying Assumptions 3.1 and 3.2 in Table 5.1. All meshes are assumed to be simplicial meshes. In Cases III and IV, T h is the barycenter refine mesh of a regular macro simplicial mesh T 2h . By saying barycenter refine mesh we mean T h is constructed by connecting the barycenter of every element in T 2h and corresponding element’s vertices to get three new elements. 5.1. Convergence study Let X be the unit square in R2 ; T ¼ 1. We use the example which the true solution is:

8 > < u1 u2 > : p

¼  cosðpxÞ sinðpyÞ expð2p2 t=ReÞ; ¼ sinðpxÞ cosðpyÞ expð2p2 t=ReÞ;

ð5:1Þ

¼ 0:25ðcosð2pxÞ þ cosð2pyÞÞ expð4p2 t=ReÞ:

This problem is the Green–Taylor vortex used to study NSEs in [33–35]. In this subsection we compare the numerical performances of Case I, Case II with the standard Crank–Nicolson Galerkin methods:

Galerkin I : V h ¼ V \ P b1 ðT h Þ;

Q h ¼ Q \ P1 ðT h Þ;

Fig. 5.1. Computational mesh.

ð5:2Þ

476

G. Chen et al. / Applied Mathematics and Computation 243 (2014) 465–481

Galerkin II : V h ¼ V \ P 2 ðT h Þ;

Q h ¼ Q \ P1 ðT h Þ:

ð5:3Þ

9

We set Re ¼ 10 and Dt ¼ 0:01 for our calculation. The type of computational mesh is presented in Fig. 5.1. The results are presented in Tables 5.2, 5.3, 5.4, 5.5, where

EmaxðuÞ ¼ max kuiþ1  uiþ1 h k0 ;

ð5:4Þ

06i6N1

EHðuÞ ¼

Dt

N1 X

!12 juiþ1  uhiþ1 j21

ð5:5Þ

;

i¼0

Table 5.2 Case I. h

1

4 8 16 32

EmaxðuÞ

EHðuÞ

ELðpÞ

r maxðuÞ

r HðuÞ

r LðpÞ

Time

0.0775465 0.0171859 0.00406204 0.0010377

1.39936 0.712044 0.357117 0.178067

0.676914 0.0129245 0.00284987 0.000674466

2.17 2.08 1.97

0.97 1.00 1.00

5.71 2.18 2.08

4.96 24.73 186.6 2392

Table 5.3 Case II. h

1

4 8 16 32

EmaxðuÞ

EHðuÞ

ELðpÞ

r maxðuÞ

r HðuÞ

r LðpÞ

Time

0.0540041 0.0111074 0.0026441 0.000645558

1.14238 0.541242 0.264627 0.131321

0.0411578 0.00719825 0.00155375 0.000371964

2.28 2.07 2.03

1.08 1.03 1.01

2.52 2.21 2.06

10.60 44.32 330.6 3670

Table 5.4 Galerkin method with P1b–P1 elements. h

1

4 8 16 32

EmaxðuÞ

EHðuÞ

ELðpÞ

r maxðuÞ

r HðuÞ

r LðpÞ

Time

0.631635 0.450803 0.262054 0.136032

12.177 14.9602 16.3592 16.526

0.117277 0.0631518 0.024551 0.00731412

0.49 0.78 0.95

-0.30 -0.13 -0.01

0.89 1.36 1.75

4.560 17.43 95.03 843.0

Table 5.5 Galerkin method with P2–P1 elements. h

1

4 8 16 32

EmaxðuÞ

EHðuÞ

ELðpÞ

r maxðuÞ

r HðuÞ

r LðpÞ

Time

0.558528 0.269947 0.188294 0.293269

8.59607 6.92452 7.7833 18.994

0.142861 0.0360357 0.0166626 0.0281435

1.05 0.52 -0.64

0.31 -0.17 -1.29

1.99 1.11 -0.76

6.034 26.78 139.2 1117

Fig. 5.2. The triangulation of the computational domain.

477

G. Chen et al. / Applied Mathematics and Computation 243 (2014) 465–481 Table 5.6 Results maximal drag cd;max , maximal lift cl;max and Dpð8sÞ for different time step size by Case II.

Dt

tðcd;max Þ

cd;max

tðcl;max Þ

cl;max

Dpð8sÞ

Time

0.04 0.02 0.01 0.005 0.0025

3.93 3.94 3.94 3.94 3.9375

2.93993 2.94023 2.94032 2.94033 2.94034

5.88 5.76 5.71 5.705 5.7

0.405718 0.469050 0.484549 0.488832 0.489383

0.10495285 0.10582080 0.10508050 0.10708350 0.10800550

3418 5403 9871 18648 33418

Fig. 5.3. Streamlines of Case II at Dt ¼ 0:0025.

478

G. Chen et al. / Applied Mathematics and Computation 243 (2014) 465–481

ELðpÞ ¼

Dt

N1 X 2 iþ1  p iþ1 kp h k0

!12 ð5:6Þ

i¼0

and rmaxðuÞ ; r HðuÞ and rLðpÞ stands for corresponding convergence rates. The results show that our methods have much better numerical performances than the standard Galerkin method. The convergence rates are as we predicted in theoretical analysis. 5.2. Flow around a cylinder We consider the benchmark problem flow around a cylinder to test our numerical schemes. This is a well-known benchmark using in [36–38]. The domain with mesh (4206 elements) is presented in Fig. 5.2, where the origin of the cylinder is at ðx; yÞ ¼ ð0:2; 0:2Þ and the radius is equal to 0.15. The time-dependent inflow profile is

(

6 pt u1 ð0; y; tÞ ¼ u1 ð2:2; y; tÞ ¼ 0:41 2 sin 8 yð0:41  yÞ;

u2 ð0; y; tÞ ¼ u2 ð2:2; y; tÞ ¼ 0:

ð5:7Þ

Nonslip conditions are prescribed at the other boundaries. Computation is performed for Re ¼ 1000 and the external force f ¼ 0 using scheme Case II. The values for the maximal drag cd;max , maximal lift cl;max and Dpð8sÞ (here DpðtÞ ¼ pð0:15; 0:2; tÞ  pð0:25; 0:2; tÞ with different time step size Dt for Case II are presented in Table 5.6. The streamlines at t ¼ 2; 4; 5; 6; 7; 8 for Case II with Dt ¼ 0:0025 are plotted in Fig. 5.3. The references values [36] for cd;max , maximal lift cl;max and Dpð8sÞ are ref cref d;max 2 ½2:93:2:97; c l;max 2 ½0:47:0:49;

Dpð8sÞref 2 ½0:115; 0:105:

ð5:8Þ

From the figure we know the streamlines of our method totally agree with the results in [36–38]. The tables show cd;max ; cl;max and Dpð8sÞ are exactly in the interval of references values. 5.3. The behavior of pressure close to initial moment We use the examples in Section 5.1 and case II for comparison. We set Re ¼ 109 , Dt ¼ 10i ; i ¼ 1; 2; . . . ; 5 and h ¼ 1=32. The initial data is approximated by Ritz-projection (3.17) and L2 -projection (3.16), the pressure at t ¼ 5  106 ðDt ¼ 105 Þ, 1h computed by two choices of initial data is presented in Figs. 5.4 and 5.5, the true pressure p 1 ¼ pðx; 5  106 Þ is prei.e. p 1h  p 1 k0 with different Dt are presented in Table 5.7. sented in Fig. 5.6. The error of kp

Fig. 5.4. Streamlines of pressure using Ritz-projection (3.17) for initial data.

G. Chen et al. / Applied Mathematics and Computation 243 (2014) 465–481

479

Fig. 5.5. Streamlines of pressure using L2 -projection (3.16) for initial data.

Fig. 5.6. Streamlines of the true pressure.

From the table and figures we see, the pressure close to initial moment is much better when initial data approximated by Ritz-projection (3.17). And the pressure close to initial moment is very inaccuracy when the initial data approximated by L2 -projection (3.16).

480

G. Chen et al. / Applied Mathematics and Computation 243 (2014) 465–481 Table 5.7 Errors of pressure at time level 1=2 with two approaches of initial data. L2 -projection

Ritz-projection

1

10

0.000373388

0.000412867

102

0.000361085

0.000462451

103

0.000650940

0.000476617

104

0.005242410

0.000476902

105

0.052118400

0.000463951

Dt

6. Conclusions The fully discrete LPS method for time-dependent Navier–Stokes equations with high Reynolds number has been proposed and analyzed in this paper. We prove the almost absolute stability of this method and obtain error estimates independent of the viscosity coefficient. At certain choice of parameters, the error estimates are absolutely comparable with the SD [4] and CIP [19] methods. Acknowledgments This research is supported by the Natural Science Foundation of China (No. 11271273). The authors would like to thank the editor and referees for their criticism, valuable comments and suggestions which helped to improve this paper. References [1] A.N. Brooks, T.J.R. Hughes, Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Eng. 32 (1–3) (1982) 199–259. [2] C. Johnson, U. Nävert, J. Pitkäranta, Finite element methods for linear hyperbolic problems, Comput. Methods Appl. Mech. Eng. 45 (1–3) (1984) 285– 312. [3] C. Johnson, J. Saranen, Streamline diffusion methods for the incompressible Euler and Navier–Stokes equation, Math. Comput. 47 (175) (1986) 1–18. [4] P. Hansbo, A. Szepessy, A velocity pressure streamline diffusion finite-element method for the incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Eng. 84 (2) (1990) 175–192. [5] T. Zhou, M. Feng, A least squares Petrov–Galerkin finite element method for the stationary Navier–Stokes equations, Math. Comput. 60 (202) (1993) 531–543. [6] C. Sun, H. Shen, The finite difference streamline diffusion methods for time-dependent convection–diffusion equations, Numer. Math. Chin. Univ. Engl. Ser. 7 (1) (1998) 72–85. [7] Q. Zhang, C. Sun, Finite difference-streamline diffusion method for nonlinear convection–diffusion equation, Math. Numer. Sin. 20 (2) (1998) 211–224. [8] T. Sun, K. Ma, The finite difference streamline diffusion methods for the incompressible Navier–Stokes equations, Appl. Math. Comput. 149 (2) (2004) 493–505. [9] Q. Zhang, Finite difference streamline diffusion method for incompressible N–S equations, Math. Numer. Sin. 25 (3) (2003) 311–320. [10] G. Chen, M. Feng, A new absolutely stable simplified Galerkin least-squares finite element method using nonconforming element for the Stokes problem, Appl. Math. Comput. 219 (2013) 5356–5366. [11] G. Chen, M. Feng, Y. He, Unified analysis for stabilized methods of low-order mixed finite elements for stationary Navier–Stokes equations, Appl. Math. Mech. 34 (8) (2013) 953–970. [12] G. Chen, M. Feng, Y. He, Finite difference streamline diffusion method using nonconforming space for incompressible time-dependent Navier–Stokes equations, Appl. Math. Mech. 34 (9) (2013) 1083–1096. [13] W. Layton, A connection between subgrid scale eddy viscosity and mixed methods, Appl. Math. Comput. 133 (2002) 147–157. [14] V. John, S. Kaya, Finite element error analysis of a variational multiscale method for the Navier–Stokes equations, Adv. Comput. Math. 28 (2008) 43–61. [15] M.A. Belenli, S. Kaya, L.G. Rebholz, A subgrid stabilization finite element method for incompressible magnetohydrodynamics, Int. J. Comput. Math. 90 (7) (2013) 1506–1523. [16] M. Feng, Y. Bai, Y. He, Y. Qin, A new stabilized subgrid eddy viscosity method based on pressure projection and extrapolated trapezoidal rule for the transient Navier–Stokes equations, J. Comput. Math. 29 (4) (2011) 415–440. [17] J. Wen, M. Feng, Y. He, Convergence analysis of a new multiscale finite element method with the P1/P0 element for the incompressible flow, Comput. Methods Appl. Mech. Eng. 258 (1) (2013) 13–25. [18] R. Codina, Analysis of a stabilized finite element approximation of the Oseen equations using orthogonal subscales, Appl. Numer. Math. 58 (2008) 264– 283. [19] E. Burman, M.A. Fernandez, Continuous interior penalty finite element method for the time-dependent Navier–Stokes equations: space discretization and convergence, Numer. Math. 107 (2007) 39–77. [20] R. Becker, M. Braack, A finite element pressure gradient stabilization for the Stokes equations based on local projections, Calcolo 38 (2001) 173–199. [21] M. Braack, E. Burman, Local projection stabilization for the Oseen problem and its interpretation as a variational multiscale method, SIAM J. Numer. Anal. 43 (6) (2006) 2544–2566. [22] G. Matthies, P. Skrzypacz, L. Tobiska, A unified convergence analysis for local projection stabilisations applied to the Oseen problem, Math. Model. Numer. Anal. 41 (4) (2007) 713–742. [23] E. Burman, A. Linke, Stabilized finite element schemes for incompressible flow using Scott–Vogelius elements, Appl. Numer. Math. 58 (11) (2008) 1704–1719. [24] Y. Qin, M. Feng, K. Luo, K. Wu, Local projection stabilized finite element method for Navier–Stokes equations, Appl. Math. Mech. 31 (5) (2010) 651–664. [25] J. Li, Z. Chen, A new local stabilized nonconforming finite element method for the Stokes equations, Computing 82 (2008) 157–170. [26] P. Bochev, C. Dohrmann, M. Gunzburger, Stabilization of low-order mixed finite elements for the Stokes equations, SIAM J. Numer. Anal. 44 (2006) 82– 101. [27] J. Li, Y. He, Z. Chen, A new stabilized finite element method for the transient Navier–Stokes equations, Comput. Methods Appl. Mech. Eng. 197 (1–4) (2007) 22–35.

G. Chen et al. / Applied Mathematics and Computation 243 (2014) 465–481

481

[28] K. Wang, M. Feng, Y. He, Two-level stabilized finite element method for the transient Navier–Stokes equations, Int. J. Comput. Math. 87 (10) (2010) 2341–2360. [29] J.L. Lions, G. Papanicolaou, R.T. Rockafellar, The finite element method for elliptic problems [M], North-Holland Publishing Company, New York, 1978. [30] M. Crouzeix, V. Thomée, The stability in Lp and W 1p of the L2 -projection onto finite element function spaces, Math. Comput. 48 (178) (1987) 521–532. [31] E. Burman, M.A. Fernández, Galerkin finite element methods with symmetric pressure stabilization for the transient Stokes equations: stability and convergence analysis, SIAM J. Numer. Anal. 47 (1) (2008) 409–439. [32] F. Hecht, New development in Freefem++, J. Numer. Math. 20 (3–4) (2012) 251–265. [33] A.J. Chorin, Numerical solution for the Navier–Stokes equations, Math. Comput. 22 (1968) 745–762. [34] D. Tafti, Comparison of some upwind-biased high-order formulations with a second order central-difference scheme for time integration of the incompressible Navier–Stokes equations, Comput. Fluids 25 (1996) 647–665. [35] V. John, W. Layton, Analysis of numerical errors in large eddy simulation, SIAM J. Numer. Anal. 40 (2002) 995–1020. [36] M. Shafer, S. Turek, Benchmark computations of laminar flow around cylinder, in: Flow Simulation with High-Performance Computers II, Vieweg, 1996. [37] V. John, Reference values for drag and lift of a two-dimensional time-dependent flow around the cylinder, Int. J. Numer. Methods Fluids 44 (2004) 777– 788. [38] W. Layton, L.G. Rebholz, C. Trenchea, Modular nonlinear filter stabilization of methods for higher Reynolds numbers flow, J. Math. Fluid Mech. 14 (2012) 325–354.