

%	 published in Mathematica Numerica Sinica, Vol.19,No.1, 1997, 83-90,

%	 and was republished in English:

%	    Chinese J. of Numer. Math. and Appl. , Vol. 19, No. 2, 1997



\documentstyle[12pt]{article}

\textwidth 16cm

\textheight 23cm

\def\pmb#1{\mbox{\boldmath $#1$}}

\def\V{\pmb V}

\def\A{\pmb A}



\begin{document}

\vskip-9em

{{\footnotesize

Chinese J. of Numerical Mathematics and Application, {\bf 19}:2, 1997, 73--81 (in English)



Mathematica Numerica Sinica, {\bf 19}:1, 1997, 83--90 (in Chinese)

}

\vskip6em

\begin{center}



{\bf 	A STAGGERED MESH COMPACT DIFFERENCE SCHEME AND A

PRESSURE--POISSON--EQUATION THAT SATISFIES THE EQUIVALENCY}\footnote{Received January 26, 1996.

\quad Copyright 1997 by Allerton Press, Inc.}



\vskip1em

                    YU Xin



(\,~\sl Institute of Mechanics, Chinese Academy of Sciences, Beijing 100080\,)



\end{center}



\vskip2em\begin{center}{\bf 	Abstract	}\end{center}\vskip1em



~~ (1) A staggered mesh compact difference scheme is

presented for solving the unsteady viscous incompressible Navier--Stokes

equations. It is fourth order accurate both in the spatial direction and

in the time direction,

at least third order accurate on the boundary; (2) Describe a

pressure--Poisson--equation that is equivalent to the discrete continuity

equation provided the discrete momentum equations remain. The discrete

continuity equation may have derivative boundary conditions, e.g., the compact

difference scheme; (3) A new ADI iterative method is proposed. The

pressure--Poisson--equation is in the discrete form. It is difficult to be

solved with a usual ADI method. We translate it to be tridiagonal

in each spatial direction of each step of the ADI iterations,

then add a pseudo time term to get a tridiagonal equation which is

easily to be solved; (4) The driven flow in a square cavity with

$Re=10000$ is simulated numerically.





\vskip2em



\centerline{\bf 1. Introduction}



\vskip1em



With the developing of the electronic computers,

more and more physical problems can be simulated numerically.

While there are still a lot of nonlinear problems

need too much CPU time and too large computer memory to calculate.

The improvement of the numerical methods can reduce the computer time and

memory greatly, e.g., for a two dimensional unsteady problem,

in order to reduce the error to $N^{-4}$,

a second order accurate scheme needs about $(N^2)^3$ spatial and time grids,

while a fourth order accurate scheme needs only $N^3$ grids.

The ratio of them is $N^3$. The ratio for the computer time is even more.

$N^3=4096$ when $N=16$, $N^3\approx 1.678\times 10^7$ when $N=256$.



The advantages of compact difference schemes over traditional methods include

the relatively high order of accuracy using a compact stencil, a better

(linear) stability, a better resolution for high frequency wave,and usually

fewer boundary points to handle.$^{[8,1,2,3,5,6,10]}$



With a three points stencil in each time level, the scheme in this paper

achieves fourth order accurate both in the spatial direction (\S 2.2) and in

the time direction (\S 2.4), at least third order accurate on the boundary

(\S 2.3). This paper is a revised version of [9].



\vskip3em

\begin{center}



{\bf 2. A Fourth Order Accurate Staggered Mesh Compact Difference Scheme}



\end{center}



The compact difference scheme in this paper improves in the

treatment of the pressure boundary differences, and the discretization of

the time derivative in [1][2].



\noindent

{\bf 2.1 The Unsteady Viscous Incompressible Navier--Stokes Equations}:

$$

\frac{\partial\V}{\partial t}+\A(\V)+\pmb\nabla p={\bf 0}, \quad

    {\rm in } \quad\Omega,\quad(\hbox{ where }

\A(\V)=(\V\cdot\pmb\nabla)\V-\nu\nabla^2\V~),              \eqno{(2.1)}

$$

$$

\hbox{div}~\V=0, \quad\quad {\rm in } \quad \Omega \quad\quad

\quad\quad\quad\quad\quad\quad\quad\quad\quad     \eqno{(2.2)}

$$

\noindent

{\bf 2.2 A Two Dimensional Staggered Mesh Compact Difference Scheme}



The compact difference scheme described in this section has the fourth order

accuracy for the spatial discretization of (2.1)(2.2), the first order accuracy

for the time derivative in (2.1). We consider the discrete form of (2.1)(2.2)

$$

\frac{\V^{n+1}-\V^{*}}{\Delta t}+\pmb\nabla_hp^{n+1}=0,

    \quad\hbox{ ( here }\V^*=\V^n-\Delta t\A_h(\V^n)~),     \eqno{(2.3)}

$$

$$

\hbox{div}_h~\V^{n+1}=0,\quad\quad\quad\quad\quad\quad\quad \eqno{(2.4)}

$$

\noindent

where $\V^n=(u^n,v^n),{\A}(\V)=(uu_x+vu_y-\nu(u_{xx}+u_{yy}),

uv_x+vv_y-\nu(v_{xx}+v_{yy})$),$\A_h,\pmb\nabla_h, $ div$_h$

are discrete forms of $\A,\pmb\nabla$, div respectively:



(1) $u_x,u_y$ in $\A_h(\V^n)$($v_x,v_y$ in $\A_h(\V^n)$ similar):



$$

\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},\eqno{(2.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},\eqno{(2.5)_2}

$$

\hskip1.5em

(2) $u_{xx},u_{yy}$ in $\A_h(\V^n)$($v_{xx},v_{yy}$ in $\A_h(\V^n)$ similar):

$$

\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}.                            \eqno{(2.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}.                            \eqno{(2.6)_2}

$$

\hskip1.5em (3) $p_x$ in $\pmb\nabla_hp$ ( $p_y$ similar):

$$

\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},                         \eqno{(2.7)}

