\documentstyle[twoside,headrule]{carticle}
%---------------------- Page Format --------------------------
\headsep 0.5 true cm
\topmargin 0pt
\oddsidemargin 0pt
\evensidemargin 0pt
\textheight 21 true cm
\textwidth 14 true cm
%\parskip 0.2cm
\parindent 2\ccwd
\nofiles
%-------------------------------------------------------------------
\def\pmb#1{\mbox{\boldmath $#1$}}
\def\V{\pmb V}
\def\A{\pmb A}
\begin{document}

%\renewcommand{\baselinestretch}{1.0}
%\renewcommand{\arraystretch}{1.5}
%\newcommand{\dfrac}[2]{\frac{\displaystyle #1}{\displaystyle #2}}
%\def\deq{\mathrel{\mathop{=}\limits^{\triangle}}}
%\def\R{\mbox{I}\!\mbox{R}}
%\def\p{\partial}
%\def\supp{\mbox{supp}}
%\def\disp{\displaystyle}
%\def\sect#1{\section*{\leftline{\large\bf\heiti #1}}}
%\def\jielun#1{\bigskip\noindent{\bf\heiti #1}\bgroup\kaishu\it }
%\def\endjielun{\egroup\par\bigskip}
%???????? I do not know \tenrm ??????????
\def\tenrm{}
%-------------------------------------------------------------------
\def\oddhead{{\protect\footnotesize
	Iterative Pressure Poisson Equation Method for Solving
	the Unsteady Incompressible N--S Equations \hfill}}
\def\evenhead{{\footnotesize \hfill YU XIN}}
%--------------------------------------------------------------------------------------------------
\thispagestyle{empty}
\vspace*{-13mm}

\noindent
\hbox to \textwidth{\footnotesize 
	Chinese Journal of Numerical Mathematics and Applications \hfill}

\hbox to \textwidth{\footnotesize 
	Vol. {\bf 23}, No. 4, Nov., 2001, 447--456 \hfill}
%\par\noindent
%\rule[2mm]{\textwidth}{0.1pt}\hspace*{-\textwidth}\rule[1.5mm]{\textwidth}{0.1pt}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\vspace*{1cm}

\begin{center}
{\Large 		Iterative Pressure Poisson Equation Method for Solving
			the Unsteady Incompressible N--S Equations
}\vskip1.7em

YU Xin	\footnote{Author's E-mail: yu@cerse.psu.edu \quad
		home page: http://www.cerse.psu.edu/yu/}

({\it LHD, Institute of Mechanics, Chinese Academy of Sciences, Beijing})\vskip1em

{\bf Abstract}		\footnote{Supported: Fund for Excellent Researchers
	Coming Back from Foreign Countries of Chinese Academy of Sciences (1999);
\qquad	China Key Laboratory for CFD (\#2000--16)}
\footnote{Further discusses,proofs and deductions will be in 
http://www.cerse.psu.edu/yu/ipp/}
	\end{center}

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^{n+1}_k$ and $p^{n+1}_{k+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$.

{\bf Key Words}: N--S equations, pressure Poisson equation method,
incompressible, compact scheme.

\vskip2em
\begin{center} {\large \bf
				1. Introduction	}\end{center}

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 div\V=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^{n+1}_{k+1}$
 (using the momentum equations with the last step pressure $p^{n+1}_k$)
\vskip-2em
$$								\eqno{(1.1)}
$$
\vskip-1em
(2) $p^{n+1}_{k+1}=p^{n+1}_k-\lambda^{-1} {\rm div}_h\V^{n+1}_{k+1}$

\noindent
where $\lambda >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.

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^{n+1}_k$ and $p^{n+1}_{k+1})$,
 instead of the pressure $p^{n+1}$ as an unknown variable, see (2.7).
It has the following advantages:

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

2. $\nabla^2_H$ in the Poisson equation (2.7) needs not adopt
a high order accurate operator, e. g., for a 2D 4th order accurate
compact scheme, $\nabla^2_H$ can adopt the 5--point central difference.
Chorin method (1.1) corresponds to (2.7)
with $\nabla^2_H=-\lambda /\Delta 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, 
	$\nabla^2_H$ 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 \S 3.4.

This paper is a revised version of [18].
	\vskip2em
\begin{center} {\large \bf
			2. The Iterative Pressure Poisson Equation Method
	} \end{center}

The unsteady viscous Navier--Stokes equations are:
$$
\frac{\partial\V}{\partial t}+\A(\V)+\pmb\nabla p={\bf 0},
    \hbox{ in } \Omega\hbox{ },\quad(\hbox{ where }
\A(\V)=(\V\cdot\pmb\nabla)\V-\nu\nabla^2\V~)			\eqno{(2.1)}
$$
\noindent
The incompressible continuity equation is:
$$
\hbox{div}~\V=0, \quad\hbox{ in }\Omega \hbox{ }\quad
\quad(\mbox{ in this paper: }\V=(u,v)^T~)\quad   		\eqno{(2.2)}
$$
\noindent
We consider an explicit discrete form of (2.1)(2.2):

$$
\frac{\V^{n+1}-\V^{*}}{\Delta t}+\pmb\nabla_hp^{n+1}=0,
    \quad\hbox{ ( where }\V^*=\V^n-\Delta t\A_h(\V^n)~),	\eqno{(2.3)}
$$

$$
\hbox{div}_h~\V^{n+1}=0,\quad\quad(~\V^{n+1}=(u^{n+1},v^{n+1})^T~
\mbox{for 2D problems })\quad				\eqno{(2.4)}
$$

\noindent
(We will describe a particular scheme --- a compact difference scheme in \S 3.1).
Then the pressure Poisson equation can be written as:

$$
	\nabla^2_h p^{n+1}=\frac 1{\Delta t}{\rm div}_h\V^*	\eqno{(2.5)}
$$

\noindent
where $\V^*=\V^n-\Delta t \A_h(\V^n),\Delta 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
\footnote
{The continued form (in contrast with the discrete form in [16]) of the equivalency is:
  the differential equation
div$(\V-\alpha \mbox{\bf NS})=0$ (corresponding to (3.3) in [16] ) and
(2.1) are equivalent to (2.2)(2.1), here {\bf NS} is the left hand side of (2.1), 
$\alpha$ is an arbitrary constant.

div$(\V-\alpha \mbox{\bf NS})=0$ leads to a pressure Poisson equation.}.

{\bf The iterative pressure Poisson equation method:}

Take $n-$step values as initial:  $\V^{n+1}_0=\V^n,p^{n+1}_0=p^n$,
 calculate $\V^{n+1}_{k+1}, p^{n+1}_{k+1},k=0,1,2,...$ iteratively:

\noindent
(1) calculate the velocity $\V^{n+1}_{k+1}$:
$$
\frac{\V^{n+1}_{k+1}-\V^{*}}{\Delta t}+\pmb\nabla_hp^{n+1}_k=0,
    \quad\hbox{ ( where }\V^*=\V^n-\Delta t\A_h(\V^n)~),     	\eqno{(2.6)}
$$
\noindent
(2) solve the pressure $p^{n+1}_{k+1}$ with an approximate Poisson equation
(first solve $p^{n+1}_{k+1}-p^{n+1}_k$ as one unknown variable):

$$
%\nabla^2_H(
\nabla ^2_H(p^{n+1}_{k+1}-p^{n+1}_k)
=\frac 1{\Delta t}{\rm div}_h\V^{n+1}_{k+1}			\eqno{(2.7)}
$$

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

$$
||{\rm div}_h \V^{n+1}_{k+1}||\leq\epsilon			\eqno{(2.8)}
$$
where $\epsilon > 0$ is a small quantity given beforehand,
 $||\cdot||$ is a norm. $\epsilon$ 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 $\pmb\nabla^2_H$
in (2.7) can adopt a simple difference operator.
e. g.,
for a 2D compact scheme, $\nabla^2_H$ 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.

%While solving (2.7), no need to calculate very carefully.
%Several iterations of an iterative method are enough.

%When we solve a steady problem with a time dependent method,
%we need only iterative one or two times for $k$.

We write the iterative method above to a general form: for linear equations
$$
			Lp=b				\eqno{(2.9)}
$$
iterate to solve $\{p_k\}$ from:
$$
	L_1(p_{k+1}-p_k)=b-Lp_k, \quad (k=0,1,2,...)		\eqno{(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			\eqno{(2.11)}
$$
the convergence condition is $\rho(L^{-1}_1(L_1-L))<1$. Comparing (2.6)--(2.10),
$L=$div$_h\pmb\nabla_h,L_1=\nabla^2_H$. Therefore $\rho(L^{-1}_1(L_1-L))<1$
should be no problem generally.

\vskip2em
\begin{center} {\large {\bf
	3. Compact Difference Schemes (2D, staggered mesh, fourth order)}}
 \end{center}

We consider the two dimensional unsteady viscous incompressible N--S equations:
$$
\frac{\partial\V}{\partial t}+\A(\V)+\pmb\nabla p={\bf 0},
    \hbox{ in } \Omega\hbox{ },					\eqno{(3.1)}
$$
$$
\hbox{div}~\V=0, \quad \quad\hbox{ in } \Omega \hbox{ }\quad\quad	 \eqno{(3.2)}
$$
where $\A(\V)=(uu_x+vu_y-\nu(u_{xx}+u_{yy}),uv_x+vv_y-\nu(v_{xx}+v_{yy}))^T$ ,
$\pmb\nabla p=(p_x,p_y)^T,$ div\V=$u_x+v_y$

\noindent
{\bf 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}{\Delta t}+\A_h(\V^n)+\pmb\nabla_hp^{n+1}=0,\quad	\eqno{(3.3)}
$$
$$
\hbox{div}_h~\V^{n+1}=0,\quad\quad\quad\quad\quad\quad		\eqno{(3.4)}
$$
\noindent
where $\V^n=(u^n,v^n)^T; \A_h,\pmb\nabla_h, $ div$_h$
are finite difference forms of $\A,\pmb\nabla$, div :

\noindent
(1) finite difference forms of the first derivatives $u_x,u_y$ in $\A_h(\V)$ are:
$$
\frac{u^{\prime}_{i-1, j-\frac{1}{2}}+4u^{\prime}_{i,j-\frac{1}{2}}
+u^{\prime}_{i+1, j-\frac{1}{2}}}{6}=
\frac{u_{i+1, j-\frac{1}{2}}-u_{i-1, j-\frac{1}{2}}}{2\Delta x},
(i\geq 1,j\geq 1),						\eqno{(3.5)_1}
$$
$$
\frac{u^{\prime}_{i, j-\frac{1}{2}}+4u^{\prime}_{i,j+\frac{1}{2}}
+u^{\prime}_{i, j+\frac{3}{2}}}{6}=
\frac{u_{i, j+\frac{3}{2}}-u_{i, j-\frac{1}{2}}}{2\Delta y},(i\geq 1,j\geq 1), \eqno{(3.5)_2}
$$

$v_y,v_x$ in $\A_h(\V)$ are similar to (3.5),

\noindent
(2) finite difference forms of the second derivatives $u_{xx},u_{yy}$ in $\A_h(\V)$:
$$
\frac{u^{\prime\prime}_{i-1, j-\frac{1}{2}}+10u^{\prime\prime}
_{i,j-\frac{1}{2}}+u^{\prime\prime}_{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}}}
{(\Delta x)^2},(i\geq 2,j\geq 1),			                            \eqno{(3.6)_1}
$$
$$
\frac{u^{\prime\prime}_{i, j-\frac{1}{2}}+10u^{\prime\prime}
_{i,j+\frac{1}{2}}+u^{\prime\prime}_{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}}}
{(\Delta y)^2},(i\geq 1,j\geq 1),			                            \eqno{(3.6)_2}
$$

$v_{yy},v_{xx}$ in $\A_h(\V)$ are similar to (3.6),

\noindent
(3) for the derivative $p_x$ in $\pmb\nabla_hp$,
$$
\frac{p^{\prime}_{i-1, j-\frac{1}{2}}+22p^{\prime}_{i,j-\frac{1}{2}}
+p^{\prime}_{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}}}
{\Delta x},(i\geq 2,j\geq 1),					\eqno{(3.7)}
$$

similar for $p_y$,

\noindent
(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}}}{\Delta x},(i\geq 1,j\geq 1),			\eqno{(3.8)}
$$

similar for $v_y$ in div$_h \V$ in $v_y$,

\noindent
(5) for $u$ in $uv_x$ in ($\V\cdot\pmb\nabla$)\V,
the interpolation is employed:
$$
u_{i-\frac{1}{2}, j} = \frac{\bar u_{i-1, j-\frac{1}{2}} +
\bar u_{i, j-\frac{1}{2}} +\bar u_{i-1, j+\frac{1}{2}} +
\bar u_{i, j+\frac{1}{2}}}{4},(i\geq 1,j\geq 1),			\eqno{(3.9)_1}
$$
$$
\bar u_{i, j-\frac{1}{2}} =u_{i, j-\frac{1}{2}}
-\frac{1}{8}[(\Delta x)^2u_{xx} + (\Delta
y)^2u_{yy}]_{i,j-\frac{1}{2}},(i\geq 0,j\geq 1),			\eqno{(3.9)_2}
$$

where $u_{xx},u_{yy}$ can adopt the result of (3.6), and for $v$ in $vu_y$,
similar to above.		\hfill				$(3.9)_3$

\noindent
(6) for $\V_t$:  (3.3) adopts the first order accurate difference.
A fourth order accurate Runge--Kutta method is given in \S 3.4.

\noindent
{\bf 3.2~~difference formulations on the boundary}

$\V|_\Gamma=\V_\Gamma=(u^\Gamma,v^\Gamma)^T$ on $\Gamma=\partial \Omega$.
Now we mainly describe difference formulations near the left boundary $x=0$.

\noindent
(1)$_1$~ $(u_x) _{0,j-\frac {1}{2}}=(-v^\Gamma_y)_{0,j-\frac {1}{2}},
~(j\geq 1)$ \qquad  (from the continuity equation (3.2)),		\hfill	(3.10)

\noindent
(1)$_2$~ $u_y$ at $y=\frac {3}{4}\Delta y$ and $v_x$ at $x=\frac {3}{4}\Delta 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}\Delta y},(i\geq 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}\Delta x},(j\geq 1),		\eqno{(3.11)}
$$
here $u_{i,0}=(u^\Gamma)_{i,0}=u^\Gamma|_{x=i\Delta x,y=0},
         v_{0,j}=(v^\Gamma)_{0,j}=v^\Gamma|_{x=0,y=j\Delta y}$

\noindent
(2)$_1$~ $u_{xx}$ at  $x=\frac{10}{9}\Delta x$:
$$
\frac{8u^{\prime\prime}_{1, j-\frac{1}{2}} +
u^{\prime\prime}_{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(\Delta
x)^2}+\frac{u^{\prime}_{0, j-\frac{1}{2}}}{3\Delta x},~(j\geq 1)		\eqno{(3.12)}
$$

where $u'_{0,j-\frac {1}{2}}=-(v^\Gamma_y)_{0,j-\frac {1}{2}}=
-v_y^\Gamma|_{x=0,y=(j-\frac{1}{2})\Delta y}$.
$v_{yy}$ at $y=\frac{10}{9}\Delta y$ is similar to (3.12),

\noindent
(2)$_2$~ $u_{yy}$ at $y=\frac {2}{3}\Delta y$:
$$
\frac{5u^{\prime\prime}_{i,\frac{1}{2}} + u^{\prime\prime}
_{i,\frac{3}{2}}}{6} 
- \frac{1}{48}(u^{\prime\prime}_{ i,\frac{1}{2}}
-2u^{\prime\prime}_{i,\frac{3}{2}}+u^{\prime\prime}_{i,\frac{5}{2}})
=\frac{4}{3}\cdot\frac{2u_{i,0}-3u_{i,\frac{1}{2}}
+u_{i,\frac{3}{2}}}{(\Delta y)^2},(i\geq 1),			\eqno{(3.13)}
$$

where $u_{i,0}=u_{i,0}^\Gamma$,~
$v_{xx}$ at $x=\frac {2}{3}\Delta x$ is similar to (3.13),

\noindent
(3)~ $p_x$ at $x=\Delta 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}}}{\Delta x},~ (j\geq 1)			\eqno{(3.14)}
$$

