Citation
Multilevel adaptive methods in computational fluid dynamics

Material Information

Title:
Multilevel adaptive methods in computational fluid dynamics
Creator:
Liu, Chaoqun
Publication Date:
Language:
English
Physical Description:
93 leaves : illustrations ; 29 cm

Subjects

Subjects / Keywords:
Fluid dynamics ( lcsh )
Multigrid methods ( lcsh )
Poisson's equation ( lcsh )
Navier-Stokes equations ( lcsh )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

Notes

Bibliography:
Includes bibliographical references.
General Note:
Submitted in partial fulfillment of the requirements for the degree, Doctor of Philosophy, Department of Mathematical and Statistical Sciences.
Statement of Responsibility:
by Chaoqun Liu.

Record Information

Source Institution:
|University of Colorado Denver
Holding Location:
Auraria Library
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
24655102 ( OCLC )
ocm24655102
Classification:
LD1190.L622 1989d .L58 ( lcc )

Full Text
Multilevel Adaptive Methods
in
Computational Fluid Dynamics
by
Chaoqun Liu
Computational Mathematics Group
Campus Box 170
University of Colorado at Denver
Denver, CO 80204
November,1989


MULTILEVEL ADAPTIVE METHODS
IN
COMPUTATIONAL FLUID DYNAMICS
by
Chaoqun Liu
B. A., Tsinghua University, Beijing, China, 1967
M. A., Tsinghua University, Beijing, China, 1981
A thesis submitted to the
Faculty of the Graduate School of the
University, of Colorado in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
Department of Mathematics
1989


This thesis for the Doctor of Philosophy degree
by
Chaoqun Liu
has been approved for the
Department of Mathematics
Steve McCormick
Williams L. Briggs
Karl Gustafson
Jan Mandel
)
Roland Sweet
Date November 27, 1989
)


Chaoqun Liu (Ph.D., Applied Mathematics)
Multilevel Adaptive Methods in Computational Fluid Dynamics
Thesis directed by Professor Steve McCormick
We first take Poisson's equation as an example to illustrate the basic
principles and schemes for the methods of multigrid (MG), fast adaptive
composite grid (FAC) and finite volume element (FVE). The computational
results show that local refinement can greatly enhance discretization
accuracy. Then we apply full multigrid (FMG) with cubic interpolation to
solve the elliptic grid generation (EGG) equations and use FAC to generate
a smooth patch grid. We also apply MG and FAC to solve a full potential flow
around the NACA0012 airfoil. For this application, discretizations are used
in the global grid that are different than in the patches. Computational
results show that MG can substantially accelerate convergence and FAC
local refinement can greatly improve accuracy, especially for problems
with shocks. We apply MG and FAC to driven cavity flow (incompressible
Navier-Stokes equations). Reviewing the weakness of classical finite
difference theory and conventional notions on artificial viscosity, we
introduce the rotated hybrid scheme and FVE scheme for medium and high
Reynolds number flows. We establish an equivalence between spatial
relaxation and time-marching schemes, then use this interpretation as a
basis for the development of an MG7 solver for pseudo-unsteady flows.
Finally, we introduce the source term correction (STC) method for
incompressible Navier-Stokes equations.
The ;form and content of this abstii
publication.
Sign
nd its


Dedicated to Mrs. Han, my grandmother
I


Acknowledgement
I am deeply in debt to Professor S. McCormick, the
founder of FAC, for his ideas, support, countless discus-
sions and helpful guidance.
Much of my work wa6 inspired by the dramatic develop-
ments in the field of multigrid, I am therefore thankful
to Professor A. Brandt for his pioneering work in the field
as well as his lectures and helpful discussions.
I am also indebted to the United States Air Force
Office of Scientific Research and the National Science
Fundation for their financial support.


CONTENTS
CHAPTER
I. INTRODUCTION
Motivation.........;........................... 1
Arrangement of the thesis...................... 2
II. POISSONS EQUATION
Equation and boundary conditions............... 4
FVE discretization........................... 5
Multigrid algorithm (FAS).................... 8
FVE accuracy.................................. 11
( Nested iteration (FMG).................... 11
FAC for nonuniform (composite) grids..... 12
III. ELLIPTIC GRID GENERATION (EGG)
EGG equations................................. 19
Finite difference equations.................... 21
Initial guess.................................. 21
Full approximation scheme (FAS)................ 22
FAC for local refinement....................... 23
Discretization error estimates..........}. 24
IV. TRANSONIC POTENTIAL FLOWS
Governing equations.......................... 26
Boundary conditions............................ 27
Discretization ............................... 27
Stability and artificial density............... 35
Multigrid for potential flows................. 36
Computational results......................... 38
V. THE ROTATED HYBRID SCHEME FOR PLANAR CAVITY
FLOWS
Differential equations and boundary
conditions.................................. 43
Conventional difference schemes for


11
the vorticity equation................. 43
Improved difference schemes for 1-D
^ vorticity equation..................... 45
Artificial viscosity revisited........47
Rotated hybrid scheme for the 2-D
vorticity equation......................... 50
Collective relaxation........................ 53
Multigrid................................... 55
Alternating direction block-line
Gauss-Seidel relaxation.................... 56
FAC for local refinement..................... 57
Computational results...............'.... 57
VI. FVE FOR DRIVEN CAVITY FLOW
Equations and boundary conditions............ 59
Difference scheme........................... 59
, Multigrid algorithm............................ 63
c FVE accuracy.............................:.. 64
Nested iteration..................^...... 65
FAC..................'....................... 65
Medium and high Re flow...................... 68
VII. MG FOR PSEUDO TIME-DEPENDENT FLOWS
)
Explicit time-marching and Jacobi
iteration.................................. 70
Explicit time^marching and Gauss-Seidel
iteration.................................. 71
ADI time-marching and line relaxation.... 72
Generalization............................. 73
MG as an explicit Euler solver............... 74
VIII. SOURCE TERM CORRECTION (STC) FOR
' INCOMPRESSIBLE NAVIER-STOKES EQUATIONS
Reformulating the governing equations
and boundary conditions.................... 79
Discretization................................ 81
The STC method................................ 82


Ill
Iteration scheme and STC algorithm....... 84
Solution for pressure......................... 85
Other flow problems.............'........ 86
Preliminary computational results........... 88
REFERENCES............................................. 90

