A rational high-order compact difference method for the steady-state stream function–vorticity formulation of the Navier–Stokes equations

Computers and Mathematics with Applications (

)

–

Contents lists available at ScienceDirect

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

A rational high-order compact difference method for the steady-state stream function–vorticity formulation of the Navier–Stokes equations P.X. Yu a,b , Z.F. Tian a,c,∗ , Hongjie Zhang c a

Department of Mechanics and Engineering Science, Fudan University, Shanghai 200433, PR China

b

School of Mechanical Engineering, Shanghai Jiao Tong University, Shanghai 200240, PR China

c

Fusion Science and Technology Center, MAE Department, UCLA, Los Angeles, CA 90095, USA

article

info

Article history: Received 28 December 2014 Received in revised form 16 January 2017 Accepted 28 January 2017 Available online xxxx Keywords: Navier–Stokes equation Rational High-order compact scheme Stream function–vorticity Driven cavity problem

abstract A rational high-order compact (RHOC) finite difference (FD) method on the nine-point stencil is proposed for solving the steady-state two-dimensional Navier–Stokes equations in the stream function–vorticity form. The resulting system of algebra equations can be solved by using the point-successive over- or under-relaxation (SOR) iteration. Numerical experiments, involving two linear and two nonlinear problems with their analytical solutions and two flow problems including the lid driven cavity and backward-facing step flows, are carried out to validate the performance of the newly proposed method. Numerical solutions of the driven cavity problem with different grid mesh sizes (maximum being 513×513) for Reynolds numbers ranging from 0 to 17500 are obtained and compared with some of the accurate results available in the literature. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction The Navier–Stokes equations are highly nonlinear and are very difficult to solve, especially when the approximate solutions are required to have a high accuracy. Therefore, accurate and efficient difference representations of the Navier–Stokes equations are of vital importance. The main purpose of this paper is to present a new fourth-order compact finite difference (FD) scheme to solve the steady two-dimensional (2D) incompressible Navier–Stokes equations in the stream function–vorticity formulation. High order finite difference methods, which are classified as wide or compact type, have been used since they have some advantages when compared with the traditional second-order central differences (CD) methods or the first-order upwind differences (UD) methods [1–4]. The high order FD methods of wide type are obtained discretizating the equations by fourth order central differences which results in a wider computational stencil. The high-order compact type (HOC) FD methods are constructed with the aims of having a narrower stencil and maintaining the high order accuracy [5–22]. These methods, which are computationally efficient and yield highly accurate numerical solutions, can be obtained for the stream function–vorticity formulations of the Navier–Stokes equations [2,15,16,18,19,21–23,3]. Gupta et al. [18] applied their compact fourth-order nine-point compact FD formula developed for partial differential equations of elliptic type [14] to the solution of the steady 2D Navier–Stokes equations in the stream function–vorticity form and solved the driven cavity

∗

Corresponding author at: Department of Mechanics and Engineering Science, Fudan University, Shanghai 200433, PR China. E-mail address: [email protected] (Z.F. Tian).

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

2

P.X. Yu et al. / Computers and Mathematics with Applications (

)

–

problem for Reynolds numbers up to Re = 2000 with a maximum of 41 × 41 grid mesh using point-SOR type of iteration. This study was followed more recently by similar research by Spotz and Carey [15], and Li et al. [16]. Spotz and Carey [15] developed a compact fourth order method on a nine-point computational stencil for the steady 2D stream function–vorticity formulation of the Navier–Stokes equations. They solved the square driven cavity flow up to Re = 1000 with a maximum of 41 × 41 grid mesh using the generalized minimal residual method and successive approximations. Li et al. [16] presented a compact fourth order scheme for the stream function–vorticity formulations, which is strictly within the 3× 3 computational stencil and had a faster convergence than that of Gupta [18]. They solved the cavity flow with a grid size of 129 × 129 for Re ≤ 7500 using point-SOR iterations and found that SOR pointwise iteration does not work when Re ≥ 9000. In [23], Zhang employed fourth-order compact finite difference schemes [14] with multigrid techniques to simulate the two-dimensional square driven cavity flow for Reynolds numbers up to Re = 7500 with a maximum of 129 × 129 grid mesh. Dennis and Hudson [21] developed the same scheme as in Ref. [14] using another approach. This method is a two-dimensional version of the methods of exponential type and uses the Numerov approximation. They solved the problem of natural convection in a square cavity with two vertical sidewalls maintained at different temperatures and obtained results up to Ra = 105 (Ra is the non-dimensional Rayleigh number). Kalita et al. [24] computed the laminar solution of the problem with a fourth-order accurate HOC scheme originally proposed by Spotz and Carey [15] for the stream function–vorticity form of the 2D steady Navier–Stokes equations. They used conjugate gradient and hybrid biconjugate gradient stabilized algorithms and solved the natural convection in a square cavity up to Ra = 107 . Recently, Erturk and Gökçöl [3] developed a new fourth-order compact formulation, whose difference with Refs. [15,16,18,21,22] is that the final form of the HOC formulation is in the same form of the Navier–Stokes equations such that any iterative numerical method used for Navier–Stokes equations, can be easily applied to this new HOC formulation, since the final form of the presented HOC formulation is in the same form with the Navier–Stokes equations. They used a fine grid mesh of 601 × 601, as it was suggested by Erturk et al. [25] in order to be able to compute and obtain a steady solution for high Reynolds numbers and solved the driven cavity flow for Reynolds number up to Re = 20 000. In this work, we introduce a new fourth-order compact (4OC) method. The main difference of the new 4OC method with the HOC methods proposed in [15,18,22] is in the way that the fourth-order compact scheme is obtained for the 2D model problems (i.e. convection–diffusion equations). The key property with the present 4OC scheme is that it allows the point-successive over- or under-relaxation (SOR) iteration for low-to-high Reynolds numbers. For the 2D model problems with constant convection coefficient, in order to determine the computational stencil coefficients appropriate for the desired order of accuracy, the new 4OC method integrates exactly any linear combination of the functions {1, x, y, xy, x2 , y2 , x2 y, xy2 , x3 , y3 , x2 y2 , x4 , y4 }. For the variable convection coefficient problem, the 4OC scheme is achieved based on the 4OC scheme proposed for the constant convection coefficient problems and a practical technique, named the reminder term modification approach [19]. As the basis of a discretization method for the incompressible, 2D, steady-state flow problems, the 4OC formulation newly proposed for the 2D convection–diffusion equation is applied to the stream function–vorticity formulation of the Navier–Stokes equations. To demonstrate the performance of this new formulation, two linear and two non-linear problems with their analytical solutions are computed at different grid mesh sizes and the lid driven cavity flow problem for Reynolds numbers up to Re = 17 500 with a maximum of 513 × 513 grid mesh are solved and compared with some of the accurate results accessible in the literature. In the next section, we first present a rational fourth-order compact (R4OC) scheme on a 3 × 3 stencil for the 2D convection–diffusion equation with constant convection coefficient; then we develop a nine-point R4OC scheme of the 2D problem with variable convection coefficient. Applications of the currently proposed 4OC formulation to the Navier–Stokes equations in the stream function–vorticity formulation follow in Section 3. Convergence studies to validate fourth-order accuracy of the proposed R4OC method for four problems with their analytical solutions are given and the numerical solutions of the lid driven cavity and backward-facing step flow problems are performed in Section 4. Finally, Section 5 is devoted to some concluding remarks. 2. Fourth-order compact FD method for model problem Consider the two-dimensional non-homogeneous convective diffusion model problem

− auxx − buyy + p(x, y)ux + q(x, y)uy = f (x, y)

(1)

where a and b are constants and p and q vary spatially, and f is a sufficiently smooth function with respect to x and y. This equation is also the principal equation for the 2D steady Navier–Stokes equations with constant viscosity in both the primitive variables formulation and the stream function–vorticity formulation. Let the step-length in the x-direction be h = 1/n and in the y-direction be k = 1/m, where n and m are the numbers of subdivisions in the x- and y-directions respectively. There are several strategies [22,2,14,15,19] to derive 4OC schemes for Eq. (1). Gupta et al. [14] derived a 4OC FD scheme based on the truncated Taylor series expansions. Their procedures give the approximate value of a function at a mesh point as a linear combination of the analytic solutions of the partial differential equation. The FD formula is obtained by collocation over a set of mesh points surrounding the given mesh point for which the difference formula is derived. The procedure to combine terms and to simplify formula is straightforward but extremely tedious, especially for derivation in higher dimensions. Another technique to obtain HOC formulae consists of considering a particular equation and

P.X. Yu et al. / Computers and Mathematics with Applications (

)

–

3

employs second-order central difference scheme repeatedly. The discretization continues by approximating the derivatives appearing in the second-order truncation error terms through the use of the original partial differential equation itself until a desired approximation order is reached. This method was outlined originally by MacKinnon [22]. Based on MacKinnon’s method [22], a 9-point 4OC difference scheme for Eq. (1) at a mesh point (xi , yj ) is given by

[−Aδx2 − Bδy2 + P δx + Q δy + E δx δy + Gδy2 δx + H δx2 δy + K δy2 δx2 ]ui,j = F

(2)

where δx and δy are the first and δx2 and δy2 are the second-order central difference operators along x- and y-directions, respectively, and the coefficients A, B, C , D, E , G, H, and F are as follows: A=a+

h2  pij

 − 2δx pi,j ,

B=b+

k2  qi,j

 − 2δy qi,j

12 a 12 b  k2  q  h2  pi,j i,j 2 δx − δx − δy − δy2 pi,j P = 1− 12 a 12 b



 Q = 1− E=−

 k2  q  h2  pi,j i,j δx − δx2 − δy − δy2 qi,j 12 a 12 b 

pi,j qi,j

h2

k2

+ 12 a b   pi,j bh2 G= + k2 , 12

 F = 1−

a

h2  pi,j 12

a



1

+ (h2 δx qi,j + k2 δy pi,j ) 6   qi,j ak2 H = + h2 , 12

K =−

b

1 12

(bh2 + ak2 )

 k2  q  i,j δx − δx2 − δy − δy2 fi,j . 12

(3)

b

If a = b = 1 and k = h, the 4OC scheme (2) with coefficients (3) is the same as one given in [15], which has been used originally by Spotz and Carey for the steady stream-function vorticity formulation on the 2D Navier–Stokes equations. Using alternative techniques, Gupta et al. [14] have put forward a similar 4OC FD scheme for Eq. (1). The idea behind their methods is to utilize the Taylor series expansions of all the functions involved in Eq. (1). All these formulas may seem to be in different forms, but there is a clear indication that they are mathematically equivalent. These discretization schemes have been shown by analytical results and numerical experiments to have better numerical stability and higher accuracy approximations than the standard second-order finite difference schemes [26,27,15]. When p and q are constants, the coefficients A, B, C , D, E , G, H, and F in Eq. (2) are reduced to P = p,

 A=a 1+

