Chinese Journal of Numerical Mathematics and Applications Vol. 23, No. 4, Nov., 2001, 447–456

Iterative Pressure Poisson Equation Method for Solving the Unsteady Incompressible N–S Equations

YU Xin ^{1}

(LHD, Institute of Mechanics, Chinese Academy of Sciences, Beijing)

Abstract ^{2} ^{3}

We present an iterative pressure Poisson equation method for solving the viscous incompressible Navier–Stokes equations. Comparing with the pressure Poisson equation method, it uses the increment of the pressure, the difference of two successive pressures (p_{k}^{n+1} and p_{k+1}^{n+1}), instead of the pressure p^{n+1} as an unknown variable. We apply it to a fourth order accurate staggered mesh compact difference scheme, simulate the driven flow in a square cavity with Re = 5000 and 10000.

Key Words: N–S equations, pressure Poisson equation method, incompressible, compact scheme.

1. Introduction

To solve incompressible fluid flow problems numerically, we can use the primitive variables velocity and pressure as well as the vortex and stream function in the control equations. A main difficulty associated with the solution of the incompressible Navier–Stokes equations in velocity–pressure formulation is the presence of the constraint divV =0,which must be satisfied at any time, does not allow the use of a simple explicit method that avoids solution of an algebraic system of equations. While using the vortex–stream function equations, the continuity equation is satisfied automatically. Therefore there is no problem of the constraint, but the boundary condition for the vortex is difficult to handle, and it is not easy to apply to three dimensional problems and problems with free surface or other fluid interface.

The methods of solving the difficulty of the velocity must satisfy the divergence–free constraint are: the artificial compressibility method^{[3][17]}; the pressure Poisson equation method^{[4][16]}; and the projection method^{[17]}; the divergence–free scheme^{[5]−[7]},and etc. For the unsteady problems, Chorin’s artificial compressibility method^{[3]} can be written as: for each time step:

(1) calculate the velocity V _{k+1}^{n+1} (using the momentum equations with the last step

pressure p_{k}^{n+1}) (1.1) (2) p_{k+1}^{n+1} = p_{k}^{n+1} − λ^{−1}div_{h}V _{k+1}^{n+1}

where λ > 0 is a small parameter that needs adjust in the computation.

The pressure Poisson equation method is^{[4][16]}: taking the divergence of the momentum equations, eliminating the divergence of the (n+1)−step velocity, we get the pressure Poisson equation which is used to replace the continuity equation.

�\underline{2 YU XIN}�

This paper improves the pressure Poisson equation method to an iterative algorithm, using the increment of the pressure, the difference of two successive pressures (p_{k}^{n+1} and p_{k+1}^{n+1}), instead of the pressure p^{n+1} as an unknown variable, see (2.7). It has the following advantages:

1Author’s E-mail: yu@cerse.psu.edu home page: http://www.cerse.psu.edu/yu/