T.
(
I


CHAPTER I
INTRODUCTION
1. Motivation^
Great progress has been achieved during the past two
decades in both computational fluid dynamics (CFD) and
applied mathematics. Some of the achievements include the
development of elliptic grid generation (EGG) [1],
multigrid [2-3], and algorithms and difference schemes
for solving Euler and Navier-Stokes equations [4-10].
However, there are still several barriers to the efficient
and accurate solution of many large scale flow problems,
including transonic, high Reynolds number, separated and
- turbulent flows. Even with5the current and foreseeable
computing machinery, substantial progress must-be made in
algorithms in order for the numerical solution of these
problems to be practically feasible.
It is well known that MG is an optimal order process
for solving many partial differential equations (PDEs). In
fact, MG schemes have proved effective for a wide number of
problems which contribute to application of MG to a variety
of flow problems [11-18]. However, there are still many
impediments to optimal order accuracy and efficiency for
most practical flow models The work presented here repre-
sents an attempt to overcome some of these difficulties.
Many impediments are caused by local phenomena, such
as shocks and boundary layers, where higher resolution and
accuracy are seriously needed. To meet this need, FAC is
used as an efficient multilevel process for local
refinement. FAC implicitly solves nonuniform grid problems
by solving a series of problems defined on uniform grids.
Therefore, FAC takes advantage of uniform grid properties
and it is easily coupled with existing codes and conven-
[


-2
tional multigrid methods. Because FAC originates from the
multigrid discipline, it can achieve enhanced accuracy at
apparently optimal order efficiency. This thesis focuses
on the use of FAC in conjunction with MG for solving
Poisson equations, EGG equations, potential flows and
cavity flows.
Accuracy and efficiency of FAC and MG solvers depend
critically on the^quality' of discretization, which can be
problematic for many highly active flows and locally
refined grids. To provide effective discretization schemes
for high Reynolds (Re) number flow equations, we developed
the rotated hybrid and finite volume element (FVE) schemes.
A significant part of this development is the answer to
the question of how to treat such aspects as Dirichlet and
Neumann boundary conditions, systems, anisotropic problems,
singularities, and discretizations and transfer at boun-
daries, interfaces and patches.
The difficulties involved with solving incompressible
Navier-Stokes equations are clarified in [27]. Our source
term correction method (STC) represents an attempt to over-
come some of these difficulties in a very convenient .way.
Several of the methods and results presented here have
already appeared in published form. See [23-26]. This
thesis summarizes this work, develops several additional
techniques (e.g., the time-marching schemes and STC), and
presents more recent numerical results.
2. Arrangement of the thesis
In Chapter II, we take the Poisson equation as an
example to illustrate the basic principles and schemes for
the methods of MG, FAC and FVE. The computational results
show that local refinement can greatly enhance discretiza-
tion accuracy and that these methods can compute an
acceptable numerical approximation very efficiently.


-3-
In Chapter III, we apply full multigrid (FMG) with cubic
interpolation to solve the elliptic grid generation (EGG)
equation and use FAC "to generate a smooth locally generated
patch grid.
In Chapter IV, we apply MG and FAC to solve a full
potential flow around a NACA0012 airfoil. For this applica-
tion, discretizations in the global grid are different than
those in the patches. Computational results presented there
show that MG can substantially accelerate convergence and
the FAC local refinement process can greatly improve
accuracy, especially for problems with shocks.
In Chapter V and VI, we apply MG and FAC to driven
planar cavity flow (incompressible Navier-Stokes
equations). Reviewing the weakness of classical finite
l
difference theory and conventional notions on artificial
viscosity, we develop a so-called rotated hybrid scheme
(Chapter V) and an FVE scheme (Chapter VI) for medium and
high Reynolds number flows.
In Chapter VII, we establish an equivalence between
spatial relaxation and time-marching schemes, then use
this interpretation as a basis for development of MG for
pseudo unsteady flow .
In Chapter VIII, we develop a so-called source term
correction method (STC) for incompressible Navier-Stokes
equations. ^
(


CHAPTER II
POISSONS EQUATION
The purpose of this chapter is to introduce the
basic concepts and illustrate some of features of
multilevel adaptive methods in a simple setting. For this
we restrict discussion to the model Poisson equation.
1. Equation and boundary conditions
f=-6(0,0)+6(l,l),
where 6(x,y) is the Delta function at point (x,y). This
is the so-called five-spot problem, which represents
the pressure equation in an idealized oil reservoir. The
term -SC0,0) corresponds to an extraction well and 6(1,1)
to an injection well.
V*Vp=f
n. 7 Jn n=o
on dfi (no flow boundary condition)
(compatibility condition). (1)
in fl
2
Here, fi is the unit square in 3? and
Control volumes
Finite elements
Injection well
Extraction well
Figure 1. Triangular elements and control volumes
for a uniform grid


-5-
2. Finite volume element (FVE) discretization
Let n. be a control volume in as depicted in Figure
1. Then, by integrating (1) over fl- and using the
divergence theorem,
JV, V*Vpda = JV, n.Vpds,
W i 1 i
we have
n.Vpds = Jq fda (2)
1 i i
Here, n is the outward unit normal defined at points of I\ .
(We use da and ds, respectively, to denote "volume" and
"surface differentials.)
The FVE discretization of (1) is now defined by
restricting the unknown ip in (2) to the subspace of
piecewise linear functions defined on triangular elements,
as depicted in Figure 1. We illustrate how this produces a
matrix equation for the nodal values of the finite element
function at the following sample points.
A. Interior points (Figure 2): The "flux" fl along the
northern half of the east segment, e of T. is given by
2 1
n.Vpds = f
Pxdy
(3)
Assuming that p is piecewise linear, then in element I we
may write
ip = ax+by+c .
Note that a =ip in I so that
x
fle =ia-2 =i('pE"^P) '
2 i
We similarly obtain flg flw flw fln fln flg ,
(4)


-6-
(
and fl Since jV, n.Vpds is the sum of the fluxes at
8 2 r _
interior points, we obtain the familiar five-point
difference formula
Figure 2. Sample interior point for IV
(Relevant points are labeled with capital letters,
elements with Roman numerals and-volume boundary segments
with subscripted lower case.)
B. Sample west side boundary points (Figure 3): Since
nV*s-V^. ".rt'W f'rWV- and flVfV
(this results from the homogeneous Neumann boundary
condition), then we have the scheme

(6)


to 11-1
-7-
N
Figure 3. West wall point for I\
Figure 4. Southwest corner point for I\
1.
C. Sample southwest corner point (Figure 4): Since flg=
f!n= 2(PN_PP) flw=fls= and Jft fdz="Jn ^(.)dz
i i
=-l, then the scheme is
PP" 1
(7)
Other grid points of 0 are treated similarly. We
write the resulting equations as


-8-
Lh'ph=fh,
where represents the vector of nodal values
approximating the solution

(8)
3. Multigrid algorithm (FAS)
To solve (8), we use the full approximation scheme
(FAS, [3]). One two-grid FAS cycle is described as
follows. (Here, uis the number of relaxation sweeps
performed before the coarse grid correction, and u the
A
number after correction.)
Step 1: Relax on using v^ Gauss-Seidel sweeps.
Step 2: Solve VTf (f,-^).
Step 3: Make the correction
Th , ,2h -
i=h V'sh'^h-h V
1
Step 4: Repeat step 1 using v sweeps.
Here, I^1 is the injection operator and 1^ represents
bilinear interpolation. The restriction operator 1^ is
derived from viewing FVE as a conservation scheme. Here,
for the uniform grid case it is the full weighting
operator [3], (See [21] for more details.) Specifically,
ijj*1 uses the stencil
rill-,
4 2 4
1 1
2 2
111
L 4 2 4 J
at interior points (Figure 5). This is derived by
considering the fraction that a fine grid volume inter-
sects the coarse one. For example, because the fine grid
(10)


point NE has one-fourth of its control volume located in-
side the coarse grid point volume at P, NEs contribution
to P in is In a similar way, the stencil at a west
side boundary point (Figure 6) is
P l_
2 4
i i
1 i
-2 4
(11)
72h
At the southwest corner, 1^ uses the stencil (Figure 7)
<- i. L
2 4
* 7
(12)
Figure 5. Restriction at an interior point


Figure 6. Restriction at west wall
I 1 /
. NE
1/ 4
p E

Figure 7. Restriction at the southwest^ point
Our experiments use V (2,1) cycles by which we mean 2,
1/2=1. and (9) itself is "solved" by a recursive use of a
/
multigrid cycle. For meshes of 33 x 33, 65 x 65, and 129
x 129, we achieve an asymptotic per-cycle convergence
factor of 0.090. (This can be improved significantly by
using a better relaxation scheme. For example, zebra-line


-11-
Gauss-Seidel achieves an asymptotic convergence factor of
0.016, while double-line Gauss-Seidel, which solves two
lines at a time, achieves a factor of 0.0083.)
4. FVE accuracy
Using the solution of the 129x129 mesh as an approxima-
tion of the exact solution of (1), we compared it with the
solutions of (8) obtained on coarser grids to get
approximations of discretization errors. As a measure, we
used the discrete energy norm defined by
1
HIThlll=2 <13)
where <,> denotes the Euclidean innerproduct. Here we
use v^=Pj-p., where is the solution of (8) on grid h and
tp is the "exact" solution (i.e., solution of (8) with
h=j-|g injected to grid h). We compute this "solution" by
performing V-cycles until the changes in the solution are
insignificant.
Table 1 shows that these errors decrease by a factor
of about 4 as h decreases by a factor of 2, suggesting that
FVE obtains 0(h ) discretization error (in the discrete
energy norm).
Table 1
FVE Discretization Errors
Mesh 17X17 33x33 65x65
Trunc. errors 0.113e0 0.391 e-1 0.994e-2
5. Nested iteration (FMG, cf. [3])
Starting from the coarsest grid (3x3) and proceeding
to the finest, we perform one V-cycle on each level in


turn, with the result cubically interpolated to the next
finer grid to start the V-cycle there.
Table 2
Accuracy of Nested Iteration
Mesh 17x17 33x33 65x65
Discret. errors 0.113e0 0.391 e-1 0.994e-2
Actual errors 0.133e0 0.400e-1 0.108e-1
This nested iteration scheme'is a semi-direct method in
the sense that its aim is to achieve an "actual error"
(discrete energy norm of the difference between the
computed solution and the solution of (1) injected to the
grid) comparable to the discretization error. Table 2
- supports the success of this objective: the actual error
in each case is within a small fraction of the
discretization error. The cost of this scheme is about
five "work units" (i.e., the cost-equivalent of about five
sweeps of relaxation).
6. FAC for nonuniform (composite) grids (Figure 8)
The singular points of

suggest the use of local patches for increased resolution
and accuracy of the solution. Figure 8 depicts the case of
two refinement levels in each corner, although for
discussion we restrict ourselves to the case of a single
patch in the NE corner (Figure 9).
A. Interface point (Figure 9 and Figure 10)
To ensure continuity, the slave points are treated by
enforcing the element functions to have values equal to the
average of values at the neighboring interface points. Now
the flux on the east segment of the surface depicted in
Figure 10 is given by


fl =fl +fl +fl +fl
6 C1 S S %
=pe^ne^se^n4Â¥>s'Ipp-
Other fluxes are flw=pw-pp, f'ln=f(PN"Pp) and fls=4^S_??P^
We' then have the equation

(14)
Figure 8. Sample three-level composite grid
H Boundary point
O Global grid point
(j) Patch point
X Slave point
PI Interface point
Finite element
Finite Volume
Figure 9. Triangulation and control volumes for a
composite grid


Figure 10. Vest side interface point
We can similarly derive schemes for other interface
points. On the south side of the interface we have
at the southwest corner we have
^V^^NE^V^S^SW1'4*?'0' ' (16)
at the northwest corner we have
(17)
and at the southeast corner we have ;
<18)
B. Restriction at interface points r
s
For the west side interface point (Figure 10), the
~2h
stencil for 1^ is given by
' 0
1
4
J_
2
_1
4
(19)
0


-15
where 1^
represents the restriction used to transfer
composite grid residuals to the global (coarse) grid.
The stencil for ij^1 is
1 i i
4 2 4
0 10.
at the south side interface,
1
0
r
(20)
(21)
at the southwest corner,
(22)
s
at the northwest corner, and
at the southeast corner.
C. FAC
For a detailed description of FAC and its motivation,
see [21]. FAC solves the nonuniform (composite) grid problem
by iterating on a series of uniform grid problems. Thus,
FAC is convenient for adapting to the existing codes, for
providing local refinement and for coupling with multigrid
as the uniform grid solver. This contributes to FACs
accuracy and efficiency. We just describe one cycle of a
two-level FAC schemes as follows:
J
Step 1: Do one FAS V(2,l) cycle on the global grid to obtain
an initial guess for the global grid and local
patches.
Step 2: Compute the composite grid residual and transfer
it to the global grid equation to FAS-correct the
right-hand side. Solve the resulting global grid


-16-
equation by one FAS V(2,l)-cycle.
Step 3: Using the global grid approximation to define boun-
dary values, solve the local patch equation by one
FAS V(2,l)-cycle. Go to step 2 if the solution has
not yet converged. J
Table 3 displays the convergence history of FAC for
the three-level example illustrated in Figure 8. Here we
use a 33x33 global grid and two nested 33x33 patches in
each corner. Displayed are the residuals computed for the
composite grid as well as the residuals computed for patches
after each patch solution. The geometrically averaged
convergence factor is computed by taking the fourth root ^
of the ratio of the residuals after the fifth and first
cycles.
Table 3
Convergence History of FAC Using Three
Levels (Two Patches) of 33x33 Grids
Cycle Composite Level 2 NE patch Level 2 SW patch Level 3 NE patch Level 3 SW patch
1 0.113e-3 0.638e-4 0.300e-4 0.468e-3 0.220e-3
2 0.158e-4 0.304e-4 0.796e-4 0.116e-3 0.106e-3
3 0.219e-5 0.342e-5 0.122e-4 0.140e-4 0.189e-4
4 0.329e-6 0.531 e-6 0.160e-5 0.205e-5 0.242e-5
5 0.513e-7 0.952e-7 0.189e-6 0.265e-6 0.181 e-6
Factor 0.146
D. FVE and FAC composite grid accuracy
To test the accuracy of the FVE discretization on a
composite grid, we compare the discrete approximation
(computed after many FAC cycles) to the 129x129 grid
solution. Again we use the discrete energy norm. Table 4
depicts the results for two- and three-level composite
grids. Note that the error for the case of three levels of


-17
17x17 grids is about half that for the 33x33 single-level
case, although the number of grid points is almost the
same. Note by observing the two-level errors that they
appear to decrease by a factor of about four when each
grid is refined by a factor of two.
Table 4
FVE Discretization Error
Mesh 9x9 17x17 33x33 65x65
Uniform grid 0.232e0 0.113e0 0.391 e-1 0.994e-2
2-level comp, grid 0.174e0 0.581 e-1 0.148e-1
3-level comp, grid 0.712e-1 0.184e-1
We also tested a nested iteration version (or full-
multigrid-like version) of FAC. Here, we used one cycle of
FAC on the composite grid consisting three 5X5 meshes
(one global, two local) and cubic interpolation to obtain
an initial guess for one cycle of FAC on the composite
grid of three 9x9 meshes. This in turn yielded an initial
guess for one FAC cycle on the grid of three 17x17 meshes,
and so on until the finest composite grid was reached.
Table 5 depicts the results. Note the analogy to Table 2
for the uniform grid case. As in that case for multigrid,
this form of FAC yields a semi-direct solver costing about
five (composite grid) work units.
Table 5
Accuracy of FAC Nested Iteration
Mesh 9x9 17x17 65x65
Discretiz. errors 0.174e0 0.580e-1 0.i48e-1
Actual errors 0.182e0 0.598e-1 0.153e-1


-18-
E. AFAC
An asynchronous version of FAC, called AFAC [28], is
designed to allow simultaneous processing of each level of
the composite grid. This is accomplished by subtracting an
error component from the patch grid 'problem that is >.
duplicated on the global grid equation, which can be done
by a simple modification of the V(2,l)-cycle.,(See [21] and
[28] for more detail.) The geometrically average convergence
factor for AFAC achieved for the 33X33 three-level composite
grid example was 0.283. Note that this factor is worse than
the corresponding factor (0.146) depicted in Table 3. This" \
is to be expected since AFAC can be interpreted as a Jacobi-
like version of the Gauss-Seidel^-like FAC process.
\


-19-
CHAPTER III
ELLIPTIC GRID GENERATION (EGG)
The purpose of this chapter is to study the behavior
of MG and FAC as fast solvers for the EGG equations. Of
course the objective of EGG is to produce grids that better
conform to.the nature of physical system by providing a more
effective grid structure on which to solve the flow
equations. Thus, MG and FAC have the dual roles of solvers
for both the EGG equations and the flow equations on the
grids that EGG produces. This chapter is concerned only with
the EGG equations themselves; the flow equations on the
EGG-generated grids will be treated in subsequent chapters.
1. EGG equations
Elliptic Grid Generation (EGG) was developed by
J. Thompson [1] in the seventies. The purpose of EGG is to
find a sufficient smooth coordinate mapping £=£(x,y), rj=rj(x,y)
which transforms the physical domain into the computational
domain. The Poisson equation with proper boundary conditions
represents this kind of mapping:
v2e=p
VVQ (24)
where P and Q are called control functions which can help
to realize grid concentration. Actually it is the inverse
mapping x=x(£,r}), y=y(^,r/) that is needed in practice. It
('
can be used to map a rectangular grid in the computational
domain to a smooth but irregular boundary-fitted grid in the
physical domain.
As a prototype problem we chose uniform inflows with
zero attack angle around the symmetric airfoil NACA0012
(Figure 11) and used EGG to generate a body-fitted grid
J
(Figure 12-13). Assuming P=Q=0 for simplicity, we may write
the EGG equation for defining this inverse mapping as


-20-
wv
ay^-2^^0
where
(25)
HX- X+yf' d ^Wyfy
These equations are nonlinear and elliptic; the
boundary conditions are Dirichlet in this case since they
specify the location of the boundary grid lines.
Although EGG is an elliptic problem, only few others
have used the MG method for EGG and noone else has used FAC
to generate smooth patch grids.
Figure 12. EGG physical grid


F
I III 1U1.11 111 II11 llLl.il 111 \ \----
A B C D E
Figure 13. EGG computational grid
2. Finite difference equations
The finite difference equation corresponding to (25)
can be obtained by using central differences on the
derivatives as follows:
XE+XW~2xP _2p XNE~XNW~XSE+XSW + _xN*XS~2xP =0 (2g)
h2 4h2 h2
where h is the mesh size in the computational domain.
Equation (26) may be written in the general form
apxp=aexe+awxw+anxn+asxs+anexne+anwxnw+asexse+aswxsw
(27)
a
,=A__=- S-.
A __
* AxmT"AOE'
where AN~AS-h2 ^^SW- 2h2 "NW"SE_2h2
and Ap=2^Q!+^ Note that (27) is a nonlinear algebraic
' r h2
equation since a, /? and 7 are functions of x and y. We can
develop difference equations for y in a similar manner.
3. Initial guess
Equation (27) can be solved by an iterative process,
but a good initial guess is usually necessary for


-22-
efficiency. We study two possible strategies for this. The
first is to use transfinite interpolation [1] from the
boundaries into the field. With indices l where i=l or I or j=l or J refers to the boundary, then we
determine the initial guess x from the boundary values x
by
x(i,D= ilj[ x(I,j)-t~|"K(l,j)+ ^ x(i,J)+j5jx(i,l)
i-1 j-1 /T i-1 J-j /r i\ '
1-1 J-l 1-1 J-l x^1-1^
The second approach we consider is full multigrid. The
central idea here is to start on very coarse levels using
a basic multigrid cycling scheme to obtain a good
approximation there, then cubically interpolate this
approximation so it~ can act as an initial guess for basic
multigrid cycles on successively finer levels.
4. Full approximation scheme (FAS)
Gauss-Seidel relaxation for solving system (27)
typically stalls after a few iterations. This is because
Gauss-Seidel, though effective for high frequency errors
components, has very little effect on low frequency
components. Multigrid capitalizes on this "smoothing"
property of Gauss-Seidel by visiting coarse grids to
resolve smooth errors. To accommodate the nonlinearities
in (27), we can use an FAS version of multigrid with
bilinear interpolation and full weighting of residuals and
approximations.
Table 6 shows the results of various V-cycles applied
to (27). In this example we used four grids. The finest was
33 X 17 and the coarser ones were 17 X 9, 9X5 and 5X3.
We depict results of six FAS cycles in terms of residual
norms (we display the residual norms, Rx, for x only
because the higher x resolution means that the y residuals


-23-
are naturally much smaller) and required "work units".
This is shown for four different V(i/ ,v ) cycles, where u
12 1
and u are the numbers of relaxation sweeps performed
before and after coarse grid correction, respectively. /
Note that V(l,l) seems most efficient, with an average per
per work unit residual reduction factor of about 0.51.
Table 6
Comparison of Various V-cycles Applied to
The EGG Equation
FAS V (1.1) V(1.2) V(2,1) V(2.2)
cycle R X i WUs R X WUs R X WUs R X WUs
1 .952e-7 2.813 .179e-7 4.063 .291e-7 4.063 .949e-8 5.313
2 .641 e-8 5.625 .134e-8 8.125 .262e-8 8.125 .678e-9 10.63
3 .819e-9 8.438 .107e-9 12.19 .230e-9 12.19 :54e-10 15.94
4 .101e-9 11.25 .19e-10 16.25 .50e-10 16.25 .17e-10 21.25
5 .26e-10 14.06 .67e-11 20.31 .19e-10 20.31 .68e-11 26.56
6 .82e-11 16.88 :20e-11 24.38 .73e-11 24.38 .27e-11 31.88
5. FAC for local refinement >
In the vicinity of the leading edge of the airfoil
and at shock waves produced by the flow problem, velocity
gradients can be very large. Local refinement can be
effective for providing more accurate results and finer
details of the solution. In the EGG context, it can be
very important for computing smooth locally generated grids.
Here we use FAC, which we test by placing a local 49 X 13
grid about the airfoil as shown in Figure 12.
Table 7 shows residual norms for the discrete EGG
equations on each grid (both the coarse and fine; here we
use a single multigrid V(l,l)-cycle as the basic grid
solver); we have depicted residual errors for each. The
y
correct error measure, however, is the composite grid


-24
residual norm which we have also shown. The work units are
measured in term of the composite grid equations, hence the
correspondences with V(l,l) in Table 6.
6. Discretization error estimates
In order to examine the accuracy of the numerical
solution, we used the following test problem:
4 2 2
0!X^-2/?x^+7x^=27r simr£cos7TT7(sin mj+cos tt£)
4 2 2
ocy^-2fiy^r+'jy^r=2v cosx^sin7rr?(cos X77+sin 7r£) (29)
£=0: x=0, y=sin?r£;
£=1: x=0, y=-sinjr£; ,
r/=0: x=sin7TT7, y=0;
77=1: x=-s in7TT7, y=0.
Its analytical solution is x=sin7r£cos7T77-and
y=cosjr^*sin5TT7. We then compared the results of applying
several V(l,l)-cycles with this solution on the three
meshes depicted in Table 8. We used L norms and average
00
deviation estimates in both x and y directions. Note the
2
very apparent 0(h ) behavior of these errors.
Table 7
Convergence History of FAC for^EGG
FAC cycle R X coarse R fine R x composite WUs
1 .683e-8 .164e-10 .720e-8 2.813 r
2 .952e-9 .152e-10 .102e-8 5.625
3 .311 e-9^ .157e-10 .351 e-9 8.438
4 .141 e-9 .154e-10 .161 e-9 11.25
5 .560e-10 .150e-10 .685e-10 14.06
6 .212e-10 .146e-10 .357e-10 16.88
j


-25-
Table 8
EGG Discretization Error Estimates
1 J hx hy max(errx) max(ery) ave(erx) ave(ery)
17 17 0.0625 0.0625 0.0065 0.0060 0.0032 0.0032
33 33 0.0313 0.0313 0.0016 0.0015 0.0008 0.0008
65 65 0.0156 - 0.0156 0.0004 0.0004 0.0002 0.0002
r\
j
i
)
V


-26-
CHAPTER IV
TRANSONIC POTENTIAL FLOWS
A
Although there is already substantial literature on
potential flows, our work is unique in its use of FAC with
different discretizations on different levels. The accuracy
and efficiency of this method seems to signify substantial
improvement.
(
1. Governing equations
The full potential equation in strong conservation
form may be written as
d(pa). dipv) n
+ ~W~ "

(30)
where p is the density of gas, 7 is the specific heat ratio
and is the Mach number at infinity (i.e., the
undisturbed fluid velocity normalized by the speed of
sound). For potential flows, we may write the velocity
components as u=y? and v=p where

x y
potential.
In computational coordinates (£,77), which are defined
in Chapter III, (30) may be rewritten as
^-(pu)+^(pV)=0
P=[l+ ^ Mja-^+Vp^/Jj}]17^0,
where
U=A ii^£ 12 22*77 12*$
2 2
x +y
A________*? .*? a -
A--~ T Aoo
2 2
11
22
(31)
r


-27
\2=

a-d Jl=
2. Boundary conditions (Figure 14)
On the wing surface and symmetric axis, the normal
velocities are zero:
V .
n an
(32)
At downstream boundaries, free flow conditions are used:
(33)
du d2

or r=0

Upstream and top boundaries are treated as far field, where
we assume zero disturbance:
u=u and v=0 . (34)
00
3. Discretization of the potential equation
The conventional finite volume method is a convenient
and accurate way to construct difference equations for
points of the interior and boundary, especially for the in-
terface points at the refinement regions. We discuss several
cases in the following, where for simplicity we assume that
I
i


-28
each grid (whether global or local) is assumed to be uniform
with equal mesh spacing in the £ and 77 direction.
A. Interior points (Figure 15)
We rewrite (31) as
f Vpwda = [ pw*nds=0 ,
Jn Jr
(35)
where w is the vector function with components U and V, (1
is a small control volume centered around an interior point
j j). r is the boundary of (1 which is assumed to be
polygonal, and n is the outward unit vector normal on F.
For any given line segment I^of T, we denote the flux
flr =1 pwnds .
0 JrA
0
(36)
By our assumption on grid uniformity (for the computation
plane (£,?)), we may take II to be a square centered at P
with sides of length h=£. .-£. =n. n. .(assume
e M+l,J 'i,J+l 'i,J
for simplicity that h=l for the finest grid) that is
aligned with the grid. Then using subscripts to denote the
east,west, north and south segments of T, we can rewrite
(35) as
fl +fl +fl +fl =0
e w n s
(37)
Let (/>U)e denote some approximation to pU along the
east boundary segment of H (e.g., (pU)g is some average
value of pU along e) and similarly for (pU) (pV) and
(pV) Then (37) may be approximated by thb equation
s
(pU)e-(pU)w+(pV)n-(pV)s=0 .
(38)


-29-
N
n
w IB e

s
S
Figure 15. Finite volume method for interior
points
Since
B=AiirAi2'1>
V-A22*l) A121<,f
then we may use central difference approximations to the
fluxes, yielding
(PAH) i+i, j1l, j >-<'* > nj, j , j+l+
-(M,,) , y>+<**> i_i, j ->1, J+i+
pi+l, .j+r^i-l. j1
~ /


-30

= o ,
(39)
where
(pAi l} i+J-, j"2[ (pAil} i+1, j+(pAi lJ i, j]
i

and so on. This can be written in the general form
ap^p=ae^e+aw^w+an^n+as^s+ane^ne+anw^nw+ase^se+asw^sw
where (40)
V^iA+i-J 4t(pAi2)i,j+i- (pAia)i,j-j1
V(pAu)i^-,j i[fpAi2)i,j-4 (pAi2)i,j^
2 l *
4tCpA12)i+|,j
ANE="iC j+i j+(PAi2 > i Jti1
2 2
ANw4[ ASE4[ *PA12 i+A, j+ 2 2
\ASVT"4[ (pAi2} i-f f j+(pAi2) i, j^-3
2 Z
and Ap= Ae+AW+AN+AS+ANE+ANW+ASE+ASW -
B. Solid wall points (Figure 16)
Because the south boundary is a solid wall, then we
have fl =0. From (37), we then get
^ s
fl +fl +fl =0
e w n
(41)


-31

which is approximated by
(pU)e-(pU)w+(pV)n=0
or
(42)
^(pA,,). 1 .(p. .-p. .)-(pA ). 1 .t(^- 1 ,+
2 ^ ll i+-,j *1+1,3 *i,j y 12 i+-,3 4 *i+l,3+1
p. ,-p. , ,-p. .)
*1,3+1*1+1.J i.J
~r(M. .) l .(p. .-p. .)+(pA). l .t(^-
2 ^ ll l-r,3 *1,3 *i-l,3 l2_;i-r-,3 4 *1,3+1
p. , . ,-p. .-p. , .)
M-l,3+lnKi,3nKi-l,3/
+(^5 i. J*(pi. . j:-''n1 i, j+^i+i. j+i+
p. , .-p. , . ,-p. , .) =0 .
*i+l, 3*1-1,3+1^1-!, 3
(43)
- f i N
n
w P e
'/////////7V77.
Figure 16. Finite volume method for solid wall points
-This we write as
appp=a#e+a#w+an^n+an#ne+anv/pnw (44)
where
1} i4, j*i[ (M12 i . 2 >'
2 2
I jl.
2J
I
/


-32-
and
ap_ae+aw+an+anepne+anw
C. Other boundary points
For upstream and top boundary points, we simply keep
the

For downstream boundary points, we implicitly embed
the free flow condition

xx v
downstream boundary).
D. Composite grid interface points (Figure 17)
For the composite grid we must be careful in our
treatment of the interface points between the coarse and
fine grid points. We use (37) to develop the following
equation for point P (Figure 18):
P=V+V
(45)
1(PN-Pp) (PA )]


-33
+ (pA12)i,j*[Â¥5SE 4('PS+PP+PSW+,PW)]
J 2 2
=0. v (46)
Here, superscript f is used to specify fine grid
quantities. Now this we write as:
appp=a#>w+an¥?n+as'ps+an#nw+aswV5sw
where
and
+AEfE+ANlftlE+ASE?SE ,
V('AX 1 > i-|. ( i. i. '
i. xi-i, x > i, j xi]
V^xx > i, j-l-71^1214, J+(^.x:1 i. i-i1'
(x>A, j) j J., j+ j j+lJ
Asr~il (,,AI x* i-i-, j+l(
^^Ix'iAj^lxl.j-i'
VWVWVV'se
(47)
/


o
X
Boundary point
Global grid point
Patch point
Interface point
Slave point
Figure 17. Composite grid
Figure 18. Difference schemes for interface points
E. Computation of the density p
In order to calculate the coefficients of

(44) and (47), we should compute the density p in advance.
For interior points, we use the expression
2rM[1-(AnprAi*V!{/Ji
l
-(A22VAi*VVJi]}>1- (48)
At the far field boundaries, p is set to unity. On the


-35
solid wall, the zero-normal velocity condition should be
used, namely,
wn=w.V[y-f(x)]=0, (49)
where y=f(x) is the airfoil contour function. Hence,
or
-Pxf(x)+py=0
[ftxJy^+x^^-tfUJy^+x^p^O,
which implies

ruK*\ tn
rj f(x)y*+x.
Vt
(50)
rrAe "22 >
Thus, at solid wall points we use the expression
(51)
4. Stability and artificial density
The discretization described above is stable and
second-order accurate for subsonic problems. For transonic
problems, however, to maintain stability we use an
artificial density concept which is common in practice
[14]. We replace p by p where
2
P=P~uoUl-e)^p+e^pl r (52)
and
r M 2,
l/0=2.0*Max|0, l-(jp) I . (53)
Here, M is the so-called cutoff Mach number (where M c c
1/=0), M is the local Mach number (i.e., the fluid
2 2l/2
speed normalized by sound speed: M=(u + v ) /c), and e
is a free parameter in the interval [0,1] ( in most of our


-36-
example, at points where u>0,s we use the approximation
5. Multigrid for potential flows
The multigrid scheme used to solve the potential flow
equations is similar to that used in solving the EGG
equations, namely, an FAS V(2,2)-scheme. We list a few
important differences:
A. Using (40) to construct the coarse grid operators
L2*1, the quantities A2*1, A2^ and A2^ which depended only
on physical node distribution, can be defined simply by
restricting the corresponding grid h quantities. On the
2ll 2Tl
other hand, p, which depends on the unknown

be calculated by (48) and (52). Initially, (48) would be
evaluated using the starting guess p2^= I^V*.
B. In the transonic flow region we use an artificial
viscosity which is related to the local Mach number M. For
compatibility in the solution process, we geometrically
2ll
align the transonic regions on all levels by defining M
to be the restriction of M*1 to grid 2h. (Careful
treatment of the Mach number as a nonlinearity could lead
to a more efficient scheme, especially for the case of
v, sharp shocks. However, this approach is potentially much
more complex with several open questions remaining.)
C. A critical aspect of the solution process is
3
(54)
where the £. are on the finest grid.
*- J
V.


-37-
careful treatment of the boundaries. For Dirichlet
2h
boundary conditions, we only need to use 1^ for all
fine-to-coarse transfers of

need to use I?*1, L2^or l\. at the boundaries. At Neumann
boundaries, we should use the difference operators L?1 and
L2\ the coarse-to-fine interplation operator 1^, and a new
for residual transfer. In particular, for ij^1 we use the
6-point stencil (Figure 19)
ll l
16 8 16
JL i
8 2 8
(55)
for the restriction of the residual to the coarse grid at
the wing surface.
C
Figure 19. 6-point stencil points for residual
transfer on solid wall
D. We use different discretizations by choosing
different cut-off Mach numbers M on different levels.
c
Actually, we choose Mc=0.85 on the global grid and Mc=0.95
\


r
-38-
on the patch. In this way we obtain a sharper shock without
a significant increase of work units.
6. Computational results
Three flow examples around the NACA0012 are investi-
gated using a global grid with relaxation only as well as
s'
with a 4-level FAS V(2f2)-cycles as the solvers. The free
stream Mach numbers used for these examples were M =0.0
00
(incompressible), 0.5 (subsonic), and 0.75 (transonic).
The finest grid was 33 X 17. Comparisons of the
performance of relaxation and multigrid are depicted in
the respective Figures 20,21, and 22. Convergence is good
for incompressible and subsonic flows where the residual
reduction factor per multigrid cycle (equivalent to about
5.625 work units) is less than 0.1. In particular, it
costs only 28.125 work units to achieve residual norms of
8 8
0.3X10 and 0.42X10 for Mw=0.0 and 0.5, respectively.
(These experiments were made without the use of FMG to
compare with earlier work.) However, our implementation
of multigrid is slower for transonic flow. Here, 45 work
units are needed to achieve a residual norm of 0.62X10
r M ,
with an artificial viscosity of i/=2.0*Max 0, 1Cj^-)2 I and
a cutoff Mach number of Mc=0.85. Increasing the cutoff
Mach number leads to a sharper shock wave at the expense,
however,, of more work units.
To test the effectiveness of FAC, we placed a 49 X 13
grid about the airfoil and ran several cycles, using the
multigrid solver described above on each grid. Note that
FAC performance on this composite grid is similar to
multigrid performance on the global grid. To test
resolution, we graphed the pressure distribution on the
wing surface in Figures 23 and 24, respectively, for the
c


-39-
global grid and for the composite grid tests. The greater
accuracy obtained by FAC is evident in the increased
sharpness of the shock wave and attenuation of the
oscillation behind the shock. ^
The complex body configuration may cause instablity in
the numerical calculation when we use conventional adaptive
grid generation methods. FAC can solve this problem by
generating smooth patch grids. Using different
discretizations in global grid and patches, we obtain
improved accuracy without significantly increasing cost.
Log(Residual)
Figure 20. Multigrid convergence history
(grid: 33x17, M =0.0)
oo


-40
Log(Residual)
Figure 21. Multigrid convergence history
(grid: 33x17, M =0.5)
00
Log(Residual)
Figure 22. Multigrid convergence history
(grid: 33X17, M =0.75)
> oo


-41-
Figure 23. Pressure distribution produced by
multigrid (grid: 33X17, M^O.75)
Figure 24. Pressure distribution produced by FAC
(global grid: 33X17, M =0.75:
oo
refind grid: 41X13)


-42-
)
CHAPTER V
THE ROTATED HYBRID SCHEME FOR
PLANAR CAVITY FLOWS
Although much progress has been made in computational
fluid dynamics, there is still a critical need to develop
better computing methods for viscous flow problems, espe-
cially in the case of high Re numbers. The vorticity-stream
function formulation is still considered as a favorite
approach for investigating 2-D incompressible viscous flows.
Unfortunately, for high Re numbers, central difference
approximation to the vorticity transport equation is
unstable. The upwind scheme is stable, but it tends to
introduce too much artificial viscosity and accuracy is
degraded [30]. Spalding [31] suggested a hybrid scheme
which tries to take the advantages of central and upwind
schemes. We developed a modified hybrid scheme [32] which
is based on a certain one-dimensional analytic solution.
However, even these schemes have trouble at high Re number
flows. Chen [33] developed the so-called finite analytic
method, with very good results even for high Re, but the
method is quite expensive. Here we introduce what we call
the rotated hybrid scheme, which achieves the accuracy of
the finite analytic method at substantially less cost.
For fast solution of flow problems, we use the
multigrid method. Here, MG uses a collective
relaxation scheme that treats the vorticity and stream
function equations in coupled^form, which is especially
important near the boundary. This is one of the key points
j 1
for achieving"full multigrid efficiency. For anisotropic
problems, which are typical of high Re number flows, we use
alternating direction block-line Gauss-Seidel. To obtain
details of the secondary vortices near corners of a square
cavity, FAC is used for efficient local grid refinement.


-43
1. Differential equations and boundary conditions
Incompressible 2-D Navier-Stokes equations in
vorticity-stream function form may be written as
Re.u.£x+Re.v.£y-£xx-£yy=0 . (56)
t +£=0 , (57)
rxx ryy *
where Re is the Reynolds number, u and v are the respective
x- and y- components of velocity, £ is vorticity and ^ is
2
the stream function. Here we consider the domain fl=[0,l] .
For a driven cavity flow assuming the zero stream line is
along the wall with the non-slip condition imposed on the
wall, we write the boundary condition as follows (Figure 25):
^M), u=0 and v=0 for x=0
^=0, u=0 and v=0 for x=l
iM), u=0 and v=0 for y=l
4*= 0, u=l and v=0 for y=0. (58) 2
Figure 25. Cavity flow
2. Conventional difference schemes for the vorticity
equation
Because equation (57) in terms of the stream function
is; just a Poisson equation, which is easily treated by
finite differences, we will focus our attention only on


-44-
(56), the vorticity equation. If we locally linearize the
coefficients and let A=Re.u/2 and B=Re.v/2, then (56) may
be written as
2A.£ +2B.£ =£ +£ (59)
sx *y sxx syy
on which we may use central differences (Figure 26):
^N^S ^E+^W-2^P ^N+^S_2^P
^E~^w
2A-|r*+ 2B-
2h
(60)
Wo
N
, (
/ { P ^
6 E
Figure 26. Difference scheme
We can write (60) in the general form
ap^p=ae^e+aw^w+an^n+as^s,
where Ag=l-Ah, A^=l+Ah, V1 -Bh, Ag=l+Bh, and Ap= A^+A^+A^
+AS=2. The central difference scheme is unstable for large
Re. In this case, we may consider upwind differencing
which, for u>0 and v>0, is given by j
^P^w rtT, ^P~^S
2A;---+ 2B*-
£e+£w 2^p ^n+^s~2^p
. 2 + ,2
(62)
(63)
h T h "
This has the general form
ap^p=ae^e+aw^w+an^n+as^s,
where Ag=max(l-2Ah,l), A^=max(l+2Ah,1), A^=max(l-2Bh,1),
Ag=max(l+2Bh,l), and Ap=AE+A^+A^+Ag. Although this upwind
scheme is1 quite stable even for large Re,Lit produces a


solution that is far away from the physical one. This can
be seen in Figure 27 for our driven cavity flow. In the
next sections, we will consider difference schemes that
offer more accuracy. We start in the next section by
developing the 1-D case.
Figure 27a. Physical
solution (Re=2000)
Figure 27b. Upwind differece
solution (Re=2000)
3. Improved difference schemes for 1-D vorticity
equation
The locally linearized 1-D vorticity equation can be
written as
ok >
*x 2A ^xx"
where A=u.Re/2. This we can solve analytically:
£--
-Ah
Ah


(64)
(65)
kP_2cosh(Ah) 2cosh(Ah)^W
which relates the values of the exact solution of (64) at
the equally-spaced mesh points P, E, and W depicted in
Figure 28. In general form we have
^P=^E^E+^W^W ^ 66')
where
-Ah
Ah
^E 2cos(Ah) an^ 2cos(Ah)


-46-
o-------oo
w \P E
Figure 28. Equally-spaced mesh points
For comparison, we can develop general forms for
central and upwind differences as well. Thus, for A^O, we
obtain
^P=Ve+(1_AE)^W (67)
where /
e-Ah
V 2cosh(Ah) for analytic scheme>
AR= for central difference scheme, and
for upwind difference scheme.
From Figure 29 we can clearly see that when Ah is small,
the central difference scheme is very close to the
analytic solution, but the accuracy degrades markedly when
Ah>l. The up-wind scheme is clearly better in this region,
but it is still not very accurate.
Spalding [31] suggested the so-called hybrid scheme,
whichvis of form (66) with
Ag=max(-2Ah, 1-Ah,0),
A^=max(2Ah,1+Ah,0),
Ap=Ag+A^. (68)
This amounts to the use of two linear segments AB and BC
(Figure 29) to approximate the analytic solution. Note,
however, the loss of accuracy in the region 0.5 < Ah < 2.
As an improvement, we developed a modified hybrid
V" 2+2Ah


scheme (MHS) which uses
A_=max{-2Ah, 1,5-5/3(Ah+0.5), 1-Ah,
0.5-l/3(Ah-0.5),0},
Aw=max{2Ah, 1.5+5/3.(Ah-0.5), 1+Ah,
W 0.5+l/3*(Ah+0.5),0},
AP=AE+AW. (69)
This amounts to using three segments AE, EC, and CD (Figure
29) to approximate the analytic solution. Note that it
closely follows the analytic solution everywhere.

Figure 29. Variety of difference schemes
4. Artificial viscosity revisited
First let us review some classical finite difference
theory. Finite differences are used to approximate deriva-
tives and truncation error is used as a measure the accu-
racy of the approximation. For example, to approximate
first-order derivatives, we have
/. ^E h
C = c c ^+,#*
*x h vxx 2
(forward difference),
or
^P
^E
2h
t h
+^xx#
h2
^xxx*6 +*##
(backward difference),
(central difference).
(70)
The forward and backward differences are first order,


-48-
and the central difference is second order. However, prac-
tical mesh sizes are not arbitrary small, so order esti-
mates may not be relevant for many applications.
To study practical cases of interest, we now introduce
artificial viscosity, which will focus on the vorticity
equation. Consider a general difference scheme for the 1-D
l
vorticity equation:
^P=AE^E+(1-AE^W (71)
*
where Ag is assumed to be a function of Ah. Define Ag so
that the analytic solution, £, of (64) satisfies
^p=aeV(1"ae)^w (72)
Now define the coefficient A* (locally) so that
Ag(Ah)=Ag(A*h), (73)
and
A*= F (AE)_1(AE(Ah))* (74)
We define:
(i) £ ~ oTr£ =0. the differential equation induced
X XX
by the difference scheme;
(ii) A*h, the mesh Re number induced by the difference
scheme; and
(iii) fij= 2^* 2^. the artificial viscosity parameter.
Although £ -rr- £ =0 is our target equation, the
effect of the difference scheme is to accurately solve the
perturbed equation
(75)
AA
This we write as
^x"[2A+(2A5r 2A)]^xx=0
^x 2A*'xx 0
(76)


-49-
or *x" 2^xx"^xx=0 (77)
This means that our difference formula has the effect of
adding artificial viscosity to the original 1-D vorticity
equation, is a measure of this additional term. In Table
10, we depict values of p/h=p^/Yi for the four difference
schemes over various values of Ah. Note that the upwind
scheme always exhibits positive artificial viscosity, while
the central scheme always has negative artificial viscosity.
The central scheme is better than upwinding when Ah but it loses accuracy completely when Ah>l. In contrast,
upwinding improves as Ah gets larger. Note that MHS has
an artificial viscosity that is small everywhere.
Table 10
Artificial Viscosity
Ah 0.0 0.25 0.5 1.0 2.0
* A*h li/h * r A*li li/h Af A*h a/* ^ e A*h m/* ^ f A*h li/h
Analytic 0.5 0.0 0.0 0.378 0.25 0.0 0.209 0.5 0.0 0.119 1.0 0.0 0.018 2.0 0.0
Modified Hybrid 0.5 0.0 0.0 0375 0.255 0.039 0.25 0.540 -0.080 , 0.121 0.05 0.015 0.0 OO 0.25
Central 0.5 0.0 0.0 0.375 0.255 0.030 0.25 0.549 -0.089 0.0 oo 0.5 -0.5
Unwind 0.5 0.0 0.0 0.400 0.203 0.466 0.333 0.34? 0.443 0.25 0.549 0.410 0.167 0.604 0.372
For completeness, we end this section by
estimating the artificial viscosity for central and upwind
differences for the 1-D vorticity equation. Using Taylor
series, we have
h^ h3 h4
^E^P+^xP*^+^xxP* 2T^xxxP* 3! +^xxxxP* 4! +
h2 h3 h4
V^P"^xP*h+^xxP# 2! ^xxxP* 3!1 ^xxxxP* TT*
Thus, the central difference formula for satisfies
.(78)


'* TTfS
2h ^xP vxxxP 3! ^xxxxxP 5!
For the 1-D vorticity equation this becomes
______ ^ 1
^E
2h
~^xP+[ 3!
2Ah (2Ah)
5j+]*h*exxp
(79)


-50
^xP ^c^xxP
(80)
where
"c r2Ah (2Ah)3
h[_3mh,,#]
2Ah -2Ah
=-(
e -e
-2Ah)/(2Ah)i
(81)
In a similar manner, for the upwind scheme (A>0) we
obtain
k*i-c
h cxP *u?xxP
/
where
=[e-2Ah-l+2Ah]/(2Ah)2.
(82)
In Figure 30, we graph terms in (81) and (82) as a function
of Ah. Note that the behavior observed in Table 10 is
verified here for all values of Ah.
Figure 30. Graph of error terms in (81) and (82)
5. Rotated hybrid scheme for the 2-D vorticity
equation
For the 1-D vorticity equation, MHS represents an im-
provement over the other approaches. In 2-D, MHS is applica-
ble only in the following situations: (a) when the Re number
is quite small or (b) when the flow direction is aligned


-51
with the grid line (x or y direction). However, when both'
Ah and Bh are large, MHS cannot be kept close to the analy-
tic solution. To overcome this, we introduce the following
rotated scheme:
A. In each subregion (Figure 31), we locally rotate the
coordinate axes to align the new x-axis (or y-axis) with
flow direction. The angle of rotation is
arctan(?)
A
arctan(^)+^
!*
!<
(83)
B. In the rotated coordinate system, we use
AP^PisAE1^E1+Aff1^W1+AN1^N1+AS1^S1
(84)
where
Ag =max{-2ARh^, 1.5-5/3(ARh^+0.5), l-ARh^,
0.5-l/3(ARh1-0.5), 0 },
Ay =max{-2ARhlf 1.5+5/3(ARh1~0.5), 1+Aj^,
and
0.5+l/3(A_h1+0.5), 0 },
It 1
A^ =max{-2BRh1P 1.5-5/3(BRh^+0.5), l-BRhlf
0.5-l/3(BRh1-0.5), 0 },
Ag -max {-2BRh1, 1.5+5 / 3 (B^-0.5), l+BRh1,
0.5+173(6^+0.5), 0 },
AR=Ac o sa+B* s i na,
Br=B* c osttAs ina,
fH/cosa
rsina
(0 £a£7r/4)
(7r/4 E
ft**1
O.H .CO
ma
cosa
(0 ^a^Tr/4)
(jr/4

-52-
C. To write (84) in terms of the original (unrota-
ted) system, we estimate £XT and £ by inter-
i Wi Ni Si
t
HE
polation: let b=l- ^ and define
^Ei=b^E+^1-b^N
^N1=b^N+(1-b^W)
£w =b-^w+(1-b)^S
and '
£g =b.^s+(l-b)£E (85)
We obtain
AP^P=AE^E+AW^W+AN^N+AS^S '
Ag=Ag *b+Ag (1-b),
Ajj=A^. b+AE *(l-b),
where


-53-
VV'b+A,
b+A^. (l-b),
and
6. Collective relaxation
Discretizations of (56) and (57) produce two
algebraic equations at each grid point. Sequential
relaxation of these individual equations would exhibit
a poor smoothing rate for multigrid: smoothing errors
in ^ could produce high-frequency error components in £
via the boundary conditions, and conversely. A better
approach consists of simultaneously solving these coupled
equations using collective relaxation. We first develop
the coupled equations at the interior points and sample
west wall boundary and southwest corner points. (Other
boundary points are analogous.) Note the special handling
of boundary points.
A. At interior points (Figure 32), we write
s
Figure 32. Interior point
B. At west wall boundary points (Figure 33), we have
N
W<^----O- -o E
P
f


-54-
1 4h -2A^ 4Ay ae£e+an^n+as^s+r£p
9 9 h4 Ap.h2 Afh2
4h Ap+2h Ay
-2h -2Ap 4Ap J
C. At the southwest corner (Figure 34), we have
^p '4h2 -2AW 2Ag 4Aw 4As
1 h4 Ap.h2 v2 Agh2
(4Ap+2Aw+Ag)h2 -2h2 -2Ap 4Ap+2A s 2As
ig. -2h2 -2Ap 2Aw

AE^E+AN^N+R^p
(89)
The basic idea behind pointwise collective
relaxation is to cycle through each point of the grid and
replace values of both unknowns there by values that


-55
satisfy the two associated equations at interior points;
this involves solving a 2X2 equation of the form (87). For
wall or corner points, this involves the 3X3 form (88) or
4X4 form (89).
Figure 34. Southwest corner )
7. Multigrid
Here we brifly describe the ingredients for the
FAS V-cycle multigrid scheme applied to cavity flow.
A full sequence of coarse grids (e.g., five for the 33x33
case) were used in our experiments, but for simplicity
we describe only the two-level case:
A. Collective relaxation on the fine grid equations:
L^p=0 ,
2
B. Collective relaxation on the coarse grid
equations:
T 2h*2h T2h,T2hch v T2hThth
Ll *P Ll Uh Jh LlfP\
_ L2h,2h L2h( j2h ,h j2h,p -J^d)
^ '2 uh Ah lJ?'2 W'
C. Coarse grid correction:
>h th xh (t2h T2h*h^
P<-P+I2h(fp-IhfP>'
<- Jrh+Ih (4-2h_I2h^h)
VP yP+i2h^P n V ,
(90)
(91)
(92)


-56
We use bilinear interpolation for and injection
for 1^ for both £ and A more sophicticated approach is
needed for residual restriction (see (10), (11) and (12)).
8. Alternating direction block-line Gauss-Seidel
relaxation
For low Re number flows, pointwise collective relaxa-
tion works well, but high Re numbersvflow tends to exhibit
anisotropies for which pointwise relaxation is no longer a
good smoother. Here we use alternating direction block-line
Gauss-Seidel relaxation. For example, the first step of a
west-to-east vertical-line sweep involves simultaneously
solving for ^ and £ at points of the first interior line and
£ at neighboring boundary points (primarily on the west wall).
See Figure 35. Except at such boundaries, this generally
involves solving a matrix equation with 2n+2 unknowns,
where n is the number of interior points in the
appropriate coordinate direction. Such an equation can
actually be solved inexpensively by multigrid, but since
it is essentially a 1-D problem, direct methods can be
efficiently used.
Figure 35. Block-line Gauss-Seidel relaxation
9. FAC for local refinement
To obtain details of the flow in the corners where


-57
secondary vortices appear [36], we use FAC for local
refinement. We have described FAC in some detail in Chapter
II and will discuss its application in cavity flow in the
next chapter. Special attention should be paid to the
residual at interface points. Here we use a finite volume
method to define interface equations.
10. Computational results
Starting with global grids, we computed four Re num-
ber flows with Re=0, 100,1000, 2000, using 33x33 and 65x65
meshes. The convergence histories are illustrated in Table
11. For low Re numbers, convergence is very fast. In fact,
the residual reduction factor of vorticity for each
multigrid V-cycle averages less than 0.05 at Re=0, and
about 0.1 at Re=100. For higher Re numbers, the convergence
speed slows, but it is still acceptable. The residual
reduction factor is, on the average, about 0.25 at Re=
1000.0 and 0.33 at Re=2000.0 for each V-cycle. Here we call
one alternating block-line Gauss-Seidel relaxation on the
finest grid one work unit. When considering total work
i
units for high Re numbers, one should take into account
the fact that typical relaxation methods presently in use
require four or more units per step. The maximum stream
function values and their position are listed in Table 12,
which shows that the computational results are quite
dependable (cf. [34][35]).
We use FAC for local refinement with a 17x17 mesh in
both the northwest and northeast corners of a 33X33 global
mesh. Table 13 gives^ the convergence history showing the
efficiency. Note that the average composite grid residual
reduction factors are less than 0.05 for Re=0, 0.25 for
Re=100, and 0.5 for Re=1000.


-58-
Table 11
Convergence History of FAS
Mesh 33x33 Mesh 65x65
Re
l 0.0 100.0 | 1000.0 | OO | 100.0 \ 1000.0 | 2000.0
Cycle WUs Residual
1 7.031 6.47E-2 4.33E-2 9.25E-2 1.18E-1 8.50E-2 4.24E-2 128E-1
2 14.063 1.09E-3 4.89E-3 1.05E-2 1.53E-3 3.53E-3 1.82E-2 5.31E-2
3 , 21.094 6.72E-5 5.12E-3 1.81E-3 1.89E-4 4.89E-3 4.17E-3 1.01E-2
4 28.125 3.23E-6 5.83E-5 3.69E-4V 6.73E-6 5.94E-5 1.65E-3 6.02E-3
5 35.156 4.18E-8 6.66E-6 1.08E-4 2.04E-7 S.54E-6 2.84E-4 1.61E-3
6 42.187 4.45E-9 8.23E-7 2.72E-5 1.67E-8 6.92E-7 7.82E-5 3.52E-4
Table 12
Computational Results
Mesh 33x33 65x65
Re 0.0 100.0 1000.0 0^ 100.0 1000.0 2000.0
9.96E-2 1.05E-1 9.64E-2 1.00E-1 1.08E-1 1.22E-1 1.13E-1
Cctntr 3.02 3.23 1.82 3.02 3.26 1.95 1.72
l*c.yc) 0.50.0.25 0.62.0.25 0.55.0.43 0.50.0.24 0.63.0.26 0.55.0.43 0.54.0.47
Table 13
Convergence History of FAC
Re = 0.0 Re = 100.0 Re"= 1000.0
x Residual
Cycle WUs Coarse Composite Coarse Composite Coarse Composite
1 7.031 6.47E-2 6.49E-2 4.16E-2 4.75E-2 6.03E-2 2.71E-1
2 14.063 1.09E-3 1.10E-3 3.02E-3 4.42E-3 3.46E-3 5.14E-2
3 21.094 6.71 E-5 6.74E-5 2.26E-4 9.06E-4 8.03E-4 1.45E-2
4 28.125 3.23E-6 3.28E-6 8.58E-5 2.69E-4 4.44E-4 6.02E-3
5 35.156 4.15E-8 4.19E-8 1.65E-5 9.79E-5 2.14E-4 1.45E-3
6 42.187 4.39E-9 4.43E-9 6.87E-6 3.11E-5 1.26E-4 9.516E-4


-59-
CHAPTER VI
FVE FOR DRIVEN CAVITY FLOW
The conventional finite volume method can be a very
effective approach to discretization of fluid flow
equations in conservative form. However, the choice of
j points and formulas to approximate the necessary fluxes is
a little ad hoc. FVE, which was introduced in Chapter II,
provides a more systematic and often more accurate form of
the finite volume method. The purpose of this chapter is
to study FVE in the context of cavity flow, especially for
i r
high Re numbers.
1. Equations and boundary conditions
We again consider cavity flow:
V.-W fc^yy*'0
$ +£=0
Txx Tyy s
^(0,y)=iKl,y)#x(0,y)=^x(l,y)=0
^(x,0)#(x,l)=^ (l,y)=0
J
*y(0,y)=l (93)
Here, £ is vorticity and ^ is stream function.
2. Difference scheme r
We first develop an FVE scheme designed for low Re
number flows. In this approach we use triangular elements
and rectangular volumes (Figure 36).


-60
OB OB BOBO oaoa bo bo oaoa BOBO OBOBOBO BOBOBOB OBOBOBO BOBOBOB OBOBOBO BOBOBOB aoaoa OBOBO aoaoa OBOBO aoaoa OBOBO
OB OB BOBO OB OB OBOBOBO BOBk'B^B OflOBOBO aoaoa OBOBO aoaoa
BOBO OB OB BOBO OB OB BOBO 08 OB BO BO BOBOBOB OBOBOBO BO BO BO B obobobo BOBOBOB obobobo BOBOBOB OBOBO aoaoa OBOBO aoaoa OBOBO aoaoa OBOBO
(.) u_10 d.o) x
-----------
Figure 36. Elements and volumes for planar cavity
flow
s
Integration of equation (93) over a control volume ft.
with surface P. yields
1 i
[ n W' ds +f £ da=0 . (94)
Jr. Jn.
1 1
A. Interior point (Figure 37): Usinga piecewise
linear basis function for £ and ^ and writing
^=ax+by+c
^=ax+/3y+')1, (95)
then in element I
fl = [ [b(/3y+7)
e2 Jr.
l
fe'W]- <">
v
Formulas for the other surface fluxes are similarly
derived. Also, it is fairly straightforward to use (95) in


-61
evaluating the volume integral of $ in (94). We are thus
led to a difference formula
AP^P=AE^e+AW^W+AN^N+AS^S+ANW^NW4;ASE^SE
4^p=^E+^w+^N+^S+^h2 * (97)
where
Re + 8(^'P+^SE_2^N)
AW= Re + 8('J'P+'I'NW_2'J'S)
V Re + 8(2^E_vV"2lW
AS= Re + 8(2^W_^,P2^SE)
ANW=^N-^W
ASE=^S"^E
Ap= AE+AW+AN+AS+ANW+ASE
^=12(^E+^W+^N+^S42^NW42^SE+7^P) /
S SE
Figure 37 .v Sample interior point
B. Sample west boundary point (Figure 38): Because of
the Dirichlet boundary condition on only one equation
is used on the volume located at the boundary. The natural


-62-
choice is
f n ds +f £ da=0
Jr. Jn.
for which we have
fl =fl =fl =0
w n s
Jn £ dxdy-48('4^E+3^N+^S+2^SE+14^P^'
i
This yields the difference equation
^vr ~ Te(4^T,+3^|J+^<;,+2^QI,+14^p) .
(98)
(99)
Figure 38.^Sample west boundary point
C. Sample southwest corner point (Figure 39): Here we
have
fl =fl =fl =0
e w s
ns=-^y !-!
This yields the formula
16(^E+^N+2^P) = 2 *
(100)


-63-
i
Figure 39. Sample southwest corner point
3. Multigrid algorithm
The basic MG algorithm is similar to that for the
Poisson problem, except that it uses collective relaxation.
Near boundary points, the equations are grouped and solved
as depicted in Figure 40. The asymptotic factors of
convergence we achieved for tests with the FAS V(2,2)-cycle
on 33x33 grids was about 0.0783 for Re=0, 0.0752 for Re=50,
and 0.142 for Re=100.
Figure 40. Boundary point groups for relaxation


-64-
I
1
\
4. FVE accuracy
Here we used individual discrete energy norms for £
and \p. Specifically, write the difference equations for
(93) as
. <10
where A^CV^) indicates the nonlinear dependence of the
y ^ ll ^ ll
first equation on Let £^=1^ £^ and ^ where
no no no no
1
h^=g^- and £^ an^ ^ are the exact solutions of (100) with
*o Ko
L.
h=hQ (computed after many V-cycles). We define -the
discrete £ and $ energy norms errors by
lvhMI = "WVV2
and
11 lwh< 11 = * v
respectively, where and w^=^-^. Table 14
depicts the results for various mesh sizes. Note that the
2
results suggest only 0(1) accuracy in vorticity, but 0(h )
accuracy in stream function, which is the physical variable
of interest.

Table 14
FVE Discretization Errors
Mesh /' 9x9 17x17 33x33
Re=0 Vorticity errors 0.788e0 0.794e0 0.526e0
Stream fn. errors 0.244e-2 0.642e-3 0.131e-3
Re-50 Vorticity errors '0.931e0 0.845e0 0.538e0
Stream fn. errors 0.257e-2 0.662e-3 0.134e-3
Re-100 Vorticity errors 0.134e+1 0.961 eO 0.566e0
Stream .fn. errors 0.357e-2 0.924e-3 0.179e-3
i


5. Nested iteration
Table 15 demonstrates that nested iteration again
achieves its goal as a semi-direct solver at a cost of about
five work units.
Table 15
Accuracy of Nested Iteration for Stream Function
Mesh 9x9 17x17, 33x33
Re=0 Discretization errors 0.244e-2 0;642e-3 0.131e-3
Actual errors 0.249e-2 0.663e-3 0.138e-3
Re=50 Discretization errors 0.257e-2 0.662e-3 0.134e-3
Actual errors 0.263e-2 0.697e-3 0.145e-3
Re=100 Discretization errors 0.357e-2 0.924e-3 0.179e-3
Actual errors i 0399e-2 0.110e-2 0.261 e-3
6. FAC
In practice, it might be appropriate to place levels
of finer patches in both the NE and NW corners of (1 to
resolve the vortices there (Figure 41). However* here we
consider only the two-level composite grid depicted in
Figure 42 for simplicity.
Figure 41. Composite grid for planar cavity flow


-66-
i'
o
X

Boundary point
Global grid point
Patch point
Slave point
Interface point
Finite element
Finite Volume
Figure 42. Elements for simplified composite grid
A. For interface points (Figure 43)
Illustrating the approach for a west side
interface point for the vorticity equation at P in Figure
43, we have
n tt-
V*S,1
2h (2^SE+8^P+8^S)+Rr
1. 1 ^SE~2^P ~2?i
!]
hr %^SE,U 1> x 1
h 2^P+?E+4^SE+Re^
^P
h
hr ^NE ^E,l,* 1*1 + ^ 1
R
i o i i ^NEl^P "2^:
Ir rNE rE,l> 1* 1* \ 1 E^p
lL h (t^P'fi,E+4^NE)+Re* h
.If )
L 2h '2^NE o^P o^N;
+s'-
(102)
Re h
Similar formulas for the other fluxes yield the equation
Ap6p=AE^E+AW^W+AN^N+AS^S+ANW^Nff+ANE^NE+ASE^SE
(103)
t. where
AE_8(^SE-'W4Re


AN_8(2^'NE'^'P ^'NW)
AS=8( 2 V^SE^P')
ANW=8(^'N"^W)
ANE=8(21|!'P^'N+^E) +2Re
(2^S'^'P_^E)+2Re
AP^(^N^S+V^NE^SE)fRe *
The difference scheme for the stream function equation is
^p=h+h^^s^Â¥m^Â¥sE+^ (104)
where
£ ^ ^ +_3_/4L*
5 96^E i2*W 1536^N 512*S 24^NW 324VNE 324^SE 1536^P.
Figure 43. West side interface point
B. FAC and AFAC
Convergence history of the FAC and AFAC algorithms
are depicted in Table 16 for the composite grid, where the
coarse grid and patch are both 33x33 grids.


-68-
Table 16
Residuals on a 33x33 Based on a Two-Level Composite Grid
Re 0 50 100
FAC AFAC FAC AFAC FAC AFAC
Cycle 1 0.61 Oe-2 0.696e-2 0.622e-2 0.709e-2 0.941 e-2 0.942E-2
Cycle 2 0.676e-3 0.963e-2 0.593e-3 0.937e-3 0.204e-2 0.190E-2
Cycle 3 0.524e-4 0.138e-2 0.486e-4 0.130e-2 0.383e-3 0.796E-3
Cycle 4 0.700e-5 0.368e-3 0.695e-5 0.389e-3 0.518e-4 0.427E-3
Cycle 5 0.114e-5 0.169e-3 0.111e-5 0.134e-3 0.998e-5 0.104E-3
Cycle 6 0.198e-6 0.798e-4 0.196e-6 0.781 e-4 0.197e-5 0.946E-4
Factor 0.127 0.409 0.125 0.406 0.184 0.398
7. Medium and high Re flow
To apply FVE to the case of medium and high Reynolds
number flow, we replace the piecewise linear elements by
the element functions that have a form to better fit the
local character of the flow. In particular, in each trian-
gular element we assume that
f=e^x+By(ax+by+c), f=ax+/3y+'y, (105)
where A=Reu/2 and B=Re*v/2 were determined by local values
that characterize the flow direction. The convergence his-
tory of multigrid for the 65X65 grid is depicted in Table
17 and Table 18, which compare our results (FVE) with those
reported in (16] It shows that our 65x65 grids produce
estimates that compare favorably with the 129x129 estimates
obtained in [16].


-69-
Table 17
Multigrid Convergence History for a 65X65 Grid
RE 400.0 1000.0
FAS cycle Work unit Rf Factor Factor
1 5.64 0.472E0 0.924E0
2 11.28 0.622E-1 0.132 0.296E0 0.321
3 16.92 0.141e-1 0.226 0.338e-1 0.114
4 22.59 0.224e-2 0.173 0.164e-1 0.485
5 28.20 0.520e-3 0.214 0.105e-1 0.637
6 33.84 0.127e-3 0.244 0.493e-2 0.471
7 39.48 0.443e-4 0.349 0.236e-2 0.479
Avg. Factor 0.213 0.369
Table 18
Comparison of FVE Results with Those in [16]
RE 400.0 1000.0
Mesh ihnax 6c Xc Yc '{'max 6c Xc Yc
Ghia 129x129 .113 2.29 .555 .395 .117 2.05 .531 .438
Ours 65x65 .106 2.10 .563 .391 .095 1.83 .538 .422


-70-
CHAPTER VII
MG FOR PSEUDO TIME-DEPENDENT FLOWS
It is common in practice to use time-dependent equa-
tions for solving time-independent problems. For example,
in computational fluid dynamics, unsteady Euler equations
are often used for steady flows. The advantages of this
pseudo time-dependent or pseudo unsteady approach are that
the unsteady Euler equations are hyperbolic even in tran-
sonic regions, explicit time-marching schemes can be easily
employed and the approach usually ensures physically correct
solutions., .
The purpose of. this chapter is to explore the equiva-
lence between time-marching schemes and relaxations on the
steady-state equations. This is admitedly a rather simple
and commonly understood notion. However, our purpose here
is to show how it can provide a simple vehicle for extending
MG to general time-dependent equations. One of the under-
lying motives for this work is to counter the misconception
that MG does not apply to hyperbolic^problems. 1
1. Explicit time-marching and Jacobi iteration
Consider, first, the simplified heat equation as an
example:
Assuming the mesh size Ax=Ay=h, we may write the forward
Euler difference equation for (106) as
v
(106)
T^-T? .
_i_J_LlI
At
(107)
or
I


-71-
For convenience, we rewrite (108) as
T?+*=(l-4At.h_2)T? .
i.J
+4Ath-2-j(T? .+T? .+T? ,+T? -) .
4 l+l,j l-l,j i,j+l 1,3-1
(109)
_2
Note that 0£4At*h £2 is a necessary and sufficient condi-
tion for stability.
Now the steady equation for (106) is
T +T =0 (110)
xx yy
and its corresponding difference equation is
T. .-2T. .+T. T. ,-2T. .+T. ,
1+1,3 1.3 1-1.3 + 1,3+1 1,3 1,3-1 _0 (111)
h2 h2
The Jacobi relaxation scheme for (111) is
T?+*=j(T? .+T? .+T? .+T? .)
i,3 4 1+1,3 1-1,3 1,3+1 i.J-l
(112)
If we include a relaxation parameter w, then the (damped)
Jacobi scheme can be written as
T?+*=(l-u>)T? .+wrj(T1? .+T? .+T? .+T1? .). (113)
1,3 i,3 4 1+1,3 i-l,3 i.J+l i.J-l
Note that the Jacobi relaxation scheme (113) is just the
explicit time-marching scheme (108) with w=4At/h This
equivalence has one important practical consequence: since
Jacobi iteration by itself is not an efficient method for
solving Mil), and especially since w=4At/h is not an
optimal relaxation factor, the explicit time-marching
scheme by itself cannot.be a very efficient method for
solving the steady equations.
2. Explicit time-marching and Gauss-Seidel iteration
In this section we illustrate that Gauss-Seidel on
steady equations can be interpreted as a time-marching
\


-72-
scheme for a certain time-dependent equation. The Gauss-
Seidel iteration scheme for the steady heat equation (110)
may be written as
T? .-2T?+*+T?+J T? 1-2T?+*+T?+l
i+1,3 1,-3_1-1,3 + 1,3+1 i,J______LtJzI =0
-9 -9
We rewrite (114) as
qjl+l^ipll jpll+l q^Il+1 ^rpH
4,_LJ_____Li hi|j i-i,j i,3-i .i.j-i
h2 h2 h2
T? .-2T? .+T? T? .-2T? .+T? .
_ l+l.l l.l l-l.l l.l+l i.l l.l-l
h2 h2 '
^ 2
Letting d=At/jh then (115) becomes
rpH+1 rrh n^l+1 rr^l mH+1 rpH
d.r4. 1.3- 1,3 _ 1-1.3 1-1,3 i.Jrl 1,1-1]
d,L3 4 At At At J
T? .-2T\ .+T? T? ,-2T? .+T? ,
= 1+113____1,3 1-1,3 . 1,3+1 1,3___1,3-1
_ 9 _ 9
(114)
(115)
(116)
Equation (116) is equivalent to the usual forward Euler timey
marching scheme for the unsteady equation
2d.Tt+dh.(Ttx+Tty)=TM+Tyy . (117)
This observation suggests that since Gauss-Seidel is
generally more efficient for solving (111) than is Jacobi,
Equation (117) may provide a better vehicle for pseudo
time-dependent problems than equation (106). This is a topic
for further exploration.
3. ADI time-marching and line relaxation
We now show how the Alternate Direction Implicit (ADI)
method can be interpreted as a relaxation scheme on steady-
state equations. ADI applied to (106) may be written as
i '
I


-73
T?+*-T? T?+J -2T?+^+T?+? T? ,-2f} .+T1? ,
1.3 i ,1 _ 1+1,3 i,3 1-1,3 _1,J+1__LzJ___i.J"1
At ,2 ,2
Tn+2_Tn+l ji+1 _2Tn+l+Tn+l ^jn+2 _2t?+2+t?+2
1.3 i,3 _ i+1,i i,i i-1,i i.i+1 i,j i.i 1
At ^2 ,2
h h
(118)
2
Letting d=At/h we can write the first equation as
T?+J .-(2+ J)T?+*+T?+J .s-It11 .-(2-.+T? .1.
1+1,3 d i,J 1-1,3 L 1,3+1 d' 1,3 1,3-lJ
: - (119)
Note that for d=At/h2=l/2, (119) becomes
T?+J .-4T?+*+Tri+J .=-[t? ,+T? J (120)
i+l,J i,3 1-1,3 L i,J+l i,3-lJ ,
which is just Jacobi line relaxation applied to the steady-
state heat equation (110).
4. Generalization
Consider a time-dependent problem
LtT=LT-f . - (121)
We may discretize the time derivative term as
(LtT)p -1 k^1+S VE (122)
k k
and the right-hand side of (121) as
(LT-f)p l AkT£+1+$; B^T^-fp . (123)
k k
Here, k ranges over the indices of point P and its appro-
priate neighbors and the superscript n marks the current
time step. Substituting terms in (121) with (122) and (123),
we have
l VF+X DkT£=S Ak^+1+S B"Tj-fp . (124)
k k k k
Putting all terms at time step n+1 on the left-hand side and
all other terms on the right-hand side, we get


-74
l cX+1-Z AkT^+1=I D^-fp . (125)
k k k k
It should be apparent that (125) .can be interpreted as a
relaxation scheme for solving a discrete version of the
steady-state equation LT=f. (125) can be used as a basis for
developing a multigrid scheme for solving the steady-state
equation. This then immediately translates to a multigrid
scheme for solving (121). To see how this is done, we
examine the case for Euler equations.
5. MG as an explicit Euler solver
Here we will use our observation of equivalence between
steady-state relaxation and pseudo time-marching scheme to
\
develop an MG solver. We take flow in a channel with a
circular arc "bump" on the lower wall as an example [39].
(See Figure 44.)
"10 0 1.0
Figure 44. Flow in a channel with a "bump
A. Governing equations
The unsteady Euler equation may be written in the
conservative form
where
(126)


-75-
q=
p ' pu 2
pu , E(q)= pu +p
pv puv
. pe . . pu(e+ £).
and
F(q)=
PV
pvu
2
PV +p
pv(e+
For homoenergetic steady flow, we have
hQ- p + 2^ (u2+v2) = const . (127)
Thus, the energy equation can be eliminated from the
solution procedure and the pressure is then determined by
the algebraic relation shown in (127).
B. Discretization
An explicit Lax-wendroff scheme can be expressed as
(see Figure 45)
^q1=q"+1-q"=(5q1)A+(^q1)B+(^q1)c+(5q1)D (128)
Here,
<{Vc= i(4<>c-!Hf tK
D i(iv &V
where all of those terms are meant to be calculated at time
n and are evaluated according to the following:
rE.+E_ E_+E. F.+F.
Aqc=[~V^y V^y+ ~L2Ax -
F+F
2 3
Ax
] At
JAxAy
AE=
dE
Wc^C
'uc=(f|)c4V
etc.
2


-76-
Here, E. is evaluated at point i and the Jacobian
,dEv lr/-5E\ /^E-. ,dE,. ,,3E-, i
(W)e=4t(55)i+W2','(5?)3+(3j)41
The condition for stability is
Ax Ay

At ^ min (
|u|+a |v,|+a
where a is the local sound speed.
),
(129)
1 2 3 \
B / fjgg V \ . k Ay
A >1 'i D h Ax

Figure 45. Difference scheme points
For transonic and supersonic flows, artificial viscous
damping should be introduced into the solution procedure
[11]. For example, the formula for (^Q^q may become
<*>,>0- ![^c- iH-
(130)
where
,At OOx
VlVWV'
At,
Ay' k
and a is an artificial damping factor taking a value between
0.0 and 0.1.
C. Multigrid
The time-marching Scheme (128) is written
qi+1=V(5qi)A+(5qi)B+(5qi)C+(5qi)D
which is equivalent to a Jacobi relaxation for the
stationary problem
(131)

(132)


-77-
We can easily design a MG scheme for the pseudo unsteady
Euler equation as follows:
Step 1: Use (131) to do time marching in the fine grid.
Step 2: Use (i32) to calculate the residual in the fine grid
and use full weighting stencil to transfer the
residual from the fine grid to the coarse grid.
Step 3: Do time marching in the coarse grid by using the
following formula:
CQ+1)=(^)+Ljh(Q^)-Ljh(lf,;)+Tf(Lhq^). (133)
Step 4: Make the correction
q q+I2h(Q-Ihhq)'
Here, q and Q are solutions of the fine and coarse grid
equations, respectively, and
L2hV<*>i V^Vb^Vc^i V
D. Computational results
Here we report some illustrative results for transonic
flbw (M =0.675). The isomach lines are depicted in Figure
' 00
46, and the mach number distribution on the upper and lower
wall is shown in Figure 47. They both coincide with.the
results reported in [11], Figure 48 shows the multigrid
convergence speed which is much faster than the single grid
relaxation scheme. However, that is still not as fast as we
can expect: a scheme that uses the more efficient block
Gauss-Seidel scheme or other more efficient implicit time-
marching schemes should make significant improvement in the
convergence rate.
j


I
-78-
Figure 46. Isomach lines
Figure 48. Convergence history


I
-79-
CHAPTER VIII
SOURSE TERM CORRECTION (STC) FOR
INCOMPRESSIBLE NAVIER-STOKES EQUATIONS
i
There are two common formulations of the Navier-Stokes
equations, vorticity-stream function (£,^) and primitive
variables (u, v, p, ...). It is usually simpler to use the
equations for 2-D but^the extension of this approach
to 3-D is not so advantageous. It is also problematic to
determine boundary conditions for vorticity. On the other
hand, primitive variable Navier-Stokes equations are easily
extended to 3-D. However, there are several difficulties
(cf. [27]) encountered in the solution of incompressible
Navier-Stokes equations in primitive form. First, central
differences for the pressure term can lead to ill-posedness
of the discrete equation. Second, there are no explicit
boundary conditions for pressure. Third^ simple explicit
time-marching schemes are usually inappropriate since they
cannot ensure VV=0 at each time step. The first two
difficulties can be avoided by using staggered grids [37],
the third by using artificial compressibility [9]. However,
staggered grids complicate programming and require specifi-
cation of an image boundary velocity, and artificial
compressibility represents a departure from the physics.
The motivation for our development of the source term
correction (STC) method is to eliminate these troubles. As
we intend to show, the STC method is simpler than the vor-
ticity-stream function approach in 2-D and it easily extends
to 3-D. We first take 2-D planar cavity flow as an example.
(
1. Reformulating the governing equations and boundary
i
conditions.
Using primitive variables, we have the following
equations and boundary conditions for cavity flow (Figure


-80
49):
V*V =0 ,
da du. 1 /d2u d2Uv dp
Udx +Vdy ~ + dy*' ~~
dv dv 1 ,d2v d2v~, dj>
l^x +vdy Re ^2 + 3y2 __ dy
u(0,y)=u(l,y)=u(x,l)=0 ,
v(0,y)=v(l,y)=v(x,l)=v(x,0)=0
u(x,0)=l .
(134)
We have three equations for three unknowns, u, v and p.
Figure 49. Planar cavity flow
To formulate the STC method, we introduce a simple
reformulation of these equations. In particular, we keep
the boundary conditions and two momentum equations, but
replace VV=0 by the following:

fe (V.v)=0 ,
d ;n-i>
in 0
v.v =0
on r.
(135)
Here, Fis the boundary of fl. Note that (135) is equivalent
to V.V=0. However, use of (135) allows us to avoid discre-


-81
tization of pressure terms, p and p as we shall see.
x y
2. Discretization
For simplicity, consider a uniform grid with Ax=Ay=h. ,
Ignoring boundary conditions, we write our formulation of
the Navier-Stokes system as
k (V-'?)=
k (V*'f)=o
2Aft!+ 2B§7
a2 a2
,d u a u x
(-T----2 *
dx2 dy2
,d2v 32v x
Su
Sv
(136)
where A=u*Re/2, B=v*Re/2, S =-Re^ and S =-Re^ Now we
u ox v ay
have four equations for'four unknowns, u, v, Su and S .
Considering first the two momentum equations and
using the standard five point scheme, we have the general
form of our difference schemes given by
Apvp_AEvE+AWvW+ANvN+ASVS+Sv - ^137')
If A,B, Su, and Sy are computed exactly, then the solution
of (137) should exactly satisfy ^ (V*^)=0 and ^ (V^)=0.
The problem is that S and S are unknown. The basic idea
U. V /
behind the STC method is that the current approximations
to Su and Sy are first used to compute the solution, un+*
and vn+1, of (137). These estimates for S and S are then
u v
corrected in an attempt to force subsequent iterations to


-82
satisfy the continuity equation.
3. The STC method
Assume for exposition that the intial guesses for u, v
and Sy are the exact solutions of (137), i.e., un=u, vn=v
and s=Sy throughout ft, where u, v and Sy are the exact so-
lutions of (137). Assume further that the initial guess for
S is also exact everywhere except at point P. Let Sn=S +AS .
u u u u
Since AS?K), then substitution of s into (137) leads to an
inexact approximate solution to u, u114^ u. Let Up+*=Up+Aup.
From (137) we have that {
' UP+1= T5(4^E+AuW*ANVASttS+Su+4Su) ' (138)
Ap
But
p= Tu^AEuE+AWuW+ANuN+ASuS+Su^ '
Ap
Subtracting (139) and (138), we then have
n+1 a 1 in
Up -Up=Aup= ASU .
h
(139)
(140)
Now approximation of the first equation in (136) yields
(see Figure 50)
Sx V; 2h
=0
or

UEE UP VNE VSE, ,UP UWW VNW_VSWsl
+ ==> 2h + -----------2h->r0- (141)
2h
2h
2h
n+1
Since Up ? Up, we cannot expect\(l41) to be satisfied by


-83-
un+* In fact,
t
dx
(Vh.^
n+l _ n+1 _
1r,UEE UP VNE VSE^ UWW NW VSW
* sc <----= +
2h
2AUp
AS..
u
4h
Note that
2h2A£
2h
* 0
2h
2h
>]
(142)
ASu=-2A^h2*^- (Vh^) .
This holds for interior and boundary points.
(143)
Figure 50. Difference scheme points for continuity equation
(143) shows that an incorrect un+1 implies that
^h ^ .
fa (V V)p7H). However, (143) also shows how to estimate the
deviation, AS of Sn from the true S Thus, the basic
idea behind STC is to choose
S£+1=S" ASu =S"+2A^.h2.^ (VM) (144)
for use in the next iteration. (144) is the source term
correction formula in the x-direction. We similarly define


-84-
the correction in y-direction:
sn+1=sn
V V
- AS.

dy
(145)
4. Iteration scheme and STC algorithm
The basic iteration scheme for STC may be written as
n+1 1 ,.u n .u n+1 Au n .u n+1 0n>
UP .u (AEUE+AWUW +ANUN+AS S Su)
vn+1_ I_ (AXvn+Alvn+1+A.vvn+AVvn+1+Sn)
P .v lTE TW +Yn+as s V
*P
n+1 n+1 ,, \ n
VP =VVP +(1-w2)vP
(146)
(147)
(148)
(149)
sJ+1=s + W4.h.[(vh.?)N-(vh.'?)g] ,
(
where u w w and w, are relaxation parameters,
l 2 3 4
that we must have an additional condition,
d(S ) d(S )
u v
(150)
(151)
Note
(152)
ay ~ ox
for (136). However, (152) could be automatically satisfied
if we choose w =w in (150) and (151), because we have
3 4
n d /r7h

i=l
n

and
b lb (7M)i >=s 4 "*>,>
(153)
!


-85-
t
J
We also have
AP=AE+AW+AN+AS (154)
where Ag ,A^ ,A^ and Ag are determined by difference
schemes left to be chosen. For example, the upwind scheme
is defined by
Ag= Max (l-2Ah,l) ,
Ay= Max (l+2Ah,l) ,
Ajj= Max (l-2Bh,l) ,
Ag= Max (l+2Bh,l) . (155)
The algorithm for STC can then be described as fol-
lows:
Step 1: Choose an initial guess for u, v, Su and Sy. For
. 0 0 o0 nO n
example, u =v =S =S =0.
Step 2:' Perform the basic iteration step (146)-(149) 3 or
4 times to obtain new un+^ and vn+*.
Step 3: Use (150) and (151) to compute new S^+* and s+1.
Step 4: Stop if
and
max i n+1 lup N pi = ei
i n+1 n |
max Ivp -Vpl < £, .
M E E l(Vh.^)?+*|< e
=1 j =1 l.J
(156)
Otherwise, go to step 2.
Here, M and N are the maximum index of grid points in
the x- and y- directions, respectively. 5
5. Solution for pressure
After obtaining the solution for, velocity u and v, we


-86
may solve a discrete Poisson equation for pressure based
on
7M( k ,(w > < k)(k >] <157>
For details see [30].
6. Other flow problems
A. Axially-symmetric Navier-Stokes equations
The governing equations for axially-symmetric flow
may be written as
k k (V^>=
(158)
dv
dv
x x 1
var +V37 r¥
d2v d2v
dv .dv ,
r r 1,
V -5- +V #-3------5(
x ox r dr Re
dx*
d2v
dri
1 ^X' dp
+ )= T*
r dr
dv
a2
d V , i/v V
____r l_____r _JL\_ dp
2 + r dr 2' dr
dr'
(159)
(159) can be rewritten as
dv dv d2v a2 d V
2A#3xX +2B.^:*- (r dx2 X 4 dr2 >=Su
dv dv a2 d v a v V
2A*s +2B*^j-1 (a/ )= S -7 V 2 r (160)
which are similar to the form in Cartesian coordinates.
Note here that
2A=Rev 2B=Re*(v ---- ) and
x r r
dv dv v
n 3 __x ____r r
VV = 3 + 3 +
dx dr r .
B. General curvature coordinates
Taking the u-equation in (136) as an example, in x-y
coordinates we have


-87
2lM +2B&
dx. dy ^2 Qyi u
(161)
dx^
In the general curvilinear coordinates
£ = £ (x,y)
n = T] (x,y),
(161) becomes
(2A£ +2B£ -£ -£ )& +(2Aj? +2B7 -t) -tj )§£
^x ^y vxx vyyo£ 'x 'y 'xx 'yy arj

]=su+(2£xV2eyy!^
(162)
(163)
(162) may be written in the form
,*du *du ,*d2u t*32u % * '
2A qt- +2B -3- -(C - + D - ) =S
d( *> d(2 ar,1
where the starred quantities are defined accordingly.
Similar representation holds for the v-equation in (136)
Also,the continuity equation can be written as
§£ (V.ib =0
% (V,^) = 0
(164)
Note that (163) and (164) have the same form as the planar
problem .
C. 3-D case
The governing equations in 3-D are written as follows:
L (V.V) = 0 ,
^ (V.?) -o',
-0 .
2k% +2^ +2cfi (! + Sfu
dx dy dz gyZ
*2
% -s
dz2 u


-88-
/
2As+2B|f +2C% -
+< £ -
a2
rO V
a2
o v
dz dy
a2
§4) =s
az2 v
,3 * 9 W
(2 + ~T +
dx2 dy dz
£24 =s
2 W
(165)
The extension of the STC algorithm to 3-D is straight-
forward.
D. Turbulent flows
Following Boussinesqs viscosity concept [40], we can
r
think of the effective viscosity parameter |i as a variable
in the flow field which is determined by some turbulent
model equations. Focusing again on the u-equation in
(136), we think of the Reynolds number as a variable:
9u du [9 1' du 1 d2u d ,_ls du 1 d2ul
udx +Vdy |dx ReJ*dx + Re ^2 "^cfy^Re^dy ^Re ^2j
= S* . (166)
(166) may be written as
SI < £-
.du.
d x d y
u
where
(167)
2A=ite- fe 2B=Re.[v- ^ C^)] . _
(167) has the same form as that for laminar flow .
7. Preliminary computational results
Example 1. Planar cavity flow
We set Re=100, e =e=e.,=0.005, w =w =0.5 and
12 3 12
w3=w4=0;4. The computational results for this case are


-89-
depicted in Table 19. It shows that our results compare
favorably with Mills experiment [34] which measured the
position of vortex center, xc and y in a closed 1
planar cavity with a moving upper wall. Table 19 displays
the position of £ as well as ^ , the maximum value of 4'-
c max
Table 19
Computational Results for Cavity Flow
Mesh Inner iters. Outer iters. 4f rmax xc yc
11x11 5 11 ,0.0921 0.58, 0.27
17x17 5 21 0.105 0.59, 0.26,
Mills' experiment 0.66, 0.25
Example 2. Suddenly expanded pipe flow
Because STC does not discretize pressure derivatives
and does not use staggered grids, programming is quite
simple. Also, convergence is very fast as shown by the
results in Table 20: compared with the well-known TEACH-T
code [38], STC saves more than half of the work units and
about 3/4 of the cpu time.
Table 20
Comparision of Vork Units with TEACH-T Code
Code Mesh Inner iters. Outer iters. Work units
TEACH-T 16x16 3 112 336
STC 17x17 5 32 160
Although we have developed the STC algorithm and
presented some preliminary computational results, we are
really only at the begining of this development. While
the results are encouraging, substantially more compu-
tation and several theoretical issues must be pursued.


-90
REFERENCES
[1] J. F. Thompson, Z. U. A. Warsi and C. W. Mastin,
Numerical Grid Generation, Foundamentals and
Applications, North-Holland (1985).
[2] A. Brandt, "Multi-level adaptive solutions to
boundary-value problems," Mathematics of Computation,
31, pp 333-390 (1977).
[3] A. Brandt, Multigrid Techniques: 1984 Guide with
Application to Fluid Dynamics. GMD studie.GMD, St.
Augustine (1984).
[4] R.W. McCormack, "The effect of viscosity in hyper
velocity impact cratering," AIAA paper, pp69-354 (1969).
[5] E. M. Murman and J. D. Cole, "Calculation of plane
steady transonic flows," AIAA J., 9, ppll4 (1970).
[6] A. Harten, "High resolution schemes for hyperbolic con-
servation laws," J. Comput. Physics, 49, pp357-393
(1983).
[7] J. L. Steger and R. F. Worming, "Flux vector splitting
j
of the inviscid gasdynamic equations with application
to finite-difference methods," J. of Comput. Phys, 40,
PP263-293 (1981).
[8] H. L. Stone, "Iterative solution of implicit
approximations of multidimensional partial
differencial equations," SIAM J. Numerical Analysis,
5, No. 3, Sept. (1968).
[9] A. J. Chorin, "A numerical method for solving
incompressible viscous flow problems," J. Comput.
Phys, 2, ppl2-26 (1967).
[10] S. V. Patankar and D. B. Spalding, "A
calculation procedure for heat, mass and momentum
transfer in 3-dimensional parabolic flows," Int. J.
Heat Mass Transfer, 15, ppl787 (1972).
[11] R. H. Ni, "A multiple grid scheme for solving the
Euler equations," AIAA J., 20, ppl565-1571 (1982).


-91
[12] A. Jameson, "Solution of Euler equations for two
dimensional transonic flow by a multigrid method,"
Appl. Math. Comput., 13>y pp327-355 (1983).
[13] D. C. Jespersen, "Design and implementation of a
multigrid code for the Euler equations," Appl. Math.
Comput., 13, pp357-374 (1983).
[14] N. L. Sankar, "A multigrid strongly implicit procedure
for 2-D transonic potential flow problems,"
AIAA-82-0931 (1982).
[15] P. W. Hemker and S. P. Spekreijse, "Multiple grid and
Oshers scheme for the efficient solution of the
steady Euler equations," Appl. Numer. Math., 2,
PP475-493 (1986).
[16] U. Ghia, K. N. Ghia and C. T. Shin, "High-Re solution
for incompressible flow using the N-S equations and a
multigrid method," J. of Comput. Fhys., 48, pp387-411
(1982).
[17] D. J. Mavriplis, A. Jameson, and L. Martinelli,
"Multigrid solution of the Navier-Stokes Equations on
triangular meshes," AIAA-89-0120 (1989).
.[18]J. Linden, G. Lonsdale, B. Steckel, and K. Stuben,
"Multigrid for the steady-state incompressible
Navier-Stokes equations," 11th International
Conference on Numerical Methods in Fluid Dynamics,
Williamsburg, Virginia, June 27-July 1 (1988).
[193J. H. Morrison and M. Napolitano, "Efficient solutions
of two-dimensional incompressible steady viscous
flows," ICASE report No 86-68 (1968).
[20] S. McCormick and J. Thomas, "The fast adaptive
composite grid method (FAC) for elliptic equations,
Math. Comp., 46, 174, pp439-456 (1986).
[21] S. McCormick, Multilevel Adaptive Methods for Partial
Differential Equations, SIAM, Phil. (1989).
[22] C. Liu and S. McCormick, "Multigrid,-elliptic grid