Pe2x 12

E = α1 q + β1 p,



Q = q ,

,

B=b 1+ G = β2 p − α1 b,

Pe2y



12

H = α2 q − β1 a,

, K = −α2 b − β2 a,

F = 1 + α1 δx + α2 δx2 + β1 δy + β2 δy fi,j



 2

where Pex = ph/a and Pey = qk/b are called as the cell Reynolds number, and α1 = − Pe12x h , β1 = −

β2 =

k2

(4) Pey k 12

, α2 =

h2 12

and

. 12 It is easy to find that the 4OC scheme (2) or (4) becomes singular for pure convective problems (a = b = 0). This means the application of the 4OC scheme (2) or (4) may be restricted at high cell Reynolds numbers. In order to circumvent the flaws associated with the above schemes, we will introduce a new method of the construction of 4OC difference schemes for Eq. (1). 2.1. 4OC difference method: constant coefficients case In this subsection, we use a new approach to derive a 4OC FD scheme for the model problem (1) with constant coefficients. Consider a finite difference approximation for Eq. (1) at a grid point (xi , yj ) as

(−αh δx2 − αk δy2 + P δx + Q δy + E δx δy + Gδy2 δx + H δx2 δy + K δy2 δx2 )ui,j   = 1 + α1 δx + α2 δx2 + β1 δy + β2 δy2 fi,j

(5)

where δx , δy , δx2 and δy2 are as defined in (2) and fi,j , δx fi,j , δx2 fi,j , δy fi,j , δy2 fi,j are the source function and its first- and second-order central difference approximations along x- and y-directions, respectively, at the nodal point (xi , yj ). In order to determine

4

P.X. Yu et al. / Computers and Mathematics with Applications (

)

–

the parameters αh , αk , P , Q , E , G, H , K , α1 , α2 , β1 and β2 , rewrite Eq. (5) as

(−αh δx2 − αk δy2 + P δx + Q δy + E δx δy + Gδy2 δx + H δx2 δy + K δy2 δx2 )ui,j    = 1 + α1 δx + α2 δx2 + β1 δy + β2 δy2 −auxx − buyy + pux + quy i,j .

(6)

Making (6) exact for 1, x, y, xy, x2 , y2 , x2 y, xy2 , x3 , y3 , x2 y2 , x4 and y4 results in a local linear system of twelve equations for the twelve unknowns A, B, P , Q , E , G, H , K , α1 , α2 , β1 and β2 . The system is solved and the parameters are obtained as P = p,

Q = q,

 αh = a  α1 = −

α2 =

1+

Pe2x 6

1+

Pe2x 12

ph2 12a

h2 12

·

·





,

Pe2 y



1+ 6

αk = b

,

Pe2 y

1+ 12

1 1+

1+

Pe2x 6

1+

Pe2x 12

Pe2x 12

2

,

qk β1 = − 12b ·

Pe2 y

,

1+ 12

,

E = α1 q + β1 p,

1



k2 12

β2 = G = β2 p − α1 b,

Pe2 y

·

H = α2 q − β1 a,

1+ 6

Pe2 y

1+ 12

K = −α2 b − β2 a,

(7)

where Pex = ph/a and Pey = qk/b. Eq. (5) with coefficients (7) is a nine-point compact FD scheme of O(h4 + k4 ) for the 2D convection–diffusion model problem (1) with constant convection coefficients. The scheme (5) may be named as the rational high-order compact (RHOC) FD scheme; i.e., the influence coefficients of the FD formulation are connected to rational functions of the coefficients of the differential operator and mesh size. It is interesting to note that the RHOC scheme (5) for the model Eq. (1) becomes the standard fourth-order Padé-like scheme for the pure convective problems (a = b = 0) or for the pure diffusive problems (p = q = 0). For example, when a = b = 0, Eq. (1) reduces to the pure convection equation (without loss of generality, here we assume p and q are constants): pux + quy = f (x, y)

(8)

and Eq. (6) reduces to

  p 1+

k2 6

      h2 h2 k2 δy2 δx + q 1 + δx2 δy ui,j = 1 + δx2 + δy2 fi,j 6

6

6

(9)

i.e.,



pδx 1+

h2 6

δx2

+



qδy 1+

k2 6

ui,j = fi,j −

δy2

h2 k2 36



1+

h2 6

δx2 δy2 

δx2

1+

k2 6

δy2

 fi,j

(10)

which, neglecting the term O(h2 k2 ), is equivalent to the following standard fourth-order Padé scheme for the pure convection equation (8):



pδx 1+

h2 6

δx2

+



qδy 1+

k2 6

δy2

ui,j = fi,j .

(11)

This fact shows that, for larger Reynolds numbers or cell Reynolds numbers, the present scheme (5) is more effective than the fourth-order compact scheme (2). In Section 4, we shall present the numerical results that verify this inference. 2.2. 4OC FD method: variable coefficients case In this subsection, an O(h4 + k4 ) compact rational scheme for the model problem (1) with variable convection coefficients will be developed by using the technique of remainder term modification proposed in [19]. Consider the FD scheme for Eq. (1) at a mesh point (xi , yj ) as

[−α¯ h δx2 − α¯ k δy2 + pi,j δx + qi,j δy + E¯ δx δy + G¯ δy2 δx + H¯ δx2 δy + K¯ δy2 δx2 ]ui,j = F¯

(12)

P.X. Yu et al. / Computers and Mathematics with Applications (

)

–

5

where

 α¯ h = a  α¯ 1 = −

α¯ 2 =

2

1+

12a

12

·



,

2

1+ pi,j h2

h2



Pex 6 Pex 12

α¯ k = b 

1

·

1+

β¯ 1 = −

2 Pex

E¯ = α¯ 1 q + β¯ 1 p,

12b

Pey

 ,

12

1

·

2

Pey

1+

,

12

2

2 Pex

6

6 2

1+ qi,j k2

Pex 12

1+

1+

, 2

2

Pey

1+

,

k2

β¯ 2 =

12

12

·

1+

Pey 6

2

Pey

1 + 12 ¯ = α¯ 2 q − β¯ 1 a, H

¯ = β¯ 2 p − α¯ 1 b, G

, K¯ = −α¯ 2 b − β¯ 2 a,

F¯ = 1 + α¯ 1 δx + α¯ 2 δx2 + β¯ 1 δy + β¯ 2 δy fi,j

 2



(13)

in which Pex = pi,j h/a and Pey = qi,j k/b. By using the Taylor series expansions and the original differential equation (1), the modified partial differential equation corresponding to the scheme (12) is obtained

− auxx − buyy + p(x, y)ux + q(x, y)uy − (2α¯ 2 qx + 2β¯ 2 py )uxy − (2α¯ 2 px )uxx − (2β¯ 2 qy )uyy − (α¯ 1 px + α¯ 2 pxx + β¯ 1 py + β¯ 2 pyy )ux − (α¯ 1 qx + α¯ 2 qxx + β¯ 1 qy + β¯ 2 qyy )uy = f + O(h4 + h2 k2 + k4 )

(14) h2

Pey k 12

k2

where α¯ 1 , α¯ 2 , β¯ 1 and β¯ 2 are given by (13), and α¯ 1 = − Pe12x h + O(h4 ), α¯ 2 = 12 + O(h4 ), β¯ 1 = − + O(k4 ) and β¯ 2 = 12 + O(k4 ). Eq. (14) shows that the local truncation error of Eq. (12) is only second-order accurate. To obtain a fourth-order compact difference scheme, we use the technique of remainder term modification [19]. Adding the term (2α¯ 2 qx + 2β¯ 2 py )uxy + (2α¯ 2 px )uxx + (2β¯ 2 qy )uyy + (α¯ 1 px + α¯ 2 pxx + β¯ 1 py + β¯ 2 pyy )ux + (α¯ 1 qx + α¯ 2 qxx + β¯ 1 qy + β¯ 2 qyy )uy to the left-hand side of Eq. (1) yields

− Af uxx − Bf uyy + Pf ux + Qf uy = Fp

(15)

where Af = a − 2α¯ 2 px ,

Pf = p + α¯ 1 px + α¯ 2 pxx + β¯ 1 py + β¯ 2 pyy ,

Bf = b − 2β¯ 2 qy ,

Qf = q + α¯ 1 qx + α¯ 2 qxx + β¯ 1 qy + β¯ 2 qyy ,

Fp = f − (2α¯ 2 qx + 2β¯ 2 py )uxy .

(16)

Eqs. (16) and (1) are the same in the form. Being similar to the scheme (5), we can obtain a fourth-order compact rational FD scheme for Eq. (1) with variable convection coefficients

  −A˜ δx2 − B˜ δy2 + Pf δx + Qf δy + E˜ δx δy + G˜ δx δy2 + H˜ δx2 δy + K˜ δx2 δy2 ui,j = F˜i,j

(17)

in which

 A˜ = Af 

α˜ 1 = −

α˜ 2 =

2

1+

2

1+ Pf h2 12Af

h2 12

˜x Pe 6

·

˜x Pe 12

 , 1

·

1+ 1+



1+

B˜ = Bf 

, 2

˜x Pe 12

β˜ 1 = −

˜ 2x Pe 6

˜ 2x Pe

,

β˜ 2 =

12Bf

12



˜y Pe

,

6 2

1+ Qf k2

k2

12

2

1+

·

˜y Pe 12

1

·

1+

2

1+

,

12

˜ 2y Pe 6 2

1+

˜y Pe

˜y Pe

,

12

E˜ = α˜ 1 Qf + β˜ 1 Pf + 2(α¯ 2 qx + β¯ 2 py ),

˜ = β˜ 2 Pf − α˜ 1 Bf , G

˜ = α˜ 2 Qf − β˜ 1 Af , H

F˜ = f + α˜ 1 fx + α˜ 2 fxx + β˜ 1 fy + β˜ 2 fyy ,

K˜ = −α˜ 2 Bf − β˜ 2 Af , (18)

˜ x = Pf h/Af and Pe ˜ y = Qf k/Bf . In (17) and (18), α¯ 2 , β¯ 2 and Af , Bf , Pf , Qf are given by (13) and (16) respectively, while where Pe Af , Bf , Pf , Qf , qx and py were written in short for (Af )i,j , (Bf )i,j , (Pf )i,j , (Qf )i,j , (qx )i,j and (py )i,j . Scheme (17) with coefficients (18) is an O(h4 + k4 ) compact rational FD scheme for the 2D model convective diffusion equation (1) on the nine-point 2D stencil.

6

P.X. Yu et al. / Computers and Mathematics with Applications (

)

–