$$

\hskip1.5em (4) $u_x$ in div$_h \V$ ( $v_y$ in div$_h \V$ similar):

$$

\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},                \eqno{(2.8)}

$$

\hskip1.5em (5) $u$ in $uu_x$ in ( $\V\cdot\pmb\nabla$ )\V: $u_{i,j+\frac 12}$,

$v$ in $vv_y: v_{i+\frac 12,j}$, $u$ in $uv_x$:

$$

u_{i+\frac{1}{2}, j} = \frac{u_{i, j-\frac{1}{2}} +

u_{i+1, j-\frac{1}{2}} + u_{i, j+\frac{1}{2}} + u_{i+1, j+\frac{1}{2}}}{4}

-\frac{1}{8}[(\Delta x)^2u_{xx} + (\Delta

y)^2u_{yy}]_{i+\frac{1}{2},j}, 				\eqno{(2.9)}

$$

where $u_{xx},u_{yy}$ can be get from the results of (2.6),

$v$ in $vu_y$ similar,



(6) for $\V_t$: (2.3) uses the first order accurate difference. A fourth

order accurate Runge--Kutta algorithm is given in \S 2.4.



\noindent

{\bf 2.3~~Difference Formulae on the Boundary}



Set $\Omega=\{ (x,y){\big |} 0<x<1,0<y<1\}, \Delta x=\Delta y=1/N$.

Consider the boundary condition  $\V|_\Gamma=\V_\Gamma=(u^\Gamma,v^\Gamma)$,

where $\Gamma=\partial \Omega$.

We describe the difference formulae on the left edge $x=0$.



(1)$_1$ $u_x|_{x=0}=-v^\Gamma_y$,\quad (since (2.2)),

    \quad (for $v:\quad v_y|_{y=0}=-u^\Gamma_x) \hfill  (2.10)$