2Supported: Fund for Excellent Researchers Coming Back from Foreign Countries of Chinese Academy of Sciences (1999); China Key Laboratory for CFD (#2000–16)

3Further discusses,proofs and deductions will be in http://www.cerse.psu.edu/yu/ipp/

1. it can ensure the discrete continuity equation satisfied as exactly as expected(see (2.8));

2. ∇_{H}^{2} in the Poisson equation (2.7) needs not adopt a high order accurate operator, e. g., for a 2D 4th order accurate compact scheme, ∇_{H}^{2} can adopt the 5–point central difference. Chorin method (1.1) corresponds to (2.7) with ∇_{H}^{2} = −λ/Δt;

3. comparing with Chorin method (1.1), it convergent much faster;

4. it can be applied to three dimensional problems directly. (for a 3D problem, ∇_{H}^{2} can adopt a 7–point central difference);

5. it can be extended to finite element schemes.

Constructing a difference scheme with a thrice spline function interpolation may improve computation accuracy without increasing the mesh points, see [9]. The compact schemes were developed on the basis of the trice spline function difference scheme. Its essential idea is using trice spline functions to make difference of the spatial derivatives. Liu Hong^{[1]} used a second order accurate scheme in [11] to calculate a driven flow in a square cavity, the results agree well with the experiment only in the case of low Re. He^{[1]} constructed a compact scheme according to [10], applied it to solving the driven flow problem, in the calculation he found oscillations near the left–upper corner. When Re much larger, there were oscillations near a lower corner too, so that the computations could not keep on. This shows that the central difference non–staggered mesh compact scheme produces non–physical numerical oscillations at where flow parameter varies acutely. He^{[1]} applied the upwind technique to the non–staggered mesh compact difference scheme to solve the incompressible flow, solved the above oscillations problem successfully. His main results in [1] are published in [19].

The staggered mesh compact scheme proposed by us does not produce the above oscillations in the same calculations of the driven flow problem without using the upwind technique. The three–point central difference scheme is of the fourth order accuracy, while the upwind three–point compact scheme is of third order accuracy. This shows a good quality of the staggered mesh schemes.

Since unsteady phenomenon is concerned, it is necessary to approximate the time derivative with high order accuracy. We use the fourth order accurate Runge–Kutta method, see §3.4.

This paper is a revised version of [18].

2. The Iterative Pressure Poisson Equation Method

The unsteady viscous Navier–Stokes equations are:

\frac{∂V}{∂t} + A(V ) + ∇p = 0, in Ω , ( where A(V ) = (V · ∇)V − ν∇^{2}V )  (2.1)

The incompressible continuity equation is:

div V = 0, in Ω  ( in this paper: V = (u, v)^{T} )  (2.2)

We consider an explicit discrete form of (2.1)(2.2):

\frac{V^{n+1} − V^{∗}}{Δt} + ∇_{h}p^{n+1} = 0, ( where V^{∗} = V^{n} − ΔtA_{h}(V^{n}) ),  (2.3)

�\underline{Iterative Pressure Poisson Equation Method for Solving the Unsteady Incompressible N–S Equations 3}�

div_{h} V^{n+1} = 0,  ( V^{n+1} = (u^{n+1}, v^{n+1})^{T} for 2D problems )  (2.4)

(We will describe a particular scheme — a compact difference scheme in §3.1). Then the pressure Poisson equation can be written as:

∇_{h}^{2}p^{n+1} = \frac{1}{Δt}div_{h}V^{∗}  (2.5)

where V^{∗} = V^{n} −ΔtA_{h}(V^{n}), Δt is the time step. Both LHS and RHS of (2.5) have various forms, corresponding to a kind of different pressure Poisson equation methods. [16] proposed a class of pressure Poisson equation schemes that satisfy the equivalency ^{4}.

The iterative pressure Poisson equation method:

Take n−step values as initial: V _{0}^{n+1} = V^{n}, p_{0}^{n+1} = p^{n}, calculate V _{k+1}^{n+1}, p_{k+1}^{n+1}, k = 0, 1, 2, ... iteratively:

(1) calculate the velocity V _{k+1}^{n+1}:

\frac{V _{k+1}^{n+1} − V^{∗}}{Δt} + ∇_{h}p_{k}^{n+1} = 0, ( where V^{∗} = V^{n} − ΔtA_{h}(V^{n}) ),  (2.6)

(2) solve the pressure p_{k+1}^{n+1} with an approximate Poisson equation (first solve p_{k+1}^{n+1} − p_{k}^{n+1} as one unknown variable):

∇_{H}^{2}(p_{k+1}^{n+1} − p_{k}^{n+1}) = \frac{1}{Δt}div_{h}V _{k+1}^{n+1}  (2.7)

(3) set V^{n+1} = V _{k+1}^{n+1} when the following inequality valid

||div_{h}V _{k+1}^{n+1}|| ≤ ϵ  (2.8)

where ϵ > 0 is a small quantity given beforehand, || · || is a norm. ϵ can be O(h^{4}) when (2.3)(2.4) are of fourth order accuracy.

div_{h} in (2.7) (2.8) is a part of the original scheme, while ∇_{H}^{2} in (2.7) can adopt a simple difference operator. e. g., for a 2D compact scheme, ∇_{H}^{2} may adopt a 5–point central difference.

The iterative algorithm above is a means of solving the original scheme, such as the compact scheme. It does not change the numerical solution. (In this point, it just likes the Gauss–Seidel method for a system of linear equations). The divergence–free scheme, (the dimensional reduction method in [5]–[7]), is also a means of solving the original scheme.

We write the iterative method above to a general form: for linear equations

Lp = b  (2.9)

iterate to solve {p_{k}} from:

L_{1}(p_{k+1} − p_{k}) = b − Lp_{k}, (k = 0, 1, 2, ...)  (2.10)

where L_{1} is an approximate form of L. Generally L_{1} is simpler than L, (2.10) is easily to be solved than (2.9). Now we briefly analyze the convergency. (2.10) can be written as:

p_{k+1} = L_{1}^{−1}(L_{1} − L)p_{k} + L_{1}^{−1}b  (2.11)

4The continued form (in contrast with the discrete form in [16]) of the equivalency is: the differential equation div(V − αNS) = 0 (corresponding to (3.3) in [16] ) and (2.1) are equivalent to (2.2)(2.1), here NS is the left hand side of (2.1), α is an arbitrary constant.

div(V − αNS) = 0 leads to a pressure Poisson equation.

the convergence condition is ρ(L_{1}^{−1}(L_{1}−L)) < 1. Comparing (2.6)–(2.10), L =div_{h}∇_{h}, L_{1} = ∇_{H}^{2}. Therefore ρ(L_{1}^{−1}(L_{1} − L)) < 1 should be no problem generally.

3. Compact Difference Schemes (2D, staggered mesh, fourth order)

We consider the two dimensional unsteady viscous incompressible N–S equations:

\frac{∂V}{∂t} + A(V ) + ∇p = 0, in Ω ,  (3.1)
div V = 0,  in Ω  (3.2)

where A(V ) = (uu_{x} + vu_{y} − ν(u_{xx} + u_{yy}), uv_{x} + vv_{y} − ν(v_{xx} + v_{yy}))^{T} , ∇p = (p_{x}, p_{y})^{T} , divV =u_{x} + v_{y} 3.1 Fourth order accurate compact schemes on a 2D staggered mesh

For this kind of schemes, the derivatives, as well as the velocity and the pressure themselves, are employed to be unknowns of the difference equations. For explicit schemes and those that have no implicit compact difference, we can solve the derivatives u_{x}, u_{y} from u first. The solution pattern is similar to the ADI method: in x and y directions, solve the partial derivatives in x and y directions respectively. For a uniform mesh: consider a finite difference scheme of (3.1)(3.2),

\frac{V^{n+1} − V^{n}}{Δt} + A_{h}(V^{n}) + ∇_{h}p^{n+1} = 0,  (3.3)
div_{h} V^{n+1} = 0,  (3.4)

where V^{n} = (u^{n}, v^{n})^{T}; A_{h}, ∇_{h}, div_{h} are finite difference forms of A, ∇, div :

(1) finite difference forms of the first derivatives u_{x}, u_{y} in A_{h}(V ) are:

\frac{u′_{i−1,j− \frac{1}{2}} + 4u′_{i,j− \frac{1}{2}} + u′_{i+1,j− \frac{1}{2}}}{6} = \frac{u_{i+1,j− \frac{1}{2}} − u_{i−1,j− \frac{1}{2}}}{2Δx} , (i ≥ 1, j ≥ 1),  (3.5)_{1}

\frac{u′_{i,j− \frac{1}{2}} + 4u′_{i,j+ \frac{1}{2}} + u′_{i,j+ \frac{3}{2}}}{6} = \frac{u_{i,j+ \frac{3}{2}} − u_{i,j− \frac{1}{2}}}{2Δy} , (i ≥ 1, j ≥ 1),  (3.5)_{2}

v_{y}, v_{x} in A_{h}(V ) are similar to (3.5),

(2) finite difference forms of the second derivatives u_{xx}, u_{yy} in A_{h}(V ):

\frac{u′′_{i−1,j− \frac{1}{2}} + 10u′′_{i,j− \frac{1}{2}} + u′′_{i+1,j− \frac{1}{2}}}{12} = \frac{u_{i−1,j− \frac{1}{2}} − 2u_{i,j− \frac{1}{2}} + u_{i+1,j− \frac{1}{2}}}{(Δx)^{2}} , (i ≥ 2, j ≥ 1), (3.6)_{1}

\frac{u′′_{i,j− \frac{1}{2}} + 10u′′_{i,j+ \frac{1}{2}} + u′′_{i,j+ \frac{3}{2}}}{12} = \frac{u_{i,j− \frac{1}{2}} − 2u_{i,j+ \frac{1}{2}} + u_{i,j+ \frac{3}{2}}}{(Δy)^{2}} , (i ≥ 1, j ≥ 1),  (3.6)_{2}

v_{yy}, v_{xx} in A_{h}(V ) are similar to (3.6),

(3) for the derivative p_{x} in ∇_{h}p,

\frac{p′_{i−1,j− \frac{1}{2}} + 22p′_{i,j− \frac{1}{2}} + p′_{i+1,j− \frac{1}{2}}}{24} = \frac{p_{i+ \frac{1}{2} ,j− \frac{1}{2}} − p_{i− \frac{1}{2} ,j− \frac{1}{2}}}{Δx} , (i ≥ 2, j ≥ 1),  (3.7)

similar for p_{y}, (4) for u_{x} in div_{h}V ,

\frac{u′_{i− \frac{1}{2} ,j− \frac{1}{2}} + 22u′_{i+ \frac{1}{2} ,j− \frac{1}{2}} + u′_{i+ \frac{3}{2} ,j− \frac{1}{2}}}{24} = \frac{u_{i+1,j− \frac{1}{2}} − u_{i,j− \frac{1}{2}}}{Δx} , (i ≥ 1, j ≥ 1),  (3.8)

similar for v_{y} in div_{h}V in v_{y},

(5) for u in uv_{x} in (V · ∇)V , the interpolation is employed:

u_{i− \frac{1}{2} ,j} = \frac{ū_{i−1,j− \frac{1}{2}} + ū_{i,j− \frac{1}{2}} + ū_{i−1,j+ \frac{1}{2}} + ū_{i,j+ \frac{1}{2}}}{4} , (i ≥ 1, j ≥ 1),  (3.9)_{1}

ū_{i,j− \frac{1}{2}} = u_{i,j− \frac{1}{2}} − \frac{1}{8}[(Δx)^{2}u_{xx} + (Δy)^{2}u_{yy}]_{i,j− \frac{1}{2}} , (i ≥ 0, j ≥ 1),  (3.9)_{2}

where u_{xx}, u_{yy} can adopt the result of (3.6), and for v in vu_{y}, similar to above. (3.9)_{3} (6) for V_{t}: (3.3) adopts the first order accurate difference. A fourth order accurate Runge– Kutta method is given in §3.4. 3.2 difference formulations on the boundary

V |_{Γ} = V_{Γ} = (u^{Γ}, v^{Γ})^{T} on Γ = ∂Ω. Now we mainly describe difference formulations near the left boundary x = 0. (1)_{1} (u_{x})_{0,j− \frac{1}{2}} = (−v_{y}^{Γ})_{0,j− \frac{1}{2}} , (j ≥ 1) (from the continuity equation (3.2)), (3.10) (1)_{2} u_{y} at y = \frac{3}{4}Δy and v_{x} at x = \frac{3}{4}Δx :

\frac{3u′_{i, \frac{1}{2}} + u′_{i, \frac{3}{2}}}{4} = \frac{u_{i, \frac{3}{2}} − u_{i,0}}{\frac{3}{2}Δy} , (i ≥ 1); \frac{3v′_{\frac{1}{2} ,j} + v′_{\frac{3}{2} ,j}}{4} = \frac{v_{\frac{3}{2} ,j} − v_{0,j}}{\frac{3}{2}Δx} , (j ≥ 1),  (3.11)

here u_{i,0} = (u^{Γ})_{i,0} = u^{Γ}|_{x=iΔx,y=0}, v_{0,j} = (v^{Γ})_{0,j} = v^{Γ}|_{x=0,y=jΔy} (2)_{1} u_{xx} at x = \frac{10}{9}Δx:

\frac{8u′′_{1,j− \frac{1}{2}} + u′′_{2,j− \frac{1}{2}}}{9} = \frac{9u_{0,j− \frac{1}{2}} − 16u_{1,j− \frac{1}{2}} + 7u_{2,j− \frac{1}{2}}}{6(Δx)^{2}} + \frac{u′_{0,j− \frac{1}{2}}}{3Δx} , (j ≥ 1)  (3.12)

where u′_{0,j− \frac{1}{2}} = −(v_{y}^{Γ})_{0,j− \frac{1}{2}} = −v_{y}^{Γ}|_{x=0,y=(j− \frac{1}{2} )Δy}. v_{yy} at y = \frac{10}{9}Δy is similar to (3.12),

(2)_{2} u_{yy} at y = \frac{2}{3}Δy:

\frac{5u′′_{i, \frac{1}{2}} + u′′_{i, \frac{3}{2}}}{6} − \frac{1}{48}(u′′_{i, \frac{1}{2}} − 2u′′_{i, \frac{3}{2}} + u′′_{i, \frac{5}{2}} ) = \frac{4}{3} · \frac{2u_{i,0} − 3u_{i, \frac{1}{2}} + u_{i, \frac{3}{2}}}{(Δy)^{2}} , (i ≥ 1),  (3.13)

where u_{i,0} = u_{i,0}^{Γ}, v_{xx} at x = \frac{2}{3}Δx is similar to (3.13),

(3) p_{x} at x = Δx:

p′_{1,j− \frac{1}{2}} + \frac{p′_{1,j− \frac{1}{2}} − 2p′_{2,j− \frac{1}{2}} + p′_{3,j− \frac{1}{2}}}{24} = \frac{p_{\frac{3}{2} ,j− \frac{1}{2}} − p_{\frac{1}{2} ,j− \frac{1}{2}}}{Δx} , (j ≥ 1)  (3.14)

p_{y} at y = Δy similar,

(4) u_{x} in divV at x = \frac{1}{2}Δx:

\frac{2u′_{0,j− \frac{1}{2}} + 15u′_{\frac{1}{2} ,j− \frac{1}{2}} + u′_{\frac{3}{2} ,j− \frac{1}{2}}}{18} = \frac{u_{1,j− \frac{1}{2}} − u_{0,j− \frac{1}{2}}}{Δx} , (j ≥ 1),  (3.15)

where u′_{0,j− \frac{1}{2}} = −(v_{y}^{Γ})_{0,j− \frac{1}{2}} . v_{y} at y = \frac{1}{2}Δy is similar to (3.15), (5) boundary values of u_{xx}, u_{yy} in (3.9)_{2}:

(u_{xx})_{0,j− \frac{1}{2}} = 2(u_{xx})_{\frac{1}{2} ,j− \frac{1}{2}} − (u_{xx})_{1,j− \frac{1}{2}} , (u_{xx})_{\frac{1}{2} ,j− \frac{1}{2}} = \frac{(u_{x})_{1,j− \frac{1}{2}} − (u_{x})_{0,j− \frac{1}{2}}}{Δx} , (3.16)_{1}
(u_{yy})_{0,j− \frac{1}{2}} = (u_{yy}^{Γ})_{0,j− \frac{1}{2}} ,  (3.16)_{2}

where (u_{x})_{1,j− \frac{1}{2}} can be got from (3.5). (u_{x})_{0,j− \frac{1}{2}} = −(v_{y}^{Γ})_{0,j− \frac{1}{2}} . 3.3 An upwind compact difference scheme

Alter the formulation (3.5)_{1} to

\frac{1}{3}u′_{i−1,j− \frac{1}{2}} + \frac{2}{3}u′_{i,j− \frac{1}{2}} = \frac{u_{i+1,j− \frac{1}{2}} + 4u_{i,j− \frac{1}{2}} − 5u_{i−1,j− \frac{1}{2}}}{6Δx} , (if v_{i,j− \frac{1}{2}} ≥ 0)( \begin{matrix}i ≥ 1 \\ j ≥ 1 \end{matrix} ) (3.17)_{1}

\frac{2}{3}u′_{i,j− \frac{1}{2}} + \frac{1}{3}u′_{i+1,j− \frac{1}{2}} = \frac{5u_{i+1,j− \frac{1}{2}} − 4u_{i,j− \frac{1}{2}} − u_{i−1,j− \frac{1}{2}}}{6Δx} , (if v_{i,j− \frac{1}{2}} < 0)( \begin{matrix}i ≥ 1 \\ j ≥ 1 \end{matrix} ) (3.17)_{2}

(where v_{i,j− \frac{1}{2}} see (3.9)_{3}). alter(3.5)_{2} similarly. Other formulations are same with (3.6)– (3.16).

The difference between (3.17) and (1.1.26)(1.1.27) in [1], i. e. (2.8)(2.9) in [19], is F (= u′) has no superscripts + and –. (I − 4 + 5I^{−1} should be I + 4 − 5I^{−1} in (2.8) of [19]). 3.4 The Runge–Kutta method (For unsteady problem, time direction)

Define f(V ) = A_{h}(V )+∇_{h}p, where p = p(V ) satisfies div_{h}(A_{h}(V )+∇_{h}p) = 0 . Thus f is a function of V . The fourth order Runge–Kutta formulation for solving V_{t} +f(V ) = 0 is:

\frac{V^{n+1} − V^{n}}{Δt} + \frac{k_{1} + 2k_{2} + 2k_{3} + k_{4}}{6} = 0,  (3.18)

where

V^{(1)} = V^{n} − \frac{Δt}{2}k_{1}, V^{(2)} = V^{n} − \frac{Δt}{2}k_{2}, V^{(3)} = V^{n} − Δtk_{3},

k_{1} = f(V^{n}), k_{2} = f(V^{(1)}), k_{3} = f(V^{(2)}), k_{4} = f(V^{(3)})

For intermediate boundary conditions, we have proposed the following formulations in the journal ”Mathematica Numerica Sinica” (Vol. 20, No. 1, 1998, page 56), (in the paper, the last ”–” in the formulation for V^{(2)} should be altered by ”+”):

V^{(1)} = V^{n}+\frac{Δt}{2}(\frac{∂V}{∂t} |_{t=nΔt}), V^{(2)} = 2(V |_{t=(n+ \frac{1}{2} )Δt})−V^{(1)}, V^{(3)} = (V |_{t=(n+1)Δt}), (on Γ)

4. Calculations of the Driven Flow in a Square Cavity

We consider the viscous incompressible fluid flow problem driven by the shearing force in a two dimensional unit square cavity. The control equations adopt the unsteady Navier– Stokes equations (3.1)(3.2). Computation area: 0 < x < 1, 0 < y < 1. Boundary conditions:

V = V_{Γ} = \{ \begin{matrix}(−1, 0)^{T} & when y = 1 \\ (0, 0)^{T} & when y = 0 or x = 0 or x = 1 \end{matrix}  (4.1)

Δx = Δy = 1/N with N × N. ∇_{H}^{2} in (2.7) adopts the five–point central difference.

To solve (2.7), use the multigrid method one loop (grid from fine to coarse, then from coarse to fine), thus the algorithm likes an integrated multigrid procedure. (for such ”multigrid method”, in the ”first–grid”, ∇_{h}^{2} is the ”fine grid”, ∇_{H}^{2} is the ”coarse grid”).

Figure 1 presents the calculation results of the iterative pressure Poisson equation method with the upwind staggered mesh compact difference scheme (3.17)(3.6)–(3.9), Re = 1/ν = 5000, N = 128. The boundary conditions and the method for time derivative adopt those in [18] (in fact, figure 1 and figure 2 are computation examples in [18]), Figure 1b and 1c show that our results agree well with those in [14]. Figure 2 presents calculation results of the iterative pressure Poisson equation method with the non–upwind staggered mesh compact difference scheme (3.5)–(3.9), Re = 1/ν = 10000, N = 256. The boundary conditions and the method for time derivative adopt those in [18]. Figure 2a1 and 2a2 are close to figure 7 in [1] (i. e., figure 5 in [19]).

Now we use the boundary conditions in §3.2 and the Runge–Kutta method in §3.4 (the intermediate boundary conditions adopt V^{(1)} = V^{(2)} = V_{Γ}|_{t=(n+ \frac{1}{2} )Δt}, V^{(3)} = V_{Γ}|_{t=(n+1)Δt}) on Γ, the iterative pressure Poisson equation method in §2, the non–upwind staggered mesh compact scheme in §3.1, Re = 10^{4}, the initial condition: V =0 on inner grids, V = V_{Γ}, N = 256 × 256, Δt = 0.6Δx, set ϵ = 10^{−5} for (2.8). Figure 3 presents the stream fields for time steps n = 73840, 86800, 121840.

To verify the equivalency of the iterative method (2.6)(2.7) and the original equations (2.3)(2.4), from t = 80640 × 0.6/256, we alter ϵ in (2.8) to 10^{−13}, the calculations continued. This shows the iterative method replaces the original equations equivalently. When the test calculations for verifying the equivalency performed 100 time steps, the procedure with ϵ = 10^{−4} finished 646 steps within the same computing time. besides, we calculated for the following three cases, all calculations continued (we calculated 10 time steps for each case):

(1)from t = 80640 × 0.6/256, we alter ϵ in (2.8) to 10^{−13}, Re to 10^{3}; (2)set N=64,Re = 10^{4}, ϵ = 10^{−14}, calculate from the beginning; (3)set N=64,Re = 10^{3}, ϵ = 10^{−13}, calculate from the beginning.

Fig.1a. Streamlines, 128 × 128 grids, Re = 5000, the upwind compact difference scheme Fig.1b. velocity u profile,x = 0.5 —-: our results, ⋯: results in [14] Fig.1c. velocity v profile,y = 0.5 —-: our results, ⋯: results in [14] Fig.2a1. Streamlines, 256 × 256 grids, the compact difference scheme, Re = 10000, t = t_{0} + 12.7490

Fig.2a2. Streamlines, 256 × 256 grids, the compact difference scheme, Re = 10000, t = t_{0} + 13.9111 Fig.2b1. velocity u profile,x = 0.5 Fig.2c1. velocity v profile,y = 0.5 Fig.3a. Streamlines, 256 × 256 grids, Re = 10000, t = 73840 × 0.6/256 values of stream function: − 0.002, −0.001, −0.0003, 0, 0.0003, 0.01, 0.02, ⋯, 0.09

Fig.3b. Streamlines, 256 × 256 grids, Re = 10000, t = 86800 × 0.6/256 values of stream function: − 0.002, −0.001, −0.0003, 0, 0.0003, 0.01, 0.02, ⋯, 0.09, 0.098 Fig.3c. Streamlines, 256 × 256 grids, Re = 10000, t = 121840 × 0.6/256 stream function values: − 0.002, −0.001, −0.0003, 0, 0.0003, 0.01, 0.02, ⋯, 0.09, 0.1

Further discusses,proofs and deductions will be put in our website: http://www.cerse.psu.edu/yu/ipp/ We would like to express our gratitude to professor WANG Famin, professor SHU Chi-wang,

professor SHEN Jie, professor CHEN Longqing, professor FU Dexun, professor MA Yanwen, professor NIU Jiayu, professor PENG Shiliu, professor HUANG Dun

References

[1]LIU Hong, Study on high order accurate difference schemes and direct numerical simulations of the unsteady vortex, Master thesis of the Institute of Mechanics, Chinese Academy of Sciences, Beijing, 1992, (main results are published in [19])

[2]P. J. Roache, Computational Fluid Dynamics, Hermosa Publishers, Albuquerque,1972

[3]A. J. Chorin, A numerical method for solving incompressible viscous flow problems, J.

Comput. Phys., 2: 1, 1967, 12–26 [4]F. H. Harlow, J. E. Welch, Numerical calculation of time–dependent viscous incompress-

ible flow of fluid with free surface, Phys. Fluids, 8:12, 1965, 2182–2189 [5]YU Xin, Divergence–free schemes for equations of incompressible fluid dynamics, Doc-

toral dissertation of the Institute of Mechanics, Chinese Academy of Sciences, Beijing, 1990

[6]C. Cuvelier, et al., Finite Element Methods and Navier–Stokes Equations, D. Reidel

Publishing Company, 1986 [7]YU Xin, A dimensional reduction method for incompressible fluid dynamics: I. The

essential algorithm, Chinese Journal of Computational Physics, 2:3, 1985, 337–346; II. The basis of V^{h} and the error estimates, Chinese Journal of Computational Physics, 3:2, 1986, 217–226; III. Second order accurate finite element schemes, Chinese Journal of Computational Physics, 5:2, 1988, 211–220; IV. Three dimensional Navier–Stokes equations, Chinese Journal of Computational Physics, 6:1, 1989, 104–116

[8]LIU Hong, FU Dexun, MA Yanwen, Hopf bifurcation of the driven flow in a square cavity, Theories, Methods and Applications of Computational Fluid Mechanics, Science Press, Beijing, 1992, 267–270

[9]S. G. Rubin, R. A. Graves,Jr, Viscous flow solutions with a cubic spline approximation,

Computers and Fluids, 3:1, 1975, 1–36

[10] N. Baba, H. Miyata, H. Kajitani, Journal of the Society of Naval Architects of Japan, 159, 1987, 33

[11] S. Abdallah, Numerical solutions for the incompressible Navier–Stokes equations in primitive variables using a non–staggered grid II, J. Comput. Phys., 70, 1987, 193–202

[12] FU Dexun, MA Yanwen, A Collection of Technical Papers, International Symposium on Computational Fluid Dynamics, 1989

[13] C. H. Bruneau, C. Jouron, An efficient scheme for solving steady incompressible Navier– Stokes equations, J. Comput. Phys., 89, 1990, 389–413

[14] U. Ghia, K. N. Ghia, C. T. Shin, High–Re solutions for incompressible flow using the Navier–Stokes equations and a multigrid method, J. Comput. Phys., 48,(1982),387–411

[15] YU Xin, Numerical solutions for the incompressible Navier–Stokes equations, Handbook for Engineering Applied Mechanics, Chapter 7 of Computational Mechanics, Shanxi People Press, Xi’an, 1996

[16] YU Xin, A staggered mesh compact difference scheme and a pressure–Poisson–equation that satisfies the equivalency, Chinese J. of Numer. Math. and Appl. 19: 2, 1997, 73–81; (Original in Chinese: Mathematica Numerica Sinica, 19:1, 1997, 83–90)

[17] R. Peyret, T. D. Taylor, Computational Methods for Fluid Flow, Springer–Verlag, 1983

[18] YU Xin, Iterative pressure Poisson equation method for solving unsteady incompressible N–S equations, Proc. of Beijing Symposium on Computational Fluid Mechanics, 5, Beijing, 1993, 127–136

[19] LIU Hong, FU Dexun, MA Yanwen, Upwind compact schemes and direct numerical simulations of the driven flow in a square cavity, Science in China, Series A, 23:6, 1993, 657–665