3. Stream function–vorticity (ψ − ω) formulation and discretization In non-dimensional form, the steady Navier–Stokes equations in the primitive variable formulation set governing twodimensional flow of viscous incompressible fluid are written as

 1 uxx + uyy + uux + v uy = −px + f1 Re  1 − vxx + vyy + uvx + vvy = −py + f2 Re ux + v y = 0 −

(19) (20) (21)

where u and v are the velocities along the x- and y-directions respectively, Re is the nondimensional Reynolds number, p is the pressure and f = (f1 , f2 ) is the non-dimensional forcing function. The difficulty in numerically solving the steady incompressible Navier–Stokes equations (19)–(21) in primitive variable form arises from the specific velocity–pressure coupling. Hence, in order to avoid the difficulty with pressure, an alternative formulation in stream function–vorticity form has been used by introducing ω = vx − uy and u = ψy , v = −ψx

−ωxx − ωyy + Re uωx + Re vωy = f¯

(22)

−ψxx − ψyy = ω

(23)

where f¯ = f1y − f2x . The fourth-order compact (4OC) scheme for the stream function Eq. (23) can be given with u = ψ , f = ω, a = b = 1 and p = q = 0 in Eqs. (5) and (7). This 4OC FD formulation is

     h2 k2  2 2 k2 h2 − δx2 + δy2 + + δx δy ψij = 1 + δx2 + δy2 ωij . 12

12

12

12

(24)

The vorticity equation (22)is of 2D convection–diffusion type, and the R4OC approximation in this case may be obtained with u = ω, f = f¯ , a = b = 1 and p = Re u, q = Re v in Eqs. (13) and (16)–(18). We note that, in Eqs. (16) and (18), the value of coefficients involves u, ux , uy , uxx , uyy , v , vx , vy , vxx and vyy , where ux (=ψxy ), uy (=ψyy ), uxx (=ψxxy ), vx (= − ψxx ), vy (= − ψxy ) and vyy (= − ψxyy ) can be approximated by the standard second-order central formulae. For a fully 4OC rational scheme for (22), the discretization of u, v with O(h4 + k4 ) and uyy (=ψyyy ), vxx (= − ψxxx ) with O(h2 + k2 ) on nine-node computational stencil must be considered, which is done as follows:

(uyy )i,j = (ψyyy )i,j = −(ωy + ψxxy )i,j   = − δy ω + δx2 δy ψ i,j + O(h2 + k2 )

(25)

and

(vxx )i,j = −(ψxxx )i,j = (ωx + ψxyy )i,j   = δx ω + δx δy2 ψ i,j + O(h2 + k2 )

(26)

while the fourth-order approximation formulations of u(=ψy ) and v(= − ψx ) at a grid point (xi , yj ) can be obtained in the following way in [19] (also see [18] if k = h): ui,j = δy ψi,j +

k2  6

δy ω + δx2 δy ψ

+ O(k4 + h2 k2 )



i ,j

(27)

and

vi,j = −δx ψi,j −

h2  6

δx ω + δx δy2 ψ



i ,j

+ O(h4 + h2 k2 ).

(28)

To validate the currently proposed R4OC method for the Navier–Stokes equations in stream function–vorticity formulation for incompressible viscous flows, we will apply it to compute a problem with an analytical solution firstly, and then obtain the numerical solutions of the lid driven cavity flow problem in Section 4. 4. Numerical experiments In this section, we perform numerical experiments to illustrate the validity and effectiveness of the R4OC schemes developed in this article. We consider six problems including two linear and two nonlinear problems which have analytical solutions, and two flow problems including the lid driven cavity and backward-facing step flow problems which have the benchmark numerical solutions. As the first four problems have their analytical solutions, Dirichlet boundary conditions are

P.X. Yu et al. / Computers and Mathematics with Applications (

)

–

7

Table 1 Comparisons of the errors and the convergence rate for different ε , Problem 1. Grid size

4OC

R4OC

L2 -error

Rate

L∞ -error

Rate

L2 -error

Rate

L∞ -error

Rate

ε=1 11 × 11 21 × 21 41 × 41 81 × 81

1.260e−5 8.229e−7 5.273e−8 3.673e−9

3.937 3.964 3.844

2.737e−5 1.704e−6 1.069e−7 7.311e−9

4.006 3.995 3.870

1.643e−5 1.072e−6 6.866e−8 4.687e−9

3.938 3.965 3.873

3.450e−5 2.179e−6 1.362e−7 9.176e−9

3.985 4.000 3.892

ε = 0.1 11 × 11 21 × 21 41 × 41 81 × 81

2.081e−4 1.349e−5 8.606e−7 5.434e−8

3.947 3.970 3.985

4.846e−4 3.092e−5 1.927e−6 1.202e−7

3.970 4.004 4.003

3.397e−5 2.263e−6 1.455e−7 9.236e−9

3.908 3.959 3.978

7.734e−5 4.885e−6 3.078e−7 1.926e−8

3.985 3.988 3.998

ε = 0.01 11 × 11 21 × 21 41 × 41 81 × 81

2.764e−3 2.147e−4 1.406e−5 8.900e−7

3.686 3.933 3.982

6.40e−3 4.953e−4 3.181e−5 1.992e−6

3.693 3.961 3.997

5.460e−5 3.588e−6 2.322e−7 1.480e−8

3.928 3.950 3.971

1.282e−4 8.104e−6 5.130e−7 3.224e−8

3.984 3.982 3.992

ε = 0.001 11 × 11 21 × 21 41 × 41 81 × 81

5.594e−3 1.263e−3 1.330e−4 9.195e−6

2.147 3.247 3.854

1.260e−2 2.841e−3 2.990e−4 2.058e−5

2.149 3.248 3.861

6.101e−5 3.967e−6 2.534e−7 1.603e−8

3.943 3.969 3.983

1.483e−4 9.469e−6 5.949e−7 3.725e−8

3.969 3.992 3.997

ε = 0.0001 11 × 11 21 × 21 41 × 41 81 × 81

5.681e−3 1.562e−3 3.931e−4 6.704e−5

1.863 1.990 2.552

1.287e−2 3.495e−3 8.854e−4 1.502e−4

1.881 1.981 2.559

6.140e−5 4.014e−6 2.567e−7 1.623e−8

3.935 3.967 3.983

1.492e−4 9.714e−6 6.186e−7 3.885e−8

3.941 3.973 3.993

used, while for the flow problem, both Dirichlet (for ψ ) and Neumann (for ψx , ψy ) boundary conditions are applied. For all the problems, a point iterative under-relaxation procedure (i.e. relaxation factor λ ≤ 1) is associated with the solution of the resulting system of algebra equations. If the numerical results were divergence, the relaxation factor which was set as 1 at the beginning, would be reduced until the convergent results were obtained. Comparisons are made between analytical solutions and numerical results for the currently proposed R4OC FD methods, as well as some previously published HOC and other approximation methods. All computations are run by our own solver using Fortran language. 4.1. Problem 1 Consider the 2D constant coefficient model equation

ε uxx + εuyy + ux + uy = f (x, y) 0 ≤ x, y ≤ 1.

(29)

The analytical solution of this equation is u(x, y) = sin π x + sin π y.

(30)

The source function f (x, y) can be obtained according to Eq. (29). Eq. (29) is numerically solved for max |um+1 − um | ≤ 10−12 at the grid points (here m is the iterative count) using the 4OC polynomial FD scheme (4) and the present R4OC FD scheme (17) on 2D uniform grids. Comparisons of computational results are listed in Table 1 for the coefficient ε ranging from 1.0 to 0.0001. The results show that the present method can obtain very accurate results and keep the theoretical convergence rate even if the diffusion coefficient is very small. However, we notice that the convergence rate of the 4OC polynomial scheme decreases as the diffusion coefficient decreases. In addition, the computational costs of the two methods are also listed in Table 2. Although the R4OC method would need a small relaxation factor λ to obtain the convergence results for small ε , its number of iterations can be saved greatly especially for ε = 0.0001. All these prove that the present R4OC method has the advantage for solving the convection-dominated diffusion problem. 4.2. Problem 2 In this subsection, we consider the variable coefficient model equation as follows, uxx + uyy + Re(p(x, y)ux + q(x, y)uy ) = f (x, y) 0 ≤ x, y ≤ 1 where the p(x, y) and q(x, y) are the functions of (x, y), which are given by p(x, y) = (1 − x2 )(2y − 1) q(x, y) = 2xy(y − 1).

(31)

8

P.X. Yu et al. / Computers and Mathematics with Applications (

)

–

Table 2 Comparisons of the number of iterations for different ε , Problem 1. Grid size

4OC

R4OC

λ

Iterations

λ

Iterations

ε=1 11 × 11 21 × 21 41 × 41 81 × 81

1.0 1.0 1.0 1.0

303 1 143 4 312 16 625

1.0 1.0 1.0 1.0

303 1 143 4 312 16 626

ε = 0.1 11 × 11 21 × 21 41 × 41 81 × 81

1.0 1.0 1.0 1.0

107 386 1 460 5 556

1.0 1.0 1.0 1.0

106 386 1 461 5 557

ε = 0.01 11 × 11 21 × 21 41 × 41 81 × 81

1.0 1.0 1.0 1.0

106 64 79 243

0.7 1.0 1.0 1.0

88 38 60 240

ε = 0.001 11 × 11 21 × 21 41 × 41 81 × 81

1.0 1.0 1.0 1.0

325 1 012 841 303

0.09 0.1 0.3 0.5

351 321 99 60

ε = 0.0001 11 × 11 21 × 21 41 × 41 81 × 81

1.0 1.0 1.0 1.0

334 1 425 5 497 11 117

0.01 0.01 0.05 0.07

3 189 3 310 1 084 459

Table 3 Maximal absolute errors and the convergence rate, Problem 2. Re

Grid size

4OC [14]

R4OC

u-error

Rate

u-error

Rate

10

65 × 65 129 × 129 257 × 257

1.9814e−3 1.2418e−4 7.8184e−6

– 3.996 3.989

1.9650e−3 1.2309e−4 7.7673e−6

– 3.997 3.986

104

65 × 65 129 × 129 257 × 257

2.1220e−1 2.4687e−2 2.0131e−3

– 3.104 3.616

2.1278e−1 8.9248e−3 6.0008e−4

– 4.575 3.895

The analytical solution of this equation is u(x, y) = sin π x + sin 13π x + cos π y + cos 13π y.

(32)