$p_y$ at $y=\Delta y$ similar,

\noindent
(4)~ $u_x$ in div\V ~ at $x=\frac 12\Delta x$:
$$
\frac{2u^{\prime}_{0, j-\frac{1}{2}}+15u^{\prime}_{\frac{1}{2}, j-\frac{1}{2}}+
u^{\prime}_{\frac{3}{2}, j-\frac{1}{2}}}{18}=
\frac{u_{1, j-\frac{1}{2}} -u_{0, j-\frac{1}{2}}}{\Delta x},(j\geq 1),	\eqno{(3.15)}
$$

where $u'_{0,j-\frac {1}{2}}=-(v^\Gamma_y)_{0,j-\frac{1}{2}}$ . 
$v_y$ at $y=\frac {1}{2}\Delta y$ is similar to (3.15),

\noindent
(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}}}{\Delta x},	\eqno{(3.16)_1}
$$
$$
(u_{yy})_{0,j-\frac{1}{2}}=(u^\Gamma_{yy})_{0,j-\frac{1}{2}},	\eqno{(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^\Gamma_y)_{0,j-\frac{1}{2}}$.

\noindent
{\bf 3.3~ An upwind compact difference scheme}

Alter the formulation (3.5)$_1$ to
$$
\frac{1}{3} u^\prime_{i-1,j-\frac{1}{2}}+\frac{2}{3}u^\prime_{i,j-\frac {1}{2}}
=\frac{u_{i+1,j-\frac{1}{2}}+4 u_{i,j-\frac{1}{2}}-5 u_{i-1,j-\frac{1}{2}}}{6\Delta x},
(\mbox{if }v_{i,j-\frac{1}{2}}\ge 0)
(\begin{array}{l} i\geq 1\\ j\geq 1\end{array})		\eqno{(3.17)_1}
$$
$$
\frac{2}{3} u^\prime_{i,j-\frac{1}{2}}+\frac{1}{3} u^\prime_{i+1,j-\frac{1}{2}}
=\frac{5 u_{i+1,j-\frac{1}{2}}-4 u_{i,j-\frac{1}{2}}- u_{i-1,j-\frac{1}{2}}}{6\Delta x},
(\mbox{if }v_{i,j-\frac{1}{2}}< 0)
(\begin{array}{l} i\geq 1\\ j\geq 1\end{array})		\eqno{(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^\prime)$ has no superscripts + and --.
($I-4+5I^{-1}$ should be $I+4-5I^{-1}$ in (2.8) of [19]).

\noindent
{\bf 3.4~The Runge--Kutta method (For unsteady problem, time direction)}

Define $\pmb f(\V)=\A_h(\V)+\pmb\nabla_hp$,
where $p=p(\V)$ satisfies div$_h(\A_h(\V)+\pmb\nabla_hp)=0$ .
Thus $\pmb f$ is a function of \V. 
The fourth order Runge--Kutta formulation for solving
$\V_t+\pmb f(\V)={\bf 0}$ is:
$$
\frac {\V^{n+1}-\V^n}{\Delta t}
+\frac {\pmb k_1+2\pmb k_2+2\pmb k_3+\pmb k_4}{6}={\bf 0},      \eqno{(3.18)}
$$
\noindent
where
$$
\V^{(1)}=\V^n-\frac {\Delta t}{2}\pmb k_1,~
\V^{(2)}=\V^n-\frac {\Delta t}{2}\pmb k_2,~
\V^{(3)}=\V^n-\Delta t\pmb k_3,~
$$
$$
\pmb k_1=\pmb f(\V^n),~~
\pmb k_2=\pmb f(\V^{(1)}),~~
\pmb k_3=\pmb f(\V^{(2)}),~~
\pmb k_4=\pmb 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{\Delta t}{2} (\frac{\partial\V}{\partial t}{\hbox{\Large $|$}}_{t=n\Delta t}),~~
\V^{(2)}=2(\V|_{t=(n+\frac{1}{2})\Delta t})-\V^{(1)},~~
\V^{(3)}=(\V|_{t=(n+1)\Delta t}),~~(\mbox{on }\Gamma\mbox{})
$$

\vskip2em
\begin{center} {\large \bf
		4. Calculations of the Driven Flow in a Square Cavity}
\end{center}

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_\Gamma=\left\{\begin{array}{ll}
(-1,0)^T  & \mbox{ when }\ \ y=1\\ (0,0)^T & \mbox{ when }\ \ y=0
\mbox{ or }x=0\mbox{ or }x=1
\end{array}\right.
						\eqno{(4.1)}
$$
$\Delta x=\Delta y=1/N$ with $N\times N$.
$\nabla_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", $\nabla^2_h$ is the "fine grid",
$\nabla^2_H$ 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/\nu=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/\nu=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 \S 3.2 and the Runge--Kutta
method in \S 3.4 (the intermediate boundary conditions adopt
$\V^{(1)}=\V^{(2)}=\V_\Gamma|_{t=(n+\frac 12)\Delta t},~~
\V^{(3)}=\V_\Gamma|_{t=(n+1)\Delta t}$) on $\Gamma$,
the iterative pressure Poisson equation method in \S 2,
the non--upwind  staggered mesh compact scheme in \S 3.1, $Re=10^4$,
the initial condition: \V={\bf 0} on inner grids, $\V=\V_\Gamma$,
$N=256\times 256, \Delta t=0.6\Delta x$,
set $\epsilon=10^{-5}$ for (2.8).
Figure 3 presents the stream fields for time steps $n=73840,86800,121840$.

To verify the {\bf equivalency} of the iterative method (2.6)(2.7) and the original
equations (2.3)(2.4), from $t=80640\times 0.6/256$,
we alter $\epsilon$ 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 $\epsilon=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\times 0.6/256$, we alter $\epsilon$ in (2.8)
to $10^{-13}, Re$ to $10^3$;

{(2)}set N=64,$Re=10^4, \epsilon=10^{-14}$, calculate from the beginning;

{(3)}set N=64,$Re=10^3, \epsilon=10^{-13}$, calculate from the beginning.

\vfill

{\small
$\begin{array}{ll}
	\begin{array}{l}\hbox{Fig.1a. Streamlines, $128\times 128$ grids},
		Re=5000,\\ \hbox{the upwind compact difference scheme}
	\end{array} &
	\begin{array}{l} \hbox{Fig.1b. velocity $u$ profile,}x=0.5 \\
		\hbox{----: our results, $\cdots$: results in [14]}
	\end{array}
\end{array}$
\eject

$\begin{array}{ll}
	&\\&\\&\\&\\&\\&\\&\\&\\&\\&\\&\\&\\&\\&\\&\\&\\&\\&\\&\\&\\&\\&\\&\\&\\
	\begin{array}{l}\hbox{Fig.1c.  velocity $v$ profile,}y=0.5 \\
		\hbox{----: our results,}\\ \cdots: \hbox{results in [14]}
	\end{array} &
	\begin{array}{l} \hbox{Fig.2a1. Streamlines, $256\times 256$ grids,}\\
		\hbox{the compact difference scheme},\\ Re=10000,t=t_0+12.7490
	\end{array}
\end{array}$
\vfill
$\begin{array}{ll}
	\begin{array}{l}\hbox{Fig.2a2. Streamlines, $256\times 256$ grids,}\\
		\hbox{the compact difference scheme},\\ Re=10000,t=t_0+13.9111
	\end{array} &
	\begin{array}{l} \hbox{Fig.2b1. velocity $u$ profile,}x=0.5 \\ \quad \\ \quad
	\end{array}
\end{array}$
\eject

$\begin{array}{ll}
	&\\&\\&\\&\\&\\&\\&\\&\\&\\&\\&\\&\\&\\&\\&\\&\\&\\&\\&\\&\\&\\
	\begin{array}{l}\hbox{Fig.2c1. velocity $v$ profile,}y=0.5 \\ \quad \\ \quad \\ \quad
	\end{array} &
	\begin{array}{l} \hbox{Fig.3a. Streamlines, $256\times 256$ grids,}\\ 
		Re=10000, t=73840\times 0.6/256 \\ 
		\hbox{values of stream function:}-0.002,-0.001, \\ 
		-0.0003,0,0.0003,0.01,0.02,\cdots,0.09
	\end{array}
\end{array}$
\vfill
$\begin{array}{ll}
	\begin{array}{l}\hbox{Fig.3b. Streamlines, $256\times 256$ grids,}\\
		Re=10000, t=86800\times 0.6/256 \\
		\hbox{values of stream function:}-0.002,-0.001, \\ -0.0003,
		0,0.0003,0.01,0.02,\cdots,0.09,0.098
	\end{array} &
	\begin{array}{l} \hbox{Fig.3c. Streamlines, $256\times 256$ grids,}\\
		Re=10000, t=121840\times 0.6/256 \\
		\hbox{stream function values:}-0.002,-0.001, \\ -0.0003,
		0,0.0003,0.01,0.02,\cdots,0.09,0.1
	\end{array}
\end{array}$
}

{\small
Further discusses,proofs and deductions will be put in our website:

http://www.cerse.psu.edu/yu/ipp/}

{\small
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}
\eject

\vskip1.5em \begin{center}{\bf 
			References		}\end{center}

\parindent=0pt

\def\toto#1#2{\centerline{\hbox to0.4cm{#1\hss}\parbox[t]{13.5cm}{#2}}}
\def\totot#1#2{\centerline{\hbox to0.6cm{#1\hss}\parbox[t]{13.5cm}{#2}}}

\toto{[1]}{ LIU Hong, { Study on high order accurate difference schemes and direct numerical
	simulations of the unsteady vortex}, {\it Master thesis of the Institute of Mechanics,
	Chinese Academy of Sciences}, Beijing, 1992, (main results are published in [19])}
\toto{[2]}{ P. J. Roache, Computational Fluid Dynamics, Hermosa Publishers, Albuquerque,1972}
\toto{[3]}{ A. J. Chorin, { A numerical method for solving incompressible viscous
	 flow problems}, {\it J. Comput. Phys.}, {\bf 2}: 1, 1967, 12--26}
\toto{[4]}{ F. H. Harlow, J. E. Welch,{ Numerical calculation of time--dependent viscous
	 incompressible flow of fluid with free surface}, 
	{\it Phys. Fluids}, {\bf 8}:12, 1965, 2182--2189}
\toto{[5]}{ YU Xin, { Divergence--free schemes for equations of incompressible fluid dynamics}, 
	{\it Doctoral dissertation of the Institute of Mechanics,
	Chinese Academy of Sciences}, Beijing, 1990}
\toto{[6]}{ C. Cuvelier, et al., Finite Element Methods and Navier--Stokes Equations, 
	D. Reidel Publishing Company, 1986}
\toto{[7]}{ YU Xin, A dimensional reduction method for incompressible fluid dynamics:
      I. The essential algorithm,
	{\it Chinese Journal of Computational Physics}, {\bf 2}:3, 1985, 337--346; 
     II. The basis of $V^h$ and the error estimates, 
	{\it Chinese Journal of Computational Physics}, {\bf 3}:2, 1986, 217--226;
    III. Second order accurate finite element schemes, 
	{\it Chinese Journal of Computational Physics}, {\bf 5}:2, 1988, 211--220;
    IV. Three dimensional Navier--Stokes equations, 
	{\it Chinese Journal of Computational Physics}, {\bf 6}:1, 1989, 104--116}
\toto{[8]}{ LIU Hong, FU Dexun, MA Yanwen, { Hopf bifurcation of the driven flow 
	in a square cavity}, {\it Theories, Methods and Applications of Computational
	Fluid Mechanics}, Science Press, Beijing, 1992, 267--270}
\toto{[9]}{ S. G. Rubin, R. A. Graves,Jr, { Viscous flow solutions with a cubic spline
	approximation}, {\it Computers and Fluids}, {\bf 3}:1, 1975, 1--36}
\totot{[10]}{ N. Baba, H. Miyata, H. Kajitani, {\it Journal of the Society of 
	Naval Architects of Japan}, {\bf 159}, 1987, 33}
\totot{[11]}{ S. Abdallah, { Numerical solutions for the incompressible Navier--Stokes
	equations in primitive variables using a non--staggered grid II},
	{\it J. Comput. Phys.}, {\bf 70}, 1987, 193--202}
\totot{[12]}{ FU Dexun, MA Yanwen, A Collection of Technical Papers, International
	Symposium on Computational Fluid Dynamics, 1989}
\totot{[13]}{ C. H. Bruneau, C. Jouron, { An efficient scheme for solving
	steady incompressible Navier--Stokes equations}, 
	{\it J. Comput. Phys.}, {\bf 89}, 1990, 389--413}
\totot{[14]}{ U. Ghia, K. N. Ghia, C. T. Shin, { High--Re solutions for incompressible
	flow using the Navier--Stokes equations and a multigrid method}, 
	{\it J. Comput. Phys.}, {\bf 48},(1982),387--411}
\totot{[15]}{ YU Xin, { Numerical solutions for the incompressible Navier--Stokes
	equations}, {\it Handbook for Engineering Applied Mechanics}, Chapter 7 of 
	{\it Computational Mechanics}, Shanxi People Press, Xi'an, 1996}
\totot{[16]}{ YU Xin, { A staggered mesh compact difference scheme and a 
	pressure--Poisson--equation that satisfies the equivalency},
	{\it Chinese J. of Numer. Math. and Appl.} {\bf 19}: 2, 1997, 73--81;
	(Original in Chinese: {\it Mathematica Numerica Sinica}, {\bf 19}:1, 1997, 83--90)}
\totot{[17]}{ R. Peyret, T. D. Taylor, Computational Methods for Fluid Flow,
	Springer--Verlag, 1983}
\totot{[18]}{ YU Xin, { Iterative pressure Poisson equation method for solving
			unsteady incompressible N--S equations},
	{\it Proc. of Beijing Symposium on Computational Fluid Mechanics},
	{\bf 5}, Beijing, 1993, 127--136}
\totot{[19]}{ LIU Hong, FU Dexun, MA Yanwen, { Upwind compact schemes and 
	direct numerical simulations of the driven flow in a square cavity},
	{\it Science in China}, Series A, {\bf 23}:6, 1993, 657--665}

\end{document}