(1)$_2$ $v_x$ at $x=\frac {3}{4}\Delta x$ and

    $u_y$ at $y=\frac {3}{4}\Delta y$:

$$

\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},\quad\quad

 \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},    \eqno{(2.11)}

$$

\hskip1.5em (2)$_1$ $u_{xx}$ at $x=\frac{10}{9}\Delta x$

    ( $v_{yy}$ at $y=\frac{10}{9}\Delta y$ similar):

$$

\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},      \eqno{(2.12)}

$$

where $u'_{0,j+\frac {1}{2}}=-(v^\Gamma_y)_{0,j+\frac {1}{2}}$, cf (2.10),



(2)$_2$ $v_{xx}$ at $x=\frac {2}{3}\Delta x$

    ( $u_{yy}$ at $y=\frac {2}{3}\Delta y$ similar):

$$

\frac{5v^{\prime\prime}_{ \frac{1}{2},j} + v^{\prime\prime}

_{\frac{3}{2},j}}{6} =\frac{4}{3}\cdot\frac{2v_{0,j}-3v_{\frac{1}{2},j}

+v_{ \frac{3}{2},j}}{\Delta x^2},           \eqno{(2.13)}

$$

\hskip1.5em (3) $p_x$ at $x=\Delta x$ ( $p_y$ at $y=\Delta y$ similar):

$$

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}, \eqno{(2.14)}

$$

(2.14) improves the only second order accurate difference formula in [1] [2]

to third order accurate,



(4) $u_x$ in div\V \quad at $x=\frac 12\Delta x$ ( $v_y$ at

    $y=\frac {1}{2}\Delta y$ similar):

$$

\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},

\hbox{(where $u'_{0,j+\frac {1}{2}}

=-(v^\Gamma_y)_{0,j+\frac{1}{2}}$ )},          \eqno{(2.15)}

$$

\hskip1.5em (5) boundary values of $u_{xx},u_{yy}$ in (2.9):

$$

u_{xx}|_{x=\frac{1}{2}\Delta x}=\frac{u_x|_{x=\Delta x}-u_x|_{x=0}}{\Delta x},

\quad\quad  u_{yy}|_{x=0}=u^\Gamma_{yy}|_{x=0},    \eqno{(2.16)}

$$

where $u_x|_{x=\Delta x}$ can be get from the results of (2.5),

 $u_x|_{x=0}=-v^\Gamma_y|_{x=0}$ since (2.2).



\noindent

{\bf 2.4~~The Fourth Order Accurate Runge--Kutta Method for the Time Derivative}



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$. To solve $\V_t+\pmb f(\V)={\bf 0}$,

we use the fourth order accurate Runge--Kutta method:

$$

\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{(2.17)}

$$

\noindent

where $\pmb k_1=\pmb f(\V^n),~\pmb k_2=\pmb f(\V^n-\frac {\Delta t}{2}\pmb k_1),

 \pmb k_3=\pmb f(\V^n-\frac {\Delta t}{2}\pmb k_2),

 \pmb k_4=\pmb f(\V^n-\Delta t\pmb k_3)$.



\vskip3em



\begin{center}

{\bf		 3. A Pressure--Poisson--Equation that Satisfies the Equivalency}

$^{[1]}$

\end{center}



\vskip1em



Let us consider the following discrete scheme,(e.g., the compact difference

scheme in \S 2, and the second order accurate MAC scheme$^{[4]}$),

$$

\frac{\V^{n+1}-\V^{*}}{\Delta t}+\pmb\nabla_hp^{n+1}={\bf 0},

    \quad\hbox{ ( where }\V^*=\V^n-\Delta t\A_h(\V^n)~),     \eqno{(3.1)}

$$

$$

\hbox{div}_h~\V^{n+1}=0.\quad\quad\quad

            \quad\quad\quad\quad                            \eqno{(3.2)}

$$

$$

\hbox{ Define}\hskip6em    \pmb{NS}_h