The source function f (x, y) also can be obtained according to Eq. (31). This problem, which has been considered by Zhang [28], is solved with the 4OC polynomial FD [14] and the present 4OC rational FD methods (17). The iterative procedure is started with zero initial data and the computations were stopped when the maximum absolute difference between two successive iteration steps was smaller than 10−10 . The maximal absolute pointwise errors and the convergence rate are listed in Table 3 for Reynolds numbers Re = 10 and 104 . Notice that the approximate solution from the present 4OC rational FD method is more accurate than that from the 4OC polynomial FD method [14] especially for larger Re case. It can be seen that the present HOC scheme attains its theoretical accuracy orders for both Reynolds numbers. However, we also notice that HOC polynomial scheme did not approximate its theoretical accuracy orders for larger Reynolds numbers when it used coarse grid size. For example, if we use 65 × 65 grid to solve this problem for Re = 104 , the HOC polynomial scheme gives the poorer results and the rate will fall to 3.104. This shows that the present scheme has the ability to solve the convection-dominated problem with coarser grid size. 4.3. Problem 3 Consider the 2D model equation for fluid flow uxx + uyy − uux − v uy = (− sin x − 2 sin y) cos x

vxx + vyy − uvx − vvy = (− sin y + 2 sin x) cos y 0 ≤ x, y ≤ π .

(33)

P.X. Yu et al. / Computers and Mathematics with Applications (

)

–

9

Table 4 Solutions of Problem 3 with h = π/20 for u(0.7π, y). y/π

Analytic solution

Chen 4OC scheme [33]

Pillai 4OC scheme [32]

Dennis 4OC scheme [21]

Present 4OC scheme (18)

0.1 0.2 0.3 0.4 0.5

0.1816356 0.3454915 0.4755283 0.5590170 0.5877852

0.1816368 0.3454934 0.4755308 0.5590201 0.5877886

0.1816368 0.3454931 0.4755299 0.5590185 0.5877868

0.1816368 0.3454932 0.4755299 0.5590185 0.5877867

0.1816366 0.3454930 0.4755301 0.5590190 0.5877873

Table 5 Absolute errors and the convergence rate, Problem 3. y/π

h = π/10

h = π/20

Rate

0.1 0.2 0.3 0.4 0.5

1.6083e−5 2.5093e−5 3.0586e−5 3.3816e−5 3.4914e−5

9.7163e−7 1.5033e−6 1.8187e−6 2.0012e−6 2.0630e−6

4.049 4.061 4.072 4.079 4.081

Fig. 1. Schematic of the vortex centers in the lid-driven cavity, Problem 5.

The analytical solution of this equation is u(x, y) = − cos x sin y,

v(x, y) = sin x cos y.

(34)

This problem, which is a modification of the one considered by Roscoe [29] as a test problem for the developed methods, has been used by Dennis et al. [21,30,31], Pillai [32] and Chen et al. [33] to demonstrate and verify the methods proposed by them. The terms on the right-hand side of Eq. (33) do not constitute all the possible pressure gradients to be truly representative of the 2D Navier–Stokes equations, but Eq. (34) satisfies the equation of continuity. So, the essential features of a model solution of the Navier–Stokes equations are preserved by this problem. Eq. (33) is solved for max |um+1 − um | ≤ 10−10 (here m is the iterative count) using the present 4OC rational FD scheme y π (17) on 2D grids with a uniform resolution of h = k = 20 . The computed results for u(0.7, y) for values of π from 0.1 to 0.5 are listed in Table 4. The results obtained by the 4OC exponential FD schemes of Chen et al. [33] and Pillai [32] are also shown in Table 4. In addition, absolute errors and the rate of the convergence of scheme (18) are given in Table 5. The rate of the convergence is computed by using the log2 (eπ/10 /eπ /20 ), where e = |uexact − ucomput | and the subscript such as π /10 stands for the grid size. We note that the accuracy order of the scheme (17) is close to 4, validating the fourth-order accuracy of the current compact rational FD scheme in this problem. 4.4. Problem 4 In this subsection, we consider the recirculating viscous flow in a square cavity driven (0 ≤ x, y ≤ 1) by combined shear and body forces. The velocity boundary conditions are no-slip conditions with u(x, 0) = u(0, y) = u(1, y) = v(x, 0) = v(x, 1) = v(0, y) = v(1, y) = 0 and u(x, 1) = 16x2 (x − 1)2 . The problem is chosen such that its analytical solution

10

P.X. Yu et al. / Computers and Mathematics with Applications (

)

–

Fig. 2. Convergence history of the residual errors for different Reynolds numbers on a mesh of 65 × 65 grids, Problem 5.

Fig. 3. Comparisons of the u-velocity along the vertical centerline and the v -velocity along the horizontal centerline at (a) Re = 100, (b) Re = 400, (c) Re = 1000 and (d) Re = 3200 with 65 × 65 grids.

is known which makes it possible to assess the method thoroughly. Following [34], the test problem, which has an exact solution for the Navier–Stokes equations (22)–(23), is given in Ω = (0, 1) × (0, 1)

ψ(x, y) = 8g1 (x)g2 (y)   ω(x, y) = −8 g1′′ (x)g2 (y) + g1 (x)g2′′ (y) u(x, y) = 8g1 (x)g2′ (y)

v(x, y) = −8g1′ (x)g2 (y)

(35)

P.X. Yu et al. / Computers and Mathematics with Applications (

)

–

11

Fig. 4. Comparisons of the u-velocity along the vertical centerline and the v -velocity along the horizontal centerline at (a) Re = 3200, (b) Re = 5000, (c) Re = 7500 and (d) Re = 10 000 with 129 × 129 grids.

(4)

(4)

f¯ (x, y) = 8 g1 (x)g2 (y) + 2g1′′ (x)g2′′ (y) + g1 (x)g2 (y) + 64Re g1 (x)g2′ (y) g1′′′ (x)g2 (y) + g1′ (x)g2′′ (y)



  − g1′ (x)g2 (y) g1′′ (x)g2′ (y) + g1 (x)g2′′′ (y)







 (36)

where g1 (x) = x2 (x − 1)2 and g2 (y) = y2 (y2 − 1). This problem, which has been considered by Shih et al. [34], is solved with the 4OC polynomial FD [15] (also see (2)) and the present R4OC FD methods for the Navier–Stokes equations (22)–(23). The iterative procedure is started with zero initial data and the calculations are carried out when max |um+1 − um | ≤ 10−10 , where m is the iterative count. The L∞ -errors for ψ and ω are given in Table 6 for Reynolds numbers Re = 10, 100 and 1000. Notice that the approximate solution from the present R4OC FD method is more accurate than that from the 4OC polynomial FD method [15]. The estimated accuracy orders of the corresponding schemes are also listed in this table. It can be seen that the present HOC scheme attains its theoretical accuracy orders for all Reynolds numbers. However, for the HOC polynomial scheme, although their estimated accuracy order is also close to their theoretical accuracy orders in most cases, the absolute errors are always higher than those of the present R4OC method especially in the cases of larger Reynolds numbers. For example, if we use 41 × 41 grid to solve this problem for Re = 1000, the HOC polynomial scheme gives very poor results. This shows that the present scheme has the ability to solve the complicated flow problem on coarser grid size. 4.5. Problem 5 Now we consider the 2D lid-driven cavity flow problem, where a square domain [0, 1] × [0, 1] contains a recirculating flow induced by the sliding motion of the top wall (y = 1) from left to right. This problem is considered as the classical problem for the assessment of numerical methods and the validation of Navier–Stokes codes [35,36,15,16,18,37,23,3,25,38–41], particularly the steady-state solution of the incompressible fluid flows governed by the stream function–vorticity formulation of the Navier–Stokes equations (22)–(23). The boundary conditions are those of no slip: on the top wall u = 1 and v = 0, on all other walls u = 0 and v = 0, as shown in Fig. 1. Further the stream function values on all four walls are

12

P.X. Yu et al. / Computers and Mathematics with Applications (

)

–

Fig. 5. Comparisons of the u-velocity along the vertical centerline and the v -velocity along the horizontal centerline at (a) Re = 2500, (b) Re = 5000, (c) Re = 7500 and (d) Re = 10 000 with 129 × 129 grids.

Table 6 Comparisons of the L∞ -errors and the convergence rate for different Re, Problem 4. Grid size

4OC

R4OC

ψ

Rate

ω

Rate

ψ

Rate

ω

Rate

Re = 10 21 × 21 41 × 41 81 × 81

4.036e−6 2.623e−7 1.820e−8

3.944 3.849

9.002e−5 5.805e−6 5.258e−7

3.955 3.465

3.230e−7 2.347e−8 1.559e−9

3.783 3.912

1.008e−5 7.740e−7 5.168e−8

3.703 3.905

Re = 100 21 × 21 41 × 41 81 × 81

1.432e−4 9.327e−6 5.888e−7

3.940 3.986

4.967e−3 3.214e−4 2.030e−5

3.950 3.985

8.087e−5 7.120e−6 4.927e−7

3.506 3.853

4.081e−3 2.508e−4 1.717e−5

4.024 3.869

Re = 1000 41 × 41 81 × 81 161 × 161

3.205e−1 6.218e−5 3.641e−6

12.33 4.094

1.314e+1 3.018e−3 1.737e−4

12.09 4.119

3.322e−4 3.916e−5 2.795e−6

3.085 3.808

1.449e−2 1.659e−3 1.469e−4

3.127 3.497

zero (ψ = 0). At high Re, several secondary and tertiary vortices begin to appear, whose characteristics depend on Re. In Fig. 1, the abbreviations BL, BR and TL refer to bottom left, bottom right and top left corners of the cavity, respectively. The number following these abbreviations refer to the vortices that appear in the flow, which are numbered according to size. This simple geometry makes it possible to compare different numerical methods in detail. Most of the available results for this problem are based on a 129 × 129 uniform grid for moderate Reynolds numbers (Re ≤ 5000) and a 257 × 257 uniform grid for larger Reynolds numbers (Re ≥ 7500). In our computations, three different grids are employed: 65 × 65 grid mesh (0 ≤ Re ≤ 3200), 129 × 129 grid mesh (3200 ≤ Re ≤ 10 000), 257 × 257 grid mesh (1000 ≤ Re ≤ 17 500) and 513 × 513

P.X. Yu et al. / Computers and Mathematics with Applications (

)

–

13

Fig. 6. Comparisons of the u-velocity along the vertical centerline and the v -velocity along the horizontal centerline at (a) Re = 5000, (b) Re = 7500, (c) Re = 10 000, (d) Re = 12 500 (e) Re = 15 000 and (f) Re = 17 500 with 257 × 257 grids.

grid mesh (5000 ≤ Re ≤ 17 500). Numerical solutions at Reynolds numbers ranging from 0 to 17 500 are obtained using the point-successive iterative procedure directly. The inner–outer iteration technique is used to solve ψ and ω. In the inner iteration, we obtain the approximate stream function and vorticity by iteratively solving Eqs. (22) and (23), respectively. The number of iterations used for this purpose is fixed at 1,3, or 5 where a choice of 1 is equivalent to a traditional iteration (no inner iteration). After that, the approximate velocities can also be computed. The whole procedures constitute a complete outer iteration. The point-successive under-relaxation iteration procedures are used for low and medium-to-high Reynolds numbers, respectively. The initial conditions for the point-successive iterative procedure are calculated from the direct solution of Stokes flow (Re = 0). Then the flows at Re = 100, 400, 1000, 2000, 2500, 3200, 5000, 7500, 10 000, 12 500, 15 000 and 17 500 are solved using the previous solution as an initial condition. The iterative procedure is repeated until the

14

P.X. Yu et al. / Computers and Mathematics with Applications (

)

–

Fig. 7. Contours of the stream function (left) and the vorticity (right) for Re = 0, 100 and 400 with a grid of 65 × 65.

Table 7 Strength and location of the centers of primary vortex for the case with Re = 1000. Reference R4OC R4OC R4OC R4OC R4OC Botella [36]a Bruneau [43] Li [16] Zhang [23] Nishida [35] Erturk [3] Wahba [41] Ghia [37] a

Grid size 129 × 129 161 × 161 193 × 193 225 × 225 257 × 257 161 × 161 1024 × 1024 129 × 129 129 × 129 129 × 129 601 × 601 501 × 501 129 × 129

Accuracy 4

h h4 h4 h4 h4 h2 h4 h4 h8 h4 h4 h2

Note: The authors used highly accurate pseudo-spectral approach.

ψmin

ω

(x, y)

−0.11885 −0.11890 −0.11892 −0.11893 −0.11893 −0.11894 −0.11892 −0.11845 −0.11881 −0.11900 −0.11894 −0.11894 −0.11793

−2.0674 −2.0678 −2.0678 −2.0679 −2.0679 −2.0678 −2.0674 −2.0588 −2.0668 −2.0686 −2.0678 −2.0677 −2.0497

(0.5313, 0.5625) (0.5313, 0.5625) (0.5313, 0.5677) (0.5313, 0.5670) (0.5313, 0.5664) (0.5308, 0.5652) (0.5313, 0.5654) – (0.5313, 0.5625) (0.5313, 0.5625) (0.5300, 0.5650) (0.5300, 0.5660) (0.5313, 0.5625)

P.X. Yu et al. / Computers and Mathematics with Applications (

)

–

15

Fig. 8. Contours of the stream function (left) and the vorticity (right) for Re = 1000, 2500 and 3200 with a grid of 129 × 129.

maximum difference between successive approximations of both ψ and ω is smaller than 10−6 . To calculate the boundary value of vorticity for the moving boundary condition, we use h 21

(6ωw + 4ω1 − ω2 ) + O(h4 ) =

1 14h

(15ψw − 16ψ1 + ψ2 ) ± Vw

(37)

given by Spotz [42], where Vw is the tangential wall velocity, 1 and 2 denote the first two neighboring internal points on normal through the boundary w . For the driven cavity problem, Vw = 0 except on the moving lid, where Vw = u = 1. The convergence history of the residual errors on a 65 × 65 grid is shown in Fig. 2. In those calculations, the number of the inner iteration is set to be 1 which is equivalent to the traditional iteration and the convergence results can be obtained up to Re = 3200. In most cases, the iteration counts are lower than 5000. In the cases of the larger Reynolds numbers, although the residual errors would exhibit oscillating phenomenon, they would also decrease to the tolerances fast. Figs. 3 and 4 present comparisons of the horizontal velocities on the vertical centerline and the vertical velocities on the horizontal centerline of the square cavity for Reynolds numbers ranging from 100 to 10 000 and compare our data with that from Ghia et al. [37]. While Ghia’s data [37] was obtained using a 129 × 129 grid (Re = 100, 400, 1000 and 3200) and a 257 × 257 grid (Re = 3200, 5000, 7500 and 10 000), our data is obtained using a 65 × 65 grid (Re = 100, 400, 1000 and 3200), a 129 × 129 grid (Re = 3200, 5000, 7500 and 10 000). In each case, our profiles are in good agreement with Ghia’s results [37]. For a future assessment of the accuracy of the present results at different Reynolds number (from Re = 2500 to 17 500), the

16

P.X. Yu et al. / Computers and Mathematics with Applications (

)

–

Fig. 9. Contours of the stream function (left) and the vorticity (right) for Re = 5000, 7500 and 10 000 with a grid of 129 × 129.

velocity components through the vertical and horizontal centerlines of the cavity are also compared with the corresponding numerical results of Erturk et al. [25], where a fine grid mesh of 601 × 601 was used, in Figs. 5 and 6. The comparison shows excellent agreement. Remark. It should be emphasized that when a transient formulation is used for Reynolds numbers above 8050, transient behavior would certainly be observed. In fact, this conclusion has been studied by some researchers. For example, Bruneau and Saad [43] used a 1024 × 1024 mesh, and determined the critical number to be between 8000 and 8050. In the present work, however, steady solutions can be obtained for Reynolds numbers above 8050, because the steady formulation is used, which is equivalent to assuming at the outset that the solution is steady. The computed streamlines and vorticity contours for 0 ≤ Re ≤ 17 500 are shown in Figs. 7–10. Notice that the streamlines and vorticity contours for Re = 0 are symmetric about the vertical centerline of the cavity and the BL and BR secondary eddies are equal in size. Figs. 7–10 also show that as the Reynolds number increases the vorticity contours move away from the cavity center towards the cavity walls and indicate that very strong vorticity gradients develop on the lid and the cavity walls at higher Reynolds numbers. The streamlines in Figs. 8–10 show as the Reynolds number increases the primary vortex moves gradually towards the cavity center. These figures depict the typical separations and secondary vortices BL and BR at the bottom corners of the cavity as well as TL at the top left of the square cavity. From Figs. 9 and 10, the evolution of tertiary vortex in the right bottom corners of the cavity is also observed. The tertiary vortex becomes visible at

P.X. Yu et al. / Computers and Mathematics with Applications (

)

–

17

Fig. 10. Contours of the stream function (left) and the vorticity (right) for Re = 12 500, 15 000 and 17 500 with a grid of 257 × 257.

Fig. 11. Schematic view of the backward-facing step flow problem.

Re = 5000 and gains a significant size for Re = 7500. If the Reynolds number is increased further, the evolution of another tertiary vortex in the left bottom corners of the cavity can be found at Re = 7500 and becomes significant for Re = 10 000.

18

P.X. Yu et al. / Computers and Mathematics with Applications (

)

–

Fig. 12. The streamlines for the backward-facing step flow problem in: (a) Re = 400, (b) Re = 500, (c) Re = 600, and (d) Re = 800.

Further increases in the Reynolds number make visible tertiary vortex in the top left of the cavity. Fig. 10 shows that the tertiary vortex in the top left of the cavity appears at Re = 12 500 and the size becomes larger for Re = 15 000. The quantitative comparison data for our solutions are tabulated in Tables 7–12. First, in Table 7, we compare the values of stream function, vorticity and location of primary vortices with different grid sizes. From the table, it is found that the value of stream function at the primary vortex obtained by the present method is −0.11892 with a coarser grid size of 193 × 193, and is in excellent agreement with the recent ‘‘benchmark’’ results which are around −0.11892 to −0.11894 [43,36,3,41] using fine grids or pseudo-spectral approach. On the other hand, even using the 129 × 129 grid, our result only exists less than 0.1% error with the ‘‘benchmark’’ solution, which is more accurate than other high-order methods such as Zhang [23], Li [16] and Nishida [35]. All these reflect that the present method has the feature of high accuracy. For future reference, we list the results for other higher Reynolds number with the grid sizes of 257 × 257 and 513 × 513. For comparison, we also exhibit some results with coarser grid. In Table 8, the stream function and vorticity values at the center of the primary vortices and those locations are listed for 0 ≤ Re ≤ 17 500. The present results are in excellent agreement with the best and most accurate solutions [36,3,25] available in the literature. In Table 9, the data on the secondary vortices in the bottom left and right corners of the square cavity are given for 0 ≤ Re ≤ 17 500 and are compared with other results in the literature. The results are in excellent agreement. For 3200 ≤ Re ≤ 17 500, the data on the secondary vortex at the top left side of the cavity are tabulated in Table 10. Our results also exhibit an excellent match with the published results. In Table 11, the data on the tertiary vortices at the bottom left and right corners for 5000 ≤ Re ≤ 17 500 along with similar results found in the literature are presented. Our computed results again agreed well with the published results. For 12 500 ≤ Re ≤ 17 500, our solutions on the tertiary vortex at the top left side of the cavity are listed in Table 12. Our results also match very well with those in the literature. Table 13 tabulates the stream function and vorticity values at the center of the primary, secondary and tertiary vortices and the corresponding location of the center of these vortices. The results in this table are in good agreement with that of Erturk et al. [3,25].

4.6. Problem 6 Finally, we consider the backward-facing step flow problem in the region of [0, L] × [−h, h], where h is the height of the step and the width of inflow and L is the channel length (see Fig. 11). This flow problem has been solved numerically by many scholars to validate their methods [45,46]. For this flow problem, the parabolic velocity profile is employed to simulate the

P.X. Yu et al. / Computers and Mathematics with Applications (

)

–

19

Table 8 Strength and location of the centers of primary vortex for the case with Re = 0, 3200, 5000, 7500, 10 000, 12 500, 15 000 and 17 500. Reference

Grid size

ψmin

ω

(x, y)

R4OC R4OC Sahin [40] Botella [36]

65 × 65 129 × 129 257 × 257 49 × 49

−0.10008 −0.10008 −0.10005 −0.10008

−3.2202 −3.2205 −3.2208 −

(0.5000, 0.7656) (0.5000, 0.7656) (0.5000, 0.7626) –

3 200

R4OC R4OC Ghia [37] Li [16] Gupta [44]

129 × 129 257 × 257 129 × 129 129 × 129 161 × 161

−0.12110 −0.12173 −0.12038 −0.12053 −0.122

−1.9531 −1.9604 −1.9886 −1.9429

(0.5156, 0.5391) (0.5195, 0.5391) (0.5165, 0.5469) – (0.5188, 0.5488)

5 000

R4OC R4OC R4OC Ghia [37] Li [16] Gupta [44] Erturk [3]

129 × 129 257 × 257 513 × 513 257 × 257 129 × 129 161 × 161 601 × 601

−0.12090 −0.12201 −0.12220 −0.11897 −0.12036 −0.122 −0.12222

−1.9243 −1.9383 −1.9406 −1.8602 −1.9243 −1.9406

(0.5156, 0.5391) (0.5156, 0.5352) (0.5156, 0.5352) (0.5117, 0.5352) – (0.5125, 0.5375) (0.5150, 0.5350)

7 500

R4OC R4OC R4OC Erturk [3]

129 × 129 257 × 257 513 × 513 601 × 601

−0.12036 −0.12197 −0.12233 −0.12234

−1.9003 −1.9219 −1.9268 −1.9265

(0.5156, 0.5313) (0.5117, 0.5313) (0.5137, 0.5313) (0.5133, 0.5317)

10 000

R4OC R4OC R4OC Erturk [3]

129 × 129 257 × 257 513 × 513 601 × 601

−0.12041 −0.12174 −0.12231 −0.12231

−1.9009 −1.9102 −1.9188 −1.9182

(0.5156, 0.5313) (0.5117, 0.5313) (0.5117, 0.5293) (0.5117, 0.5300)

12 500

R4OC R4OC Erturk [3]

257 × 257 513 × 513 601 × 601

−0.12149 −0.12223 −0.12220

−1.9017 −1.9132 −1.9123

(0.5117, 0.5273) (0.5117, 0.5293) (0.5117, 0.5283)

15 000

R4OC R4OC Erturk [3]

257 × 257 513 × 513 601 × 601

−0.12151 −0.12214 −0.12206

−1.9021 −1.9089 −1.9077

(0.5117, 0.5273) (0.5098, 0.5273) (0.5100, 0.5283)

17 500

R4OC R4OC Erturk [3]

257 × 257 513 × 513 601 × 601

−0.12098 −0.12201 −0.12189

−1.8885 −1.9035 −1.9037

(0.5117, 0.5273) (0.5098, 0.5273) (0.5100, 0.5283)

Re 0

–

–

incoming flow, and the function for 0 ≤ y ≤ 0.5 is given by [46] as u = 12y − 24y2 ,

(38)

v = 0,

(39)

ψ = 6y − 8y . 2

3

(40)

For the outlet (the right side), we apply Neumann boundary for stream function ψ and vorticity ω because the total length of the tunnel L is assumed to be long enough. All other boundaries used the wall boundaries. Here, we set ψ = 0 on the left and bottom wall and ψ = 0.5 on the upper wall considering the continuity of ψ , and the same scheme (37) is employed to calculate the boundary value of the vorticity on the walls. Now, we deduce the boundary scheme for ω at the inlet. Referring the wall boundary scheme (37), we can easily obtain the following equation using the Taylor series expansion for Vw = 0. 1

(15ψw − 16ψ1 + ψ2 ) = −6hψxx −

4h2

ψxxx + O(h4 ).

(41)

ψxx = −ω − ψyy = −ω + (−12 + 48y) , −ω + g (y), ψxxx = −ωx − ψxyy = −ωx + vyy = −ωx .

(42)

h For the inlet boundary, we notice

3

(43)

Using the second order scheme for ωx , we have 1 (3ωw − 4ω1 + ω2 ) + O(h2 ). 2h Substituting Eqs. (42) and (44) to Eq. (41), the boundary scheme for vorticity at the inlet is

ψxxx =

1 h

(15ψw − 16ψ1 + ψ2 ) =

2h 3

(6ωw − 9g (y) + 4ω1 − ω2 ) + O(h4 ).

(44)

(45)

20

P.X. Yu et al. / Computers and Mathematics with Applications (

)

–

Table 9 Strength and location of the centers of the secondary vortex BL-2 and BR-2 for the case with Re = 0, 1000, 3200, 5000, 7500, 10 000, 12 500, 15 000 and 17 500. Re

Reference

Grid size

BL-2

BR-2

ψmax

(x, y)

ψmax

(x, y)

R4OC R4OC Sahin [40]

65 × 65 129 × 129 257 × 257

2.20e−6 2.20e−6 2.22e−6

(0.0313, 0.0469) (0.0391, 0.0391) (0.0369, 0.0378)

2.20e−6 2.20e−6 2.22e−6

(0.9687, 0.0469) (0.9609, 0.0391) (0.9630, 0.0378)

1 000

R4OC R4OC Ghia [37] Botella [36] Gupta [44] Erturk [3]

129 × 129 257 × 257 129 × 129 161 × 161 81 × 81 601 × 601

2.31e−4 2.33e−4 2.31e−4 2.33e−4 2.02e−4 2.33e−4

(0.0859, 0.0781) (0.0820, 0.0781) (0.0859, 0.0781) (0.0833, 0.0781) (0.0875, 0.0750) (0.0833, 0.0783)

1.73e−3 1.73e−3 1.75e−3 1.73e−3 1.70e−3 1.73e−3

(0.8594, 0.1094) (0.8633, 0.1133) (0.8594, 0.1094) (0.8640, 0.1118) (0.8625, 0.1125) (0.8633, 0.1117)

3 200

R4OC R4OC Ghia [37] Gupta [44]

129 × 129 257 × 257 129 × 129 161 × 161

1.08e−3 1.11e−3 9.78e−4 1.03e−3

(0.0781, 0.1250) (0.0820, 0.1172) (0.0859, 0.1094) (0.0813, 0.1188)

2.81e−3 2.83e−3 3.14e−3 2.86e−3

(0.8281, 0.0859) (0.8242, 0.0859) (0.8125, 0.0859) (0.8125, 0.0875)

5 000

R4OC R4OC R4OC Ghia [37] Gupta [44] Erturk [3]

129 × 129 257 × 257 513 × 513 257 × 257 161 × 161 601 × 601

1.33e−3 1.37e−3 1.38e−3 1.36e−3 1.32e−3 1.38e−3

(0.0703, 0.1406) (0.0742, 0.1367) (0.0723, 0.1367) (0.0703, 0.1367) (0.0750, 0.1313) (0.0733, 0.1367)

3.02e−3 3.07e−3 3.07e−3 3.08e−3 2.96e−3 3.07e−3

(0.7969, 0.0703) (0.8047, 0.0742) (0.8047, 0.0723) (0.8086, 0.0742) (0.8000, 0.0750) (0.8050, 0.0733)

7 500

R4OC R4OC R4OC Erturk [3]

129 × 129 257 × 257 513 × 513 601 × 601

1.46e−3 1.52e−3 1.53e−3 1.53e−3

(0.0625, 0.1563) (0.0625, 0.1563) (0.0645, 0.1523) (0.0650, 0.1517)

3.11e−3 3.22e−3 3.23e−3 3.23e−3

(0.7891, 0.0625) (0.7930, 0.0664) (0.7891, 0.0645) (0.7900, 0.0650)

10 000

R4OC R4OC R4OC Erturk [3]

129 × 129 257 × 257 513 × 513 601 × 601

1.51e−3 1.60e−3 1.62e−3 1.61e−3

(0.0625, 0.1563) (0.0586, 0.1641) (0.0586, 0.1621) (0.0583, 0.1633)

3.05e−3 3.17e−3 3.19e−3 3.19e−3

(0.7734, 0.0625) (0.7734, 0.0586) (0.7734, 0.0586) (0.7750, 0.0600)

12 500

R4OC R4OC Erturk [3]

257 × 257 513 × 513 601 × 601

1.64e−3 1.66e−3 1.66e−3

(0.0547, 0.1680) (0.0547, 0.1680) (0.0550, 0.1683)

3.07e−3 3.10e−3 3.10e−3

(0.7617, 0.0547) (0.7598, 0.0547) (0.7600, 0.0550)

15 000

R4OC R4OC Erturk [3]

257 × 257 513 × 513 601 × 601

1.65e−3 1.68e−3 1.67e−3

(0.0547, 0.1719) (0.0526, 0.1719) (0.0533, 0.1717)

2.97e−3 3.00e−3 3.00e−3

(0.7461, 0.0508) (0.7480, 0.0508) (0.7450, 0.0500)

17 500

R4OC R4OC Erturk [3]

257 × 257 513 × 513 601 × 601

1.64e−3 1.66e−3 1.65e−3

(0.0508, 0.1758) (0.0508, 0.1777) (0.0500, 0.1767)

2.85e−3 2.90e−3 2.90e−3

(0.7344, 0.0469) (0.7344, 0.0469) (0.7333, 0.0467)

0

For the outlet boundary, the one-side four-order scheme is used for both ψ and ω.

ψw = ωw =

1 25 1 25

(48ψ1 − 36ψ2 + 16ψ3 − 3ψ4 ) + O(h4 )

(46)

(48ω1 − 36ω2 + 16ω3 − 3ω4 ) + O(h4 ).

(47)

Like the lid driven cavity case, the iterative procedure of this case is also repeated until the maximum difference between successive approximations of both ψ and ω is smaller than 10−6 . According to the results of the grid independence study in Table 14, it is found that the differences between the values of the stream function with 1201 × 41 and 1501 × 51 grids are lower than 0.5%. Then, we decide to choose the grid size of 1201 × 41 for other Reynolds numbers. Fig. 12 displays the streamlines of the backward-facing step flow problem ranging from Re = 400 to Re = 800. With the increase of Reynolds number, it can be seen that the eddy on the bottom wall becomes large. For the larger Re, another eddy would appear on the top wall. All these are identical to the flow patterns in the preceding publications. Quantitative results including the value and the location of the stream function in the top and bottom wall eddies center, the recirculation point x1 of the bottom wall eddy and the separation point x2 and recirculation point x3 of the top wall eddy of the backward-facing step flow problem are tabulated in Tables 15 and 16. Compared with the other available publications, our numerical results are in good agreement with those previous data again. Overall, with the proper numerical boundary for the vorticity, our present method can be employed to solve the different flow problems rather than the cavity flows.

P.X. Yu et al. / Computers and Mathematics with Applications (

)

–

21

Table 10 Strength and location of the centers of the secondary vortex TL-2 for the case with Re = 3200, 5000, 7500, 10 000, 12 500, 15 000 and 17 500. Reference

Grid size

ψmax

(x, y)

3 200

R4OC R4OC Ghia [37] Gupta [44]

129 × 129 257 × 257 129 × 129 161 × 161

6.80e−4 7.03e−4 7.28e−4 7.33e−4

(0.0547, 0.8984) (0.0547, 0.8984) (0.0547, 0.8984) (0.0563, 0.9000)

5 000

R4OC R4OC R4OC Ghia [37] Gupta [44] Erturk [3]

129 × 129 257 × 257 513 × 513 257 × 257 161 × 161 601 × 601

1.38e−3 1.43e−3 1.44e−3 1.46e−3 1.54e−3 1.45e−3

(0.0625, 0.9063) (0.0625, 0.9102) (0.0625, 0.9102) (0.0625, 0.9102) (0.0688, 0.9125) (0.0633, 0.9100)

7 500

R4OC R4OC R4OC Erturk [3]

129 × 129 257 × 257 513 × 513 601 × 601

1.99e−3 2.11e−3 2.13e−3 2.13e−3

(0.0625, 0.9141) (0.0664, 0.9102) (0.0664, 0.9121) (0.0667, 0.9117)

10 000

R4OC R4OC R4OC Erturk [3]

129 × 129 257 × 257 513 × 513 601 × 601

2.36e−3 2.59e−3 2.62e−3 2.62e−3

(0.0703, 0.9141) (0.0703, 0.9102) (0.0703, 0.9102) (0.0700, 0.9100)

12 500

R4OC R4OC Erturk [3]

257 × 257 513 × 513 601 × 601

2.94e−3 3.00e−3 2.99e−3

(0.0742, 0.9102) (0.0742, 0.9102) (0.0733, 0.9100)

15 000

R4OC R4OC Erturk [3]

257 × 257 513 × 513 601 × 601

3.17e−3 3.29e−3 3.29e−3

(0.0742, 0.9102) (0.0762, 0.9102) (0.0767, 0.9100)

17 500

R4OC R4OC Erturk [3]

257 × 257 513 × 513 601 × 601

3.44e−3 3.53e−3 3.52e−3

(0.0781, 0.9102) (0.0781, 0.9121) (0.0783, 0.9117)

Re

Table 11 Strength and location of the centers of the tertiary vortex BL-3 and BR-3 for the case with Re = 5000, 7500, 10 000, 12 500, 15 000 and 17 500. Re

Reference

Grid size

BL-3

BR-3

ψmin

(x, y)

ψmin

(x, y)

5 000

R4OC R4OC Ghia [37] Gupta [44] Erturk [3]

257 × 257 513 × 513 257 × 257 161 × 161 601 × 601

−6.61e−8 −6.65e−8 −7.09e−8 −5.15e−8 −6.55e−8

(0.0078, 0.0078) (0.0078, 0.0078) (0.0117, 0.0078) (0.0063, 0.0063) (0.0083, 0.0083)

−1.40e−6 −1.42e−6 −1.43e−6 −1.70e−6 −1.43e−6

(0.9805, 0.0195) (0.9785, 0.0195) (0.9805, 0.0195) (0.9750, 0.0188) (0.9783, 0.0183)

7 500

R4OC R4OC Ghia [37] Gupta [44] Erturk [3]

257 × 257 513 × 513 257 × 257 161 × 161 601 × 601

−1.93e−7 −2.02e−7 −1.83e−7 −1.64e−7 −2.02e−7

(0.0117, 0.0117) (0.0117, 0.0117) (0.0117, 0.0117) (0.0063, 0.0125) (0.0117, 0.0117)

−3.20e−5 −3.27e−5 −3.28e−5 −1.89e−5 −3.27e−5

(0.9531, 0.0430) (0.9512, 0.0410) (0.9492, 0.0430) (0.9500, 0.0375) (0.9517, 0.0417)

10 000

R4OC R4OC Ghia [37] Erturk [3]

257 × 257 513 × 513 257 × 257 601 × 601

−1.01e−6 −1.11e−6 −7.76e−7 −1.10e−6

(0.0156, 0.0195) (0.0176, 0.0195) (0.0156, 0.0195) (0.0167, 0.0200)

−1.36e−4 −1.40e−4 −1.31e−4 −1.40e−4

(0.9336, 0.0664) (0.9355, 0.0684) (0.9336, 0.0625) (0.9350, 0.0683)

12 500

R4OC R4OC Erturk [3]

257 × 257 513 × 513 601 × 601

−5.95e−6 −6.62e−6 −6.45e−6

(0.0273, 0.0313) (0.0273, 0.0313) (0.0267, 0.0317)

−2.46e−4 −2.55e−4 −2.55e−4

(0.9297, 0.0820) (0.9277, 0.0820) (0.9283, 0.0817)

15 000

R4OC R4OC Erturk [3]

257 × 257 513 × 513 601 × 601

−2.05e−5 −2.33e−5 −2.23e−5

(0.0352, 0.0430) (0.0391, 0.0410) (0.0383, 0.0417)

−3.34e−4 −3.40e−4 −3.40e−4

(0.9258, 0.0898) (0.9258, 0.0879) (0.9267, 0.0883)

17 500

R4OC R4OC Erturk [3]

257 × 257 513 × 513 601 × 601

−4.70e−5 −5.10e−5 −4.83e−5

(0.0508, 0.0469) (0.0508, 0.0488) (0.0500, 0.0483)

−3.86e−4 −4.04e−4 −4.05e−4

(0.9258, 0.0938) (0.9277, 0.0957) (0.9283, 0.0967)

5. Conclusions In this paper, a new HOC, referred to as RHOC, finite difference method has been formulated for solving the twodimensional steady stream function–vorticity formulation of the Navier–Stokes equations. As a basis of the discretization of the 2D incompressible steady stream function–vorticity form of the Navier–Stokes equations, the R4OC scheme for the

22

P.X. Yu et al. / Computers and Mathematics with Applications (

)

–

Table 12 Strength and location of the centers of the tertiary vortex TL-3 for the case with Re = 12 500, 15 000 and 17 500. Re

Reference

Grid size

ψmin

(x, y)

12 500

R4OC R4OC Erturk [3]

257 × 257 513 × 513 601 × 601

−1.09e−6 −1.56e−6 −1.63e−6

(0.0078, 0.8320) (0.0078, 0.8301) (0.0067, 0.8300)

15 000

R4OC R4OC Erturk [3]

257 × 257 513 × 513 601 × 601

−1.23e−5 −1.58e−5 −1.58e−5

(0.0156, 0.8281) (0.0156, 0.8262) (0.0150, 0.8250)

17 500

R4OC R4OC Erturk [3]

257 × 257 513 × 513 601 × 601

−3.74e−5 −4.15e−5 −4.12e−5

(0.0195, 0.8242) (0.0215, 0.8223) (0.0200, 0.8217)

Table 13 Properties of the strength and the location of the vortex. Re

3200

5000

7500

10 000

12 500

15 000

17 500

ψ ω (x , y )

−0.12178 −1.96067

−0.12220 −1.94060

−0.12233 −1.92677

−0.12231 −1.91875

−0.12223 −1.91325

−0.12214 −1.90894

−0.12201 −1.90353

(0.5175,0.5391)

(0.5156,0.5352)

(0.5137,0.5313)

(0.5117,0.5293)

(0.5117,0.5293)

(0.5098,0.5273)

(0.5098,0.5273)

BR-2

ψ ω (x , y )

2.8285e−3 2.1974 (0.8242,0.0840)

3.0726e−3 2.7514 (0.8047,0.0723)

3.2256e−3 3.2682 (0.7891,0.0645)

3.1892e−3 3.8015 (0.7734,0.0586)

3.0972e−3 4.3628 (0.7598,0.0547)

3.0009e−3 4.9812 (0.7480,0.0508)

2.9020e−3 5.5672 (0.7344,0.0469)

BR-3

ψ ω (x, y)

−1.9339e−7 −1.4063e−2

−1.4205e−6 −3.6865e−2

−3.2670e−5 −0.1546

−1.4007e−4 −0.3069

−2.5496e−4 −0.4067

−3.4011e−4 −0.4595

−4.0499e−4 −0.5113

(0.9883,0.0098)

(0.9785,0.0195)

(0.9512,0.0410)

(0.9355,0.0684)

(0.9277,0.0820)

(0.9258,0.0879)

(0.9277,0.0957)

BL-2

ψ ω (x, y)

1.1157e−3 1.1620 (0.0820,0.1191)

1.3759e−3 1.5009 (0.0723,0.1367)

1.5347e−3 1.8724 (0.0645,0.1523)

1.6162e−3 2.1531 (0.0586,0.1621)

1.6627e−3 2.3632 (0.0547,0.1680)

1.6772e−3 2.5206 (0.0527,0.1719)

1.6612e−3 2.7113 (0.0508,0.1777)

BL-3

ψ ω (x, y)

−3.4759e−8 −9.2004e−3

−6.6546e−8 −9.9701e−3

−2.0161e−7 −1.6773e−2

−1.1143e−6 −3.1569e−2

−6.6243e−6 −7.5045e−2

−2.3268e−5 −0.1393

−5.1067e−5 −0.2056

(0.0059,0.0078)

(0.0078,0.0078)

(0.0117,0.0117)

(0.0176,0.0195)

(0.0273,0.0313)

(0.0391,0.0410)

(0.0508,0.0488)

TL-2

ψ ω (x, y)

7.0875e−4 1.6902 (0.0527,0.8984)

1.4445e−3 2.1240 (0.0625,0.9102)

2.1292e−3 2.2431 (0.0664,0.9121)

2.6226e−3 2.3184 (0.0703,0.9102)

2.9950e−3 2.3772 (0.0742,0.9102)

3.2887e−3 2.4272 (0.0762,0.9102)

3.5327e−3 2.4721 (0.0781,0.9121)

TL-3

ψ ω (x , y )

– – –

– – –

– – –

– – –

−1.5662e−6 −0.2702

−1.5844e−5 −0.5692

−4.1512e−5 −0.8243

(0.0078,0.8301)

(0.0156,0.8262)

(0.0215,0.8223)

Primary

Table 14 Results of the grid independence study for the backward-facing step flow problem with Re = 800. Grid

601 × 21 901 × 31 1201 × 41 1501 × 51

Bottom vortex

Top vortex

ψ

x

y

ψ

x

y

−0.03520 −0.03461 −0.03442 −0.03434

3.300 3.333 3.325 3.340

−0.200 −0.200 −0.200 −0.200

0.5072 0.5068 0.5067 0.5066

7.250 7.333 7.475 7.460

0.300 0.300 0.325 0.320

model convection–diffusion equation was first established. And, it has also been found that the present scheme is actually the standard fourth-order Padé scheme for the pure convective problems. The algebraic system of equations associated with the currently proposed finite difference approximations can be solved by using the point-successive iterative procedure directly. To demonstrate that fourth-order convergence, four problems which include two linear and two nonlinear problems with their exact solutions are first solved by using the present method. The newly proposed method is then applied to the lid driven square cavity and the backward-facing step flow problems. The computational results show that the present method has the advantage of better scale resolution with smaller number of grid nodes.

P.X. Yu et al. / Computers and Mathematics with Applications (

)

–

23

Table 15 Comparison of the results of bottom vortex for the backward-facing step flow problem. Re

Reference

ψ

x

y

x1

400

Present Yu and Tian [46] Kalita and Gupta [45]

−0.03430 −0.03330 −0.03326

1.800 1.672 1.750

−0.200 −0.181 −0.200

4.275 4.394 4.200

500

Present Yu and Tian [46] Kalita and Gupta [45]

−0.03439 −0.03332 −0.03328

2.200 2.377 2.150

−0.200 −0.222 −0.200

4.850 5.046 4.750

600

Present Yu and Tian [46] Kalita and Gupta [45]

−0.03442 −0.03337 −0.03370

2.575 2.760 2.438

−0.200 −0.222 −0.200

5.325 5.591 5.313

800

Present Yu and Tian [46] Kalita and Gupta [45] Gartling [47] Keskar and Lyn [48]

−0.03442 −0.03342 −0.03381 −0.03420 −0.03420

3.325 3.558 3.525 3.350 –

−0.200 −0.222 −0.218 −0.200

6.025 6.190 6.031 6.100 6.097

–

Table 16 Comparison of the results of top vortex for the backward-facing step flow problem. Re

Reference

ψ

x

y

x2

x3

500

Present Yu and Tian [46] Kalita and Gupta [45]

0.5007 0.5004 0.5006

5.450 5.404 5.400

0.425 0.431 0.400

4.275 4.244 4.163

6.675 6.621 6.613

600

Present Yu and Tian [46] Kalita and Gupta [45]

0.5025 0.5017 0.5023

6.200 6.190 6.175

0.375 0.386 0.375

4.475 4.550 4.468

8.025 8.090 7.969

800

Present Yu and Tian [46] Kalita and Gupta [45] Gartling [47] Keskar and Lyn [48]

0.5067 0.5058 0.5065 0.5065 0.5065

7.475 7.570 7.375 7.400 –

0.325 0.330 0.313 0.300 –

4.950 5.046 4.876 4.850 4.850

10.375 10.530 10.312 10.480 10.480

Acknowledgments This work was supported in part by the National Natural Science Foundation of China under Grant 11372075, 91330112 and 11502054, China Postdoctoral Science Foundation under Grant 2014M550211 and Special Program for Applied Research on Super Computation of the NSFC-Guangdong Joint Fund (the second phase). References [1] D.W. Zingg, Comparison of high-accuracy finite-difference methods for linear wave propagation, SIAM J. Sci. Comput. 22 (2000) 476–502. [2] J.Y. Choo, D.H. Schultz, A stable high order method for the heated cavity problem, Internat. J. Numer. Methods Fluids 15 (1992) 1313–1332. [3] E. Erturk, C. Gökçöl, Fourth-order compact formulation of Navier–Stokes equations and driven cavity flow at high Reynolds numbers, Internat. J. Numer. Methods Fluids 50 (2006) 421–436. [4] J.C. Strikwerda, High-order accurate schemes for incompressible viscous flow, Internat. J. Numer. Methods Fluids 24 (1997) 715–734. [5] A. Mohebbia, M. Dehghan, High-order compact solution of the one-dimensional heat and advection-diffusion equations, Appl. Math. Model. 34 (2010) 3071–3084. [6] Y. Ge, Multigrid method and fourth-order compact difference discretization scheme with unequal meshsizes for 3D Poisson equation, J. Comput. Phys. 229 (2010) 6381–6391. [7] Y. Ge, F. Cao, Multigrid method based on the transformation-free HOC scheme on nonuniform grids for 2D convection diffusion problems, J. Comput. Phys. 230 (2011) 4051–4070. [8] Y. He, J. Li, Convergence of three iterative methods based on the finite element discretization for the stationary Navier–Stokes equations, Comput. Methods Appl. Mech. Engrg. 198 (2009) 1351–1359. [9] Y. He, L. Mei, Y. Shang, J. Cui, Newton iterative parallel finite element algorithm for the steady Navier–Stokes equations, J. Sci. Comput. 44 (2010) 92–106. [10] H. Xu, Y. He, Some iterative finite element methods for steady Navier–Stokes equations with different viscosities, J. Comput. Phys. 232 (2013) 136–152. [11] Y. He, Y. Zhang, Y. Shang, H. Xu, Two-level Newton iterative method for the 2D/3D steady Navier–Stokes equations, Numer. Methods Partial Differential Equations 28 (2012) 1620–1642. [12] R. Steijl, H.W.M. Hoeijmakers, Fourth-order accurate compact-difference discretization method for Helmholtz and incompressible Navier–Stokes equations, Internat. J. Numer. Methods Fluids 46 (2004) 227–244. [13] I. Singer, E. Turkel, High-order finite difference methods for the Helmholtz equation, Comput. Methods Appl. Mech. Engrg. 163 (1998) 343–358. [14] M.M. Gupta, R.P. Manohar, J.W. Stephenson, A single cell high order scheme for the convection–diffusion equation with variable coefficients, Internat. J. Numer. Methods Fluids 4 (1984) 641–651. [15] W.F. Spotz, G.F. Carey, High-order compact scheme for the steady stream-function vorticity equations, Internat. J. Numer. Methods Engrg. 38 (1995) 3497–3512. [16] M. Li, T. Tang, B. Fornberg, A compact fourth-order finite difference scheme for the incompressible Navier–Stokes equations, Internat. J. Numer. Methods Fluids 20 (1995) 1137–1151.

24

P.X. Yu et al. / Computers and Mathematics with Applications (

)

–

[17] K.K.Q. Zhang, B. Shotorban, W.J. Minkowycz, F. Mashayek, A compact finite difference method on staggered grid for Navier–Stokes flows, Internat. J. Numer. Methods Fluids 52 (2006) 867–881. [18] M.M. Gupta, High-order solution on incompressible Navier–Stokes equations, J. Comput. Phys. 93 (1991) 343–359. [19] Z.F. Tian, S.Q. Dai, High-order compact exponential finite difference methods for convection–diffusion type problems, J. Comput. Phys. 220 (2007) 952–974. [20] Z.F. Tian, A rational high-order compact ADI method for unsteady convection–diffusion equations, Comput. Phys. Comm. 182 (2011) 649–662. [21] S.C.R. Dennis, J.D. Hudson, Compact h4 finite-difference approximations to operators of Navier–Stokes type, J. Comput. Phys. 85 (1989) 390–416. [22] R.J. MacKinnon, R.W. Johnson, Differential equation based representtaion of truncation errors for accurate numerical simulation, Internat. J. Numer. Methods Fluids 13 (1991) 739–757. [23] J. Zhang, Numerical simulation of 2D square driven cavity using fourth-order compact finite difference schemes, Comput. Math. Appl. 45 (2003) 43–52. [24] J.C. Kalita, D.C. Dalal, A.K. Dass, Fully compact higher-order computation of steady-state natural convection in a square cavity, Phys. Rev. E 64 (2001) 066703(1–13). [25] E. Erturk, T.C. Corke, C. Gökçöl, Numerical solutions of 2-D steady incompressible driven cavity flow at high Reynolds numbers, Internat. J. Numer. Methods Fluids 48 (2005) 747–774. [26] J. Kouatchou, Asymptotic stability of a nine-point multigrid algorithm for convection diffusion equations, Electron. Trans. Numer. Anal. 6 (1997) 153–161. [27] M.M. Gupta, J. Kouatchou, J. Zhang, A compact multigrid solver for convection diffusion equations, J. Comput. Phys. 132 (1997) 123–129. [28] J. Zhang, Preconditioned iterative methods and finite difference schemes for convection–diffusion, Appl. Math. Comput. 109 (2000) 11–30. [29] D.F. Roscoe, New methods for the derivation of stable difference representations for differential equations, J. Inst. Math. Appl. 16 (1975) 291–301. [30] S.C.R. Dennis, Finite difference associated with second-order difference equations, Quart. J. Mech. Appl. Math. 13 (1960) 487–507. [31] S.C.R. Dennis, J.D. Hudson, Accurate representations of partial differential equations by finite difference schemes, J. Inst. Math. Appl. 23 (1979) 43–51. [32] A.C.R. Pillai, Fourth-order exponential finite difference methods for boundary value problems of convective diffusion type, Internat. J. Numer. Methods Fluids 37 (2001) 87–106. [33] G.Q. Chen, Z. Gao, Z.F. Yang, A perturbational h4 exponential finite difference scheme for convection diffusion equation, J. Comput. Phys. 104 (1993) 129–139. [34] T.M. Shih, C.H. Tan, B.C. Hwang, Effects of grid staggering on numerical schemes, Internat. J. Numer. Methods Fluids 9 (1989) 193–212. [35] H. Nishida, N. Satofuka, Higher-order solutions of square driven cavity flow using a variable-order multigrid method, Internat. J. Numer. Methods Engrg. 34 (1992) 637–653. [36] O. Botella, R. Peyret, Benchmark spectral results on the lid-driven cavity flow, Comput. & Fluids 27 (1998) 421–433. [37] U. Ghia, K.N. Ghia, C.T. Shin, High Re-solution for incompressible Navier–Stokes equation and a multigrid method, J. Comput. Phys. 48 (1982) 387–411. [38] E. Erturk, Discussions on driven cavity flow, Internat. J. Numer. Methods Fluids 60 (2009) 275–294. [39] S. Hou, Q. Zou, S. Chen, G.A.Cogley Doolen, Simulation of cavity flows by the lattice Boltzmann method, J. Comput. Phys. 118 (1995) 329–347. [40] M. Sahin, R.G. Owens, A novel fully-implicit finite volume method applied to the lid-driven cavity problem. Part I. High Reynolds number flow calculations, Internat. J. Numer. Methods Fluids 42 (2003) 57–77. [41] E.M. Wahba, Steady flow simulations inside a driven cavity up to Reynolds number 35000, Comput. & Fluids 66 (2012) 85–97. [42] W.F. Spotz, Accuracy and performance of numerical wall boundary conditions for steady, 2D, incompressible streamfunction vorticity, Internat. J. Numer. Methods Fluids 28 (1998) 737–757. [43] C.H. Bruneau, M. Saad, The 2D lid-driven cavity problem revisited, Comput. & Fluids 35 (2006) 326–348. [44] M.M. Gupta, J.C. Kalita, A new paradigm for solving Navier–Stokes equations: streamfunction-velocity formulation, J. Comput. Phys. 207 (2005) 52–68. [45] J.C. Kalita, M.M. Gupta, A streamfunction-velocity approach for 2D transient incompressible viscous flows, Internat. J. Numer. Methods Fluids 62 (2010) 237–266. [46] P.X. Yu, Z.F. Tian, A compact streamfunction-velocity scheme on nonuniform grids for the 2D steady incompressible Navier–Stokes equations, Comput. Math. Appl. 66 (2013) 1192–1212. [47] D.K. Gartling, A test problem for outflow boundary conditions-flow over a backward facing step, Internat. J. Numer. Methods Fluids 11 (1990) 953–967. [48] J. Keskar, D.A. Lyn, Computations of a laminar backward-facing step flow at Re = 800 with a spectral domain decomposition method, Internat. J. Numer. Methods Fluids 29 (1999) 411–427.