=\frac{\V^{n+1}-\V^{*}}{\Delta t}+\pmb\nabla_hp^{n+1},\hskip10em \eqno{(3.1)_A}

$$

then (3.1) can be written to $\pmb{NS}_h$={\bf 0}.



To avoid div$_h(\pmb{NS}_h)$ invalid, suppose (3.1) is defined on the

difference mesh of \V, e.g., for the staggered mesh scheme in \S 2:

$$

\pmb{NS}_h \sim((NS_h)^1_{i,j+\frac{1}{2}},(NS_h)^2_{i+\frac{1}{2},j}),

\quad\quad\V^n_h\sim(u^n_{i,j+\frac{1}{2}},v^n_{i+\frac{1}{2},j}).

$$

We write div$_h\V$ in the form containing the boundary conditions:

$D_h(\V,\V_\Gamma,\pmb\nabla\V_\Gamma)$,

where \V\, is defined on the inner difference points,

$\V_\Gamma=\V|_\Gamma,\pmb\nabla \V_\Gamma=\pmb\nabla\V|_\Gamma$.

(The compact scheme in \S 2 does not use all components of \V.

It uses only the components which can be get from $\V_\Gamma$ and

the continuity equation:

$\pmb\nabla(\V_\Gamma\cdot \pmb n)$

and $\pmb\nabla(\V_\Gamma\cdot\pmb \tau)\cdot\pmb \tau$, here $\pmb n$

and $\pmb \tau$ are the unit normal and tangent vectors on $\Gamma$.

the MAC scheme $^{[4]}$ does not use $\pmb\nabla\V_\Gamma$, i.e.,

there is no derivative boundary condition in the continuity equation).

Thus the discrete continuity equation (3.2) can be written to

$$

D_h(\V^{n+1},\V^{n+1}_\Gamma,\pmb\nabla\V^{n+1}_\Gamma)=0.    \eqno{(3.2)'}

$$

From (3.1) (3.1)$_A$,\quad hence \quad

$\pmb{NS}_h={\bf 0},\V^{n+1}-\alpha \pmb{NS}_h=\V^{n+1}$,

(where $\alpha$ is an arbitrary constant). Therefore from $(3.2)'$ we get

$$

D_h(\V^{n+1}-\alpha\pmb{NS}_h,

\V^{n+1}_\Gamma,\pmb\nabla\V^{n+1}_\Gamma)=0.               \eqno{(3.3)}

$$

{\bf Theorem 3.1} For any arbitrary constant $\alpha$, (3.3)(3.1) are

equivalent to (3.2)(3.1).



Set $\alpha=\Delta t$.

Then (3.3) with (3.1) leads to the pressure--Poisson--equation

$$

D_h(\pmb\nabla _h p^{n+1},{\bf 0},{\bf 0})

=\frac{1}{\Delta t}D_h(\V^*,\V^{n+1}_\Gamma

,\pmb\nabla\V^{n+1}_\Gamma)                   \eqno{(3.4)}

$$

\noindent

provided $D_h(\V,\V_\Gamma,\pmb\nabla\V_\Gamma)$ is a linear function of

$(\V,\V_\Gamma,\pmb\nabla\V_\Gamma)$. RHS of (3.4) is\\

$\frac {1}{\Delta t}{\rm div}_h \V^*$ (cf (2.8)(2.15)) with

$\V^*|_\Gamma=\V^{n+1}_\Gamma,(\pmb\nabla\V^*)_\Gamma=\pmb\nabla\V^{n+1}_\Gamma$

(see  $u'_{0,j+\frac{1}{2}}$ in (2.15)).\\

LHS of (3.4) is div$_h(\pmb\nabla_hp^{n+1})$ (see (2.7)(2.14) (2.8)(2.15))

with $(\pmb\nabla p^{n+1})_\Gamma={\bf 0}$,

(see $u_{0,j+\frac {1}{2}}$ in (2.15)),

( $\pmb\nabla(\pmb\nabla p^{n+1}))_\Gamma={\bf 0}$,

(see $u'_{0,j+\frac {1}{2}}$ in (2.15)). These {\bf numerical boundary\\

 conditions} do not affect the result of

$p^{n+1}$ by theorem 3.1. See [4] \S 6.3.1.



So we can solve the pressure--Poisson--equation (3.4) for $p^{n+1}$,

i.e., use (3.4)(3.1) or (3.3)(3.1) instead of (3.1)(3.2)

to solve the Navier--Stokes equations (2.1)(2.2).



\vskip3em

\begin{center}



{\bf   4. A New ADI Iterative Method for Solving the Pressure--Poisson--Equation}



\end{center}



It is difficult to solve the discrete pressure--Poisson--equation (3.4) with

a usual ADI method. We translate it to be tridiagonal in each spatial direction

of each step of the ADI iterations ((4.5)--(4.7)), then add a pseudo time

term to get a tridiagonal equation which is easily to be solved (see (4.7)),

(a usual ADI method adds the pseudo time term at the beginning).

(3.4) can be written to

$$

{\rm div}_h(\pmb\nabla_h p)=\frac {1}{\Delta t}{\rm div}_h \V^*,

        \quad\quad({\hbox{+ boundary condition}})               \eqno{(4.1)}

$$

The new ADI method in the $x$ direction (similar in the $y$ direction):



\noindent

\hskip1.5em 

(1) calculate $p_y$ in $\pmb\nabla_h p$ using the initial $p$ or the

last step $p$, and the formulae corresponding to (2.7) (2.14):

$$

\displaylines{\hspace {0.5cm}

\frac {(p_y)_{i+\frac {1}{2},j-1}+22(p_y)_{i+\frac

{1}{2},j}+(p_y)_{i+\frac {1}{2},j+1}}{24}=\frac {p_{i+\frac {1}{2},j+\frac

{1}{2}}-p_{i+\frac {1}{2},j-\frac {1}{2}}}{\Delta y},(j=2,3,\cdots),

\hfill (4.2)_1

\cr\noalign{\vskip2mm}

\hspace {0.5cm}

(p_y)_{i+\frac {1}{2},1}+\frac {(p_y)_{i+\frac {1}{2},1}-2(p_y)_{i+\frac

{1}{2},2}+(p_y)_{i+\frac {1}{2},3}}{24}=\frac {p_{i+\frac {1}{2},\frac

{3}{2}}-p_{i+\frac {1}{2},\frac {1}{2}}}{\Delta y}, \hfill (4.2)_2\cr}

$$

\hskip1.5em (2) calculate $p_{yy}$ in div$_h(\pmb\nabla_h p)$ using

the formulae corresponding to (2.8) (2.15):

$$

\displaylines{\hspace {0.5cm}

\frac {(p_{yy})_{i+\frac {1}{2},j-\frac {1}{2}}+22(p_{yy})_{i+\frac

{1}{2},j+\frac {1}{2}}+(p_{yy})_{i+\frac {1}{2},j+\frac {3}{2}}}{24}=

\frac {(p_y)_{i+\frac {1}{2},j+1}-(p_y)_{i+\frac {1}{2},j}}{\Delta y},

(j=1,2,\cdots)

\cr\noalign{\vskip2mm}

\hspace {0.5cm}

\frac {2(p_{yy})_{i+\frac {1}{2},0}+15(p_{yy})_{i+\frac {1}{2},\frac

{1}{2}}+(p_{yy})_{i+\frac {1}{2},\frac {3}{2}}}{18}=\frac

{(p_y)_{i+\frac {1}{2},1}-(p_y)_{i+\frac {1}{2},0}}{\Delta y},

\hfill{\hbox {\raise3ex\hbox{ (4.3)}}}

\cr\noalign{\vskip2mm}\hfill (\,(p_y)_{i+\frac {1}{2},0}=0,

(p_{yy})_{i+\frac {1}{2},0}=0\,)\hspace {5cm}\cr}

$$

\hskip1.5em

(3)\quad\quad $ p_{xx}=\frac {1}{\Delta t} {\rm div}_h \V^*-p_{yy} \hfill (4.4)$



(4) solve the equations corresponding to (2.8) (2.15)

for $p_x$ with the recursion



--correction method:

$$

 \displaylines{\hspace {0.5cm}

\frac {(p_{xx})_{i-\frac {1}{2},j+\frac {1}{2}}+22(p_{xx})_{i+\frac

{1}{2},j+\frac {1}{2}}+(p_{xx})_{i+\frac {3}{2},j+\frac {1}{2}}}{24}=

\frac {(p_x)_{i+1,j+\frac {1}{2}}-(p_x)_{i,j+\frac {1}{2}}}{\Delta x},

(j=1,2,\cdots )

\cr\noalign{\vskip2mm}

\hspace {0.5cm}

\frac {2(p_{xx})_{0,j+\frac {1}{2}}+15(p_{xx})_{\frac {1}{2},j+\frac

{1}{2}}+(p_{xx})_{\frac {3}{2},j+\frac {1}{2}}}{18}=\frac

{(p_x)_{1,j+\frac {1}{2}}-(p_x)_{0,j+\frac {1}{2}}}{\Delta x},

\hfill{\hbox {\raise3ex\hbox{ (4.5)}}}

\cr\noalign{\vskip2mm}\hfill (\,(p_x)_{0,j+\frac {1}{2}}=0,

(p_{xx})_{0,j+\frac {1}{2}}=0\,)\hspace {5cm}\cr}

$$

\hskip1.5em (5) calculate the LHS of (2.7) and (2.14):

$$

 \displaylines{\hspace {0.5cm}

(p^*_x)_{i,j+\frac {1}{2}}=\frac {(p_x)_{i-1,j+\frac

{1}{2}}+22(p_x)_{i,j+\frac {1}{2}}+(p_x)_{i+1,j+\frac{3}{2}}}{24},

(i=2,3,4,\cdots) \hfill

{\hbox {\raise-4ex\hbox{ (4.6)}}}

\cr\noalign{\vskip0.1mm}

\hspace {0.5cm}

(p^*_x)_{1,j+\frac {1}{2}}=\frac {25(p_x)_{1,j+\frac

{1}{2}}-2(p_x)_{2,j+\frac {1}{2}}+(p_x)_{3,j+\frac{1}{2}}}{24}, \hfill\cr}

$$

\hskip1.5em

(6) differencing (2.7) (2.14), adding a pseudo time term, the resulting

equation on $p$ is tridiagonal as follows, solve it to get $p$

$$

 \displaylines{\hspace {0.5cm}

\frac {{(p^*_x)}_{i+1,j+\frac {1}{2}}-(p^*_x)_{i,j+\frac {1}{2}}}{\Delta

x}=\frac {p_{i-\frac {1}{2},j+\frac {1}{2}}-2p_{i+\frac {1}{2},j+\frac

{1}{2}}+p_{i+\frac {3}{2},j+\frac {1}{2}}}{(\Delta x)^2}-\frac

{p_{i+\frac {1}{2},j+\frac {1}{2}}}{\Delta t_*},(i=1,2\cdots)

\cr\noalign{\vskip2mm}

\hspace {0.5cm}

\frac {{(p^*_x)}_{1,j+\frac {1}{2}}}{\Delta x}=\frac {p_{\frac

{3}{2},j+\frac {1}{2}}-p_{\frac {1}{2},j+\frac {1}{2}}}{(\Delta

x)^2}-\frac {p_{\frac {1}{2},j+\frac {1}{2}}}{\Delta t_{*1}}, \hfill

{\hbox {\raise3ex\hbox{ (4.7)}}}\cr}

$$

The corresponding formulae on the right and the upper boundaries can be

obtained by translating $(4.2)_2, (4.3)_2, (4.5)_2, (4.6)_2, (4.7)_2$ by:



\centerline{$i\rightarrow N-i,\Delta x\rightarrow -\Delta x$ and

$j\rightarrow N-j,\Delta y\rightarrow -\Delta y$.}



\noindent

The computations in this paper and in [1] were performed by the above

ADI method.



\vskip3em



\centerline{\bf   5. Computations of the Driven Flow in a Square Cavity }



\vskip1em



Consider the Navier--Stokes equations (2.1) (2.2) on $0<x<1,\, 0<y<1$

with boundary conditions $\V=(-1,0)$, (on $y$=1), \V=(0,0),

(on $y=0$, $x=0$, $x=1$), $Re=1/\nu=$10000.

We simulate the driven flow in a square cavity with the staggered mesh

compact difference scheme with $\Delta x=\Delta y=1/N$, where $N=256$.



To verify the equivalency (theorem 3.1) further, in [1] we used a coarse

grid ($N=8$), ADI iterations, reduced the error to $|{\rm div}_h\V|<10^{-15}$

in the calculation.

Here we adopt a multigrid method, $N=256, Re=10000$,

reduce the error to $|{\rm div}_h\V|<6\times 10^{-14}$,

This shows the pressure--Poisson--equation (4.1) (or (3.4)) can replace

the continuity equation div$_h \V=0$

(suppose retaining the momentum equation (3.1)).



The streamlines (in figure 1) are contours of the stream function $\psi$.

In order to obtain $\psi$, here we solve a Poisson equation deduced from

$$

\frac{\psi_{i,j+1}-\psi_{i,j}}{\Delta y}

=u_{i,j+\frac{1}{2}}+\frac{(\Delta y)^2}{24}u_{yy},\quad

\frac{\psi_{i+1,j}-\psi_{i,j}}{\Delta x}

=-(v_{i+\frac{1}{2},j}+\frac{(\Delta x)^2}{24}v_{xx}),	\eqno{(5.1)}

$$

\noindent

where $u_{yy},v_{xx}$ are at least second order accurate,

(then (5.2) with (5.1) is a fourth order accurate approximation of

the continuity equation (2.2),

$$

[(\psi_{i+1,j+1}-\psi_{i+1,j})-(\psi_{i,j+1}-\psi_{i,j})]

-[(\psi_{i+1,j+1}-\psi_{i,j+1})-(\psi_{i+1,j}-\psi_{i,j})]=0.	\eqno{(5.2)}

$$

As a result, (2.3)--(2.7),(5.1),(5.2),(2.9)--(2.14),(2.16) is another fourth

order accurate compact difference scheme for solving the Navier--Stokes

equations (2.1)(2.2)).



$\| \cdot\|$ in figure 3 is defined by

$$

\|\V\|=\frac{1}{2N}\hbox{\huge (}\,\sqrt{\sum_{i,j}(u_{i,j+\frac{1}{2}})^2}

+\sqrt{\sum_{i,j}(v_{i+\frac{1}{2},j})^2}\,\hbox{\huge )},

\quad\quad	\hbox{\large (}\V=(u,v)\hbox{\large )},		\eqno{(5.3)}

$$

$$

\V_4=\V|_{t=t_1+536.655},\quad \V_5=\V|_{t=t_1+542.340}.\quad\quad\quad\quad\quad\quad\quad\quad

$$



The computing time of the results in [2] was relatively short. Here we

employ the multigrid method for computation of the pressure.

Longer calculation shows that the time period of the numerical solution is

about $T=6.727$, approximating $T=6.36$ (the period given in [3]).

It can be seen that the numerical solution is still not completely periodic.

In figure 3 near $t=t_1+549$

$$

\min\limits_n \|(\V|_{t=t_1+n\Delta t})-\V_5 \|

= \|(\V|_{t=t_1+n_1\Delta t})-\V_5 \|  \approx 1.13 \times 10^{-5}.

$$

With careful calculations we can get

$$

\min\limits_{0\leq n\leq 100} \|(\V|_{t=t_1+n_1\Delta t + n\Delta t /100})

-\V_5 \|\approx 7.15 \times 10^{-6}.

$$



The results in [3] demonstrates that the period shown by the figures

of the streamlines is a quarter of that shown by the figure of the velocity

at the geometric center. [7] gave an analogous result.

This coincides with the results in figure 3 and figure 2 here.

From figure 2, $u(\frac 12,\frac 12,t)$ and $v(\frac 12,\frac 12,t)$

can be approximated by

$$

|A(t)|\sin(2\pi t/T+\theta_1)

+|B(t)|\sin(2\pi t/(T/4)+\theta_2)+C(t),	\eqno{(5.4)}

$$

where $T$ also changes with time $t$,

$|B|$ becomes positive in figure 2a,\quad

$|B| \approx |A|$ in figure 2b,\quad $|B| \gg |A|$ in figure 2c.



\eject

\centerline{\bf		 References	}\vskip1em



\small

\def\hang{\hangindent\parindent}

\def\textindent#1{\indent\llap{#1\enspace}\ignorespaces}

\def\re{\par\hang\textindent}



\re{[1]} YU Xin, An Equivalent Pressure--Poisson--Equation for

N--S Equations and Staggered Mesh Compact Difference Schemes,\quad

{\sl First Asian Computational Fluid Dynamics Conference},



1995, Hong Kong, 937--942.



\re{[2]} YU Xin, An Iterative--Pressure--Poisson--Method for Solving

Unsteady Incompressible N--S Equations,\quad

{\sl Papers Presented at the Beijing Workshop on Computational Fluid Dynamics},

{\bf 5}, 1993, 127--136.



\re{[3]} LIU Hong, FU Dexun, MA Yanwen, Upwind Compact Method and

Direct Numerical Simulation of Driven Flow in a Sqare Cavity,

{\sl Science in China}, (Series {\bf A}), {\bf 23}, 6,1993, 657--665.



\re{[4]} R. PEYRET, T. D. TAYLOR,

Computational Methods for Fluid Flow, Springer--Verlag, 1983.



\re{[5]} S. G. RUBIN,\quad R. A. GRAVES, Jr.,\quad Viscous Flow Solutions

with a Cubic Spline



Approximation,\quad {\sl Computers and Fluids},{\bf 3}, 1975, 1--36.



\re{[6]} M. CIMENT,\quad S. H. LEVENTHAL,\quad Higher Order Compact Implicit

Schemes for the Wave Equation,\quad {\sl Mathematics of Computation},

{\bf 29}, 1975, 985--994.



\re{[7]} Jie SHEN, {\sl Journal of Computational Physics}, {\bf 95},

1991, 228--.



\re{[8]} Bernardo COCKBURN, Chi-Wang SHU, Nonlinearly Stable Compact Schemes

for Shock



Calculations, {\sl SIAM J. Numer. Anal.}, {\bf 31}, 3, 1994, 607--627



\re{[9]} YU Xin, A Staggered Mesh Compact Difference Scheme and a

Pressure--Poisson--Equation that Satisfies the Equivalency,

{\sl Papers Presented at the Beijing Workshop on Computational Fluid Dynamics},

{\bf 7}, 1995, 137--145.



\re{[10]} FU Dexun, MA Yanwen, LIU Hong, Upwind Compact Schemes and

Applications, \quad\sl Proc. of the 5th Int. Symp. on Computational Fluid

Dynamics\rm, {\bf I}, Sendai, JSCFD, 1993, 184--190.



\vskip2em

%(This paper was published in

% Mathematica Numerica Sinica, Vol.19,No.1, 1997, 83-90,(in Chinese),

% and published in:

%  Chinese J. of Numer. Math. and Appl. , Vol. 19, No. 2, 1997, (in English)



% if you want to view Figures, you can also email me: yu@cerse.psu.edu,

%I will mail them to you)

\end{document}

