A COMPARISON OF MULTILEVEL ADAPTIVE METHODS
FOR INCOMPRESSIBLE FLOW IN GROOVED CHANNELS
by
Michael Patrick O'Brien
B.S., Boston University, 1988
A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado at Denver
in partial fulfillment
of the requirements for the degree of
Master of Science
Applied Mathematics
1992
iP r
This thesis for the Master of Science
degree by
.Michael Patrick O'Brien
has been approved for the
Department of
Mathematics
t W
Date
O'Brien, Michael Patrick (M.S., Applied Mathematics)
A Comparison of Multilevel Adaptive Methods for Incompressible Flow in Grooved
Channels
!
Thesis directed by Professor Stephen F. McCormick
ABSTRACT
Incompressible flow in a periodically grooved channel is investigated by
direct numerical simulation, using the finite volume method on a staggered grid. A
secondorder, fullyimplicit time marching scheme is used together with the
multigrid full approximation scheme (FAS). A comparison of convergence rates has
i
been made between the Fast Adaptive Composite Grid Method (FAC) and the
PseudoAsynchronous Fast Adaptive Composite Grid Method (PAFAC). A local fine
grid is placed about the grooved cavity to achieve Increased accuracy, without the
[
I
need for a global fine grid. PAFAC. and therefore AFAC, have been shown to be a
i
viable solution algorithm for the general geometry fluid flow problem.
This abstract accurately represents the content of the candidate's thesis. I
i
recommend its publication.
j Signed
i
Stephen F. McCormick
m
ACKNOWLEDGMENTS
i
I would like to thank Dr. Stephen McCormick for his guidance, support and
humor. I am also grateful to Dr. Chaoqun Liu for his guidance and helpful
.1
discussions.! Without the patience of both of these men many things would have
i
gone unlearned. A thank you is also In order to Dr. Thomas Manteuffel for sitting In
! 1
on my review committee. Appreciation Is given to Martin Marietta for providing the
. 'I i
funding necessary to complete the degree. I would also like to thank everyone at
home who put up with me during the time I was engaged in this project.
I
I
I
I
IV
CONTENTS
1.
2.
2.1
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
References
Introduction...........................................1
Geometry and Computational Domain......................3
Staggered Grids........................................5
Initial Flow Conditions................................6
Governing Equations and Boundary Conditions............8
Discretization of Governing Equations..................10
Intergrid Transfers....................................18
Multilevel Adaptive Methods............................22
Fast Adaptive Composite Grid Method:...................26
PseudoParallel Asynchronous Fast Adaptive
Composite Grid Method..................................29
Distributive Relaxation................................32
Computational Results..................................38
Conclusions............................................45
.......................................................46
v
1. Introduction
I r
The determination of fluid properties in dynamic flow situations is a topic
that has received Increasing attention from both the private and the academic
i :
sectors. Arguably the most popular form of solution for a fluid dynamics problem Is
some form of numerical flow simulation.
To accurately simulate a flow using a numerical method, the full nature of
the flow must be captured In the mathematical model. Regions of the flow
containing large gradients In the variables of Interest should naturally receive more
attention than the regions containing relatively benign flow. One technique that
addresses this purpose is a redistribution of the grid points Into the regions of
Interest, keeping the total number of points constant. While this method works for
many cases, 'it ofi:en results In grid elements with Inappropriate aspect ratios
(element length /element width), leading to Inadequate resolution. Another
numerical technique, which will be considered In this thesis, that fully captures the
regions of interest Is the multilevel adaptive method. In this scheme, multiple
coarse grids are combined with an Imbedded fine grid patch to Increase the
resolution to;a desired level in the regions of Interest, with a relatively small Increase
l 1
In the total number of points. Further, utilizing a uniform grid spacing for both
coarse grid and patch grid levels simplifies the construction of an algorithm which
Incorporates multilevel solvers. Taking into account the simple uniform grid
1
structure Implemented, a composite grid may be constructed using points from both
the global grid and the patch. The discrete flow equation may then be posed on the
composite grid arid equations on the constituent uniform grids may then be used to
accelerate the solution process.
One possible Iterative multilevel scheme solves each level of a nested
sequence of grids!in series during a socalled V cycle. While this method works well
and yields excellent results, there Is still room for Improvement given the options
and power available In today's computing environments. One such Improvement Is
based upon a simultaneous processing of all grid levels, yielding a parallel version of
the multilevel method. A pseudoparallellzed version Is used here. The code does
I
not actually solve the grids In parallel, but rather simulates parallel computation by
reserving and using previous Iteration values to allow Independent solution of the
current grid levels.
This thesis utilizes standard and pseudoparallel versions of a multilevel
i
adaptive method to obtain comparative solutions for flow transition In a periodically
grooved channel. The residuals of the solutions are then used to compare the
i
convergence factors of the two methods, to verify the efficiency of the pseudo
parallel multilevel method.
2
2. Geometry and Computational Domain
j
The model used to illustrate the multilevel methods presented is flow in a
channel thatis periodically grooved in the streamwlse flow direction (Figure 2.1).
I
An Infinite Periodically Grooved Channel
Here, after nondlmenslonallzlng with respect to the channel halfwidth h, a is the

groove depth; Cis the groove length, and L is the separation distance between
j
grooves. L will also be used as the wavelength of the specified mode of disturbance
added to the flow The channel is assumed to be infinite in extent in the streamwlse
(x) and spanwise (z) directions. This corresponds to a situation in which the flow Is
i
i
fully developed, implying an infinite number of grooves.
I
3
1
For the tests conducted here, one period of the Infinite length grooved
channel has been used as the computational domain (Figure 2.2). Since the regions
of interest In this flow field lie along the south wall and In the groove, a fine grid
patch has been Included (Figure 2.3) to obtain a greater resolution In these areas.
The overall grid Is structured to contain 5 levels of coarse grids coupled with a fine
grid patch and a global grid. All of the grid levels use a uniform structure to
facilitate discretization and multilevel Implementation.
1
* 2(h)
V
^ a
l
Figure 2.2
Computational Domain Within the Grooved Channel
Figure 2.3
Channel Groove With Fine Grid Patch
2.1 Staggered Grids
A staggered grid scheme has been used on all grid levels. The staggered grid
cell (Figure 2.4) Is defined such that the flow parameters sire calculated at the
I
midpoints of1 the sides of the volume boundaries, rather than at the grid points, with
the pressure being calculated at the center of the volume. The use of this type of
element geometry alleviates the need for pressure boundary conditions at the walls
of the channel.
Figure 2.4
Staggered Grid Element
5
3. Initial Flow Conditions
To begin the simulation, an assumption must be made about the state of the
flow. Here, the Initial flow conditions assumed to be present are defined as uo, vo,
and Po. the solution to the steady state Poiseuille flow equations:
i
u0(x,y,t) = u(x,y,0) = 1 y2
v0(x,y,t) = v(x,y,0) = 0
2
Po(x,y,t) = P(x,y,0) = x + const
! Re , (l)
which represent the flow situation depicted in (Figure 3.1). I
Figure 3.1
Planar Channel Flow
I
6
I
The Initial boundary conditions may then be assumed as If on the walls of a smooth
channel:
u(x,^l,t) = v(x,l,t) = 0
u(x,lU) = v(x,l,t) = 0 (2)
i i
j
The evolution of the flow in the grooved channel then progresses from the
:
Initial Polseuflle flow conditions. The flow description previously discussed defines a
I '
steady state solution, not a terribly Interesting problem to study. However, If a
; i
disturbance is Introduced Into the stationary Polseuflle flow, the flow field will then
I '
generally become time dependent. The perturbation of the primitive variables that
. ij
will be used here to develop the flow Is the model of surface roughness based upon
periodic grooves.,
I
i
i
j
i
7
4. Governing Equations and Boundary Conditions
The solution of the flow field In a periodically grooved channel (Figure 2.1) Is
l
of Interest since it is one possible model of the effect of surface roughness on flow
transition.
The twodimensional, time dependent Incompressible NavlerStokes
equations have been utilized as the mathematical flow field model:
dt dx dy Re ^dx2 ' Ac"0
dv <9uv " dw i 1 r d2v > a j 1 *0 (A\
dt dx dy Re
da dv ii n (5)
dx dy u,
where u and v are the velocity components In the (x) and (y) component directions, P
Is the pressure and Re Is the Reynolds number based on the mean centerline
velocity, u, the channel halfwidth h, and the dynamic viscosity of the fluid, V.
Re
uh
v
(6)
8
I
I
I I
i
i
The boundary conditions augmenting equations (3), (4), and (5) are
i i'
i
u = 0 and v = 0 (on solid walls)
u(x + L,y,t) = u(x,y,t)
v(xfL,y,t) = v(x,y,t)
P(x + L,y,t) = P(x, y, t) c(t)L t
where c(t) Is defined by the volume flowrate relation:
Q= ifu(0,y,t)dy = ^
i 3
i
9
i
5. Discretization of Governing Equations
An Implicit secondorder finite difference scheme is used to develop the
[
discrete equations on the interior of the globed gild and the patch. Since both
i
regions are of a uniform grid structure, they will be handled In a similar manner.
The finite volume method (Figure 5.1) given by Liu, Liu, and McCormick in [1] has
i
been used to. obtain the discretization of the governing equations (3). (4), and (5).
1 >u11,j+10 4Vi1J+1 \ >V.o 4Vi.M
p u H1,1 > 11,J c 
uii,jic $ ,Vmj t i >u t
U,
1+1.J
Figure 5.1
Finite Volume Method for xMomentum Equation
The x direction momentum equation (3) may be rewritten in conservative form as
du ; d
"a+ *
i
f
1 du
/
uu+ P +
Re dx ) dy
uv
1
Re dy
\
= 0.
10
I
Integrating the x momentum equation over the control volume of (Figure 5.1) and
applying the jdtvergence theorem yields the discretized form:
AxAy +
"+1 "+1
uw, w
V
f
Re
n+1 ^
n+1
n+1 n+1
us 1 T s'
Re
da
\fyj
Re
(
n+1, A 'in i ,,nl j /^uY+1 ^
iy+
i
3un+1_4un +un
JU1,J ' ; i.J U1.J
; At
ydxj
+ Rn+1
e y
Ay+
y
un+lvn+l
n n
v
Re
Ay
'daV^
Ax
i
Ax = 0
s y
(7)
The terms at the volume boundaries may be found by averaging or central
differencing as
u"+1 4 Ki+<)
1+1 Ww (ci+O
un+l IS' (C. + "u)
u+1 % (>C+C)
n+i n fvn+1 +vn+1 (: i,j+i v 1i.j+i
Vn+1 T s 1
i
j
11
$*!& 41 a1 a a*i a*
)n+l n+l nn+1
_ Ui+l.J UU
Ax
vn+1 n+l ______..n+l
' II; I U_ ;
y
Ax
Y+1 nn+1 _nn+1
_ ui.j+l ul,j
y
Ay
\n+1 n+l n+l
_ Ui.j Ui.Jl
. Ay
In a more general form (Figure 5.2) the xmomentum equation may then be written
as
^SUS ^.pUp ^P^P ^uf
12
where
, 1 ! 1 ^Ug+1 + Up+1 N
Re Ax2 2Ax 2 J
1 1 'u^ + u"+1'
Re Ax2 2 Ax 2 v
; i 1 I yll+1 ^ vil,j+l ^vi,j+l
Re Ay2 2Ay 2 V >
1 1 'vn+1 + vn+1 ^ vil,j T vi,j
Re Ay2 2Ay 2
Ap 2^t +(^E "*^w An +As)
C =C =
^ w Ax
s _ug~14ug
0 o  UN o o
 o 0 1 Â£ o Pp lmc 0 1 _
o o  i? o o
Figure 5.2
Generalized Grid Neighbor Points of Up
13
Generalizing the neighboring grid points of vp in a like fashion (Figure 5.3)
the ymomeritum equation may be rewritten as
l .1
I
BevJ + Bwvw + Bnvn + Bsvs Bpvp DpPp + DsPs = Sv. (9)
i
t i f i t t
1 O t O f PNo tVN o t
1 1 1 o i t 1 o fVw 1 PpO tv' 1 o fE
1 1 o ; ^ 1 o f 1 ps 1 o
Figure 5.3 Generalized Grid Neighbor Points of vp
For south wall points (Figure 5.4) equation (8) must be modified by substituting in
I
the expressions ,
I
ur=t
n i.1
V i
n+1 = 0
i
i
'du\
dy J.
3uf; un+
i.j
i.j+1
Ay
(10)
I i!
which are based on polynomial interpolation of u, Jt u, J+1, and u,.
14
U;
1+1.J
South Wall Point Neighbors
The general form (10) for the south wall point equation is derived through a
substitution of the terms
1 1 1 1 :i 1 ' Ug+1 + Up+1 N
Re Ax 2 1 2 Ax 2 J
_ 1  1
Re A*2 2Ax 2 ,
i 1 f vn+1 4v"+1 ^ v i1, j+1 vi,j+l
Re Ay2 2Ay 2 j
As = 0i
AP 2 At"^^E + ^w + ^N + ^s^ +
i
1
1 2
Re Ay1
C =C = 
p Ax
_ Up"1 4Up
; 2At
(ll)
15
Following; the same method, the required continuity equation may be
I
discretized as
Fe^E' ^pup"^^nvn ^pvp
(12)
where
FE=FP=
E ;p 'Ax
gn = Gs =
: Ay
S =0.
Due to the nonuniform geometry around the cavity, a special method, given In [1] Is
I
required for the cavity lip points (Figure 5.5), resulting In a discretization of the form
3umtl 44u*;. +U*'1 f
11^AxAy+
At
ur.nr_J.rÂ§f+P;
RevdxJe
n+I N
n+1
U
[ un+1u+1f1 + Pn+1
l w ReUJ. ^
^ (
Ay +
u"+1v:+1 
D D
V
Re
Ay
'jkY"
&). ,
Ax
1 ["I^uVAx 1 (
Re
\By) 2 Re v ),2 2
= 0,
(13)
I
I
16
St S' 4'IS'
where
' vn+l 3"+l_______
_ >.j 3 .j+l
M ~ ^
sin+l
Jil
un+1unt\
_ l.J l.J1
Ay
V 
V
\
V
V
u______________
Figure 5.5
Cavity Lip Point Geometry
To obtain the general form In (10) for the cavity lip points, we substitute
AT = A?!d 4
1 1
Re 6Ay2
A
new
P
+
_1___1_
Re Ay2
(14)
17
I
I
6. Intergrid Transfers
The desire to use nested gilds and fine grid patches to Implement the
multilevel adaptive methods complicates the solution process due to the handling of
the grid Interface points. Combining the Interface adjustments Into this algorithm
requires steps to be taken to ensure that the Interfaces between the patch and
global grid are handled correctly. The construction of the grid scheme and
Interfaces between the gilds Is essential to the success of the method.
The algorithm used to solve the flow field equations employs the lntergrld
transfers defined In [2] everywhere except at the cavity lip points. Relaxation
requires a restriction of the approximations and residuals from fine grid levels to
coarse gild levels. For restriction of the approximations, the following stencil
! ,i
operators should be used:
C(v):
If (P>:
i
_i
2
i
.2.
~j_ r
.2 2.
i I
4 4
J_ l
A 4.
(15)
I
18
For the restriction of residuals, a full weighting stencil Is used, determined from the
socalled area law developed by Liu [3]:
i^b (R):
I?( R):
:i 1 r
8 4 8
,1 1 1
.8 4 8.
'1 r
8 8
1 1
4 4
1 1
.8 8.
'1 r
4 4
;l 1
:4 4.
(16)
I
19
When Interpolating the corrections from a coarse grid to a fine grid, bilinear
Interpolation Is used, which Is given by the following stencils:
3
4
1_
4
(17)
The cavity lip: points are handled using separate stencils, given in [1], due to the
peculiar geometry surrounding the points (Figure 6.1):
iLCu)i
3
I
b
2b
16
0
10'
16
_3_
16.
(18)
20
O ph
t
1 2h
2h
2h
Interpolation Operator at Cavity Lip Points
21
7. Multilevel Adaptive Methods
I :
i
j i
The purpose In using a numerical Iterative scheme Is to begin with an Initial
guess and, through a relatively simple process, reduce the error In that guess to
i 'i
within some acceptable tolerance of the actual solution. A natural way to augment
i 3
such an Iterative method Is to formulate an Initial guess that Is closer to the actual
I :
solution than the'original guess was, at a low relative cost. This can be done
naturally when using a multilevel adaptive method by first approximating the
solution on the coarse grid levels and then, through successive Interpolation and
standard multilevel cycles, develop approximations on Increasingly finer levels.
Using the solution to the coarsest grid as an approximation to the next finer grid
!
and repeating this up to the finest grid level provides a good estimate to the solution
on the finest grid If the solution function is smooth, from which It may be concluded
i 
that the error, after a few standard multilevel cycles must then be comparable to the
i
discretization' error itself (which is all that may be expected of a practiced scheme).
I .
At thej heart of the algorithm, this form of embedded iteration employs a
i ii
multigrid basic cycling scheme. This scheme Is heuristically explained, by
i ;
McCormick [4], as follows:
j
1) Perform relaxations on the finest grid level.
2) Perform relaxations on progressively coarser grid levels, down to the
! j
coarsest level, to reduce smooth error components.
22
I
3) After sufficiently good approximations are achieved. Interpolate
approximations to the finer grid levels.
Many relaxation schemes have a "smoothing property" In that they selectively
eliminate high frequency modes In the solution. By selecting a relaxation scheme
that Incorporates such a smoothing property, an algorithm which employs the
relaxation coupled with the multlgrld basic cycling scheme then addresses both the
i
smooth and oscillatory components of the solution.
While! this method works well, so far it addresses only linear models. Most
problems require an algorithm capable of addressing nonlinear models. The Pull
Approximation Scheme (FAS) Is such an algorithm. FAS uses a nonlinear
adaptation of the imultlgrid basic cycling scheme. A nonlinear relaxation Is used,
which maintains a good approximation to the equations. Moreover, the coarse grid
corrections formulated by FAS are also nonlinear, further enhancing the method.
When' attempting a solution of the original equation with an uncertain Initial
guess, whether the model Is linear or nonlinear, embedded Iteration Is a useful
approach. In! the worst case, the real solution would be purely oscillatory In nature
and no Improvement In the trial solution would occur. But whether or not an
Improvement In the accuracy of the trial solution occurs after application of the
I
embedded iteration, the error eh = uh uh should then be oscillatory. If any
i j(
smooth modes had existed In the solution, they would have been well approximated
on the coarser grids. A properly designed relaxation scheme should be able to
resolve the oscillatory modes of the solution efficiently, so applying such a
23
relaxation scheme Is the next logical step. Using relaxation here may yield quick
Initial advances, but the problem of slow asymptotic convergence will manifest Itself,
often before the convergence criterion is met. But slow convergence Indicates that
the error Is smooth. Implying that the time is right to reapply the embedded
iteration scheme. that Is, multilevel coarse grid correction.
An Important point in the development Is that, at this point, the embedded
iteration scheme should be used on the residual error ? = f Au to estimate the
correction to the solution, rather than continuing to estimate the solution Itself. By
estimating the error. It Is then possible to transfer the results between grids and
update the current approximation to Uh. Additional relaxation sweeps may then be
performed to reduce the error further. Repetition of this standard multilevel cycle
suggests that applying embedded Iteration to the residual equation to reduce the
smooth error components, and then applying the relaxation scheme on the original
equations to reduce the oscillatory error components, would be an excellent solution
technique. This Is of course the basis of the multilevel methods.
Often slow convergence rates plague each level of the grid refinement, up to
and Including the fine grid. Coarsening the grid repeatedly until a level is reached
with acceptable relaxation convergence factors leads to the next level of complexity
, i
In the multilevel method. Each coarse grid would then provide a good
approximation to the next finer grid, reducing the solution effort on the fine grids.
24
The multilevel scheme utilized for this study Is based upon an Iteration of the V
I
Cycle, where a single VCycle Is given by
Uh
I
I
and Is defined In [4] as:
Step 1: Relax v, times on Ahub = fh starting with initial guess uh.
Step 2: If h = h,, go to Step 4. Otherwise, set rh = fh Ahuh,
f2h = I2hr\u2h <0,u2h <MV2h(u2\f2h).
Step 3: Set uh < uh +Ibhu2h.
Step 4: Relax v2 times on Ahuh = fh starting with initial guess uh. (2q)
This suggests an iterative scheme In which successively coarser grid levels are
i
employed until an acceptable convergence rate for relaxation is found. Then, from
the coarsest grid level, the algorithm would proceed back through the fine grid
levels, making corrections to the solution at each level.
I
25
8. Fast Adaptive Composite Grid Method
I
The above, solution recipe will be augmented by the addition of a composite
grid to further enhance the procedure. The formation of a composite grid that
: i
; 'I
utilizes the global coarse grid, the local line grid patch, and the Interface points
. ' 11
provides a more accurate discretization to guide the solution process. The definition
and Incorporation of the composite grid into the calculations will be accomplished
I
via the Fast Adaptive Composite Grid (FAC) method. FAC is an Iterative mesh
i
refinement method applied to the composite grid residual equations, which are used
to guide Improvements to the classical process. This method, adapted from [5], may
be described as follows:
I ;
1) Solve the problem beginning with the coarsest grid down to the level of
the global grid, obtaining the boundary values and initial guesses for
i ,
the patches by Interpolation from the global grid solution
Lhuh=/h
(21)
2) Perform distributive relaxation on the patch
Lhuh = A
.22 2
(22)
26
I
3) Form the composite grid defect, which consists of residuals from three
regions: the coarsegrid Interior points, the patchgrid Interior points,
and the coarsepatchgrld Interface points:
I
I
R
, v composite
f Lu on global grid
/! L1 u1 on Interface
// ~ LhuI In patch.
(23)
4) Restrict the composite grid residuals to the global grid.
;i I
5) Solve the defect equation on global and patch grid levels:
LX =1
LJX +IhRcompo8ite
2 2
A
In patch
elsewhere.
(24)
6) Interpolate the globalgrid solution to define patch values:
for patch Interior points
' 2 2
^h^Uh for patch boundary points.
1 2 (25)
The FAC procedure has been coupled with a recursive multilevel VCycle method,
i I
resulting In a jveiy efficient overall process.
i
i
l
27
I I
Given the current availability and relatively low cost of multiprocessor
computers, at natural extension to FAC would allow for the parallel computation of
the various gild levels. The parallel version of FAC has been termed the
I
Asynchronous Fast Adaptive Composite Grid Method (AFAC). AFAC performs each
level of refinement simultaneously. The AFAC method inserts a new grid level
between each pair of adjacent grid levels that has the resolution of the coarser grid
and the domain of the finer grid. These new grids are then used to eliminate
correction components duplicated by the various grid levels, avoiding the need for a
i
regressive underrelaxation scheme.
The FAC method and an adaptation of the AFAC method have been used in
this thesis to compare results obtained from the periodically grooved channel flow
problem.
28
I
I
i
9. PseudoParallel Asynchronous Fast Adaptive
Composite Grid Method
At this point in its genesis, the development of the FAC multilevel procedures
1 i
described so far could be used to solve many flow problems effectively. However,
given today's capabilities and the expectations of future computing environments, a
version of this code that would solve all grid levels in parallel, with no appreciable
! i
loss of accuracy or convergence rate, could dramatically reduce throughput time for
I j
the solution and would be highly effective in large scale simulation. For Instance,
on a hypercube computer, an upperlevel host processor to an AFAC code could be
used to optimize the solution time required. The processor would be used to
evaluate the number of points on a grid, thereby obtaining the work required to
solve each level. Then, based on the CPU availability of the hypercubes constituent
i 
processors, the host processor could farm out the grid levels to the machines having
the necessary amount of available computational power. This type of method would
i ;
help to insure that the most is made of the available computing resources, while
finishing the simulation in a minimum of time. While this seems compelling, is
there a way to devise such a multilevel scheme and Implement it in code so that it
will run efficiently in this maimer and continue to exhibit fast convergence?
A useful step to be taken before answering this question, which requires
substantial effort to build a truly parallel version of the code, is to create artificially
29
a parallel version of the FAC code which would run In single processor mode, a
pseudoparallel AFAC version (PAFAC). This code could then be used to determine
changes In the convergence factors between the FAC and the PAFAC
Implementations. Any deficiencies In the thought process or theoretical
development would then be forced out with a lesser undertaking than that of
creating a new code.
A drawback to creating such a pseudoparallel version stems from being
forced to base PAFAC upon a scheme which uses previous Iteration values for the
i
various grids In order to allow for simultaneous processing. The first Vcycle
Iteration through the algorithm is made using the initial trial solution. All
successive Vcycle passes are made using the values obtained from the previous
cycle. This procedure does not affect the implicit nature of the algorithm, since as
updated values produced on the current cycle become available they are Included in
the calculations. The values in current use are obtained simply from the previous
iteration, while the newly generated values are held In storage for application In the
next Vcycle.
During  the conversion to the PAFAC code, many elements of the FAC
algorithm were updated to reflect the holdover of values. After the Initialization and
residual calculation, the values used in the Iterative VCycle scheme are delayed by
one cycle and utilized to begin the next Iteration. Following the relaxation phase,
the delayed residuals are then passed to the restriction routines, where the
corrections are! made from the finer grid to the global grid, via the composite grid.
30
It is Important to note that, by using the previous iteration values to generate
the PAFAC algorithm, a degradation of the true convergence rate occurs. By forcing
the delayed values to be used for calculation In the current Iteration, some of the
! ;
effectiveness Is lost since the corrections being made In the current Iteration were
generated for the last iteration. This algorithmic Introduction of a sort of retardation
Into the solution process serves to slow the overall convergence. However, this
realty Is the basis for our ability to perform AFAC. It enables the simultaneous
processing Of all refinement levels. The loss of speed is not significant as the
following results will show.
i
31
10. Distributive Relaxation
Throughout the previous development, many references have been made to
"relaxation techniques." Many different types of relaxation techniques have been
developed with varying degrees of success, depending on the application. During
the process of formulating the discretized model equations, care must be exercised
to ensure that the relaxation scheme chosen will be an effective smoother for the
discrete equations.
Once the Incompressible NavlerStokes equations have been put Into a
discretized matrix form, it Is often noted that many conventional iterative solvers fall
to converge due to a lack of diagonal dominance In the resulting matrix. Several
attempts have been made to Improve upon this system of equations to allow an
Iterative solution. Chorin [6] and Patankar and Spalding [7] detail two such
i
examples. Unfortunately, such methods often yield Inconsistent physics or very
slow convergence. A multigrid approach based on distributive relaxation that does
not suffer from these troubles was Introduced by Brandt [8] and Liu [31. This
method In effect uses a preconditioner to Improve the characteristics of the
discretized system. Implementation of distributive relaxation may be accomplished
as follows.
I
32
First, define the new variables wl, w2, and w3 as
u j 1 o" 1 o ^4 "wl"
V I 0; 1 5y w2
p 1 o o 0 3 1 w3_
where
Q(h=5;+u5]l+Sv<5y^Ah,
S~ is the backward Euler operator,
6, and <5 are central difference operators,
and
Ah is a discrete Laplacian.
(26)
Next, define the discretized matrix form of the NavlerStokes equations as
O B O 1 u fl'
0 Qh, V f2
1 J* & o 1 p _f3
j
I
Applying (26) to (27) yields
i
Q* 0 S, u Q o S, 1 0 1 1 'wl'
o Q, s, V = 0 Q, S, 0 1 8y w2
1 p i i o 0 3 1 w3_
,1
I
(27)
(28)
I
33
I
Equation (27) is then related to the linear system
Q 0 5, ' "wl"
o :q* s, w2 = f2
l 1 & _w3. _f3_
i
(29)
Even though the system In (29) only applies In special cases, It is easily
solved by conventional methods, and may be used as a basis to develop an
applicable method. Adapting this method to the flow problem results In the
I
following scheme.
Begin; with the system
P P
AEuE+Awuw+ANuN+AsusApup+*L = St
Ax
Beve + Bwvw + Bnvn + Bsvs BpVp + 5* PjL = S..
Ay
uE uP vn ~ vp = s
Ax
Ay
(30)
I
34
Then perform the distributive relaxation according to the following recipe:
I
1) Freeze P
Perform point GaussSeldel relaxation on the momentum equations (30)
l
to obtain a new u and v
2) For each grid cell (Figure 10.1) modify the velocities to
l
satisfy,the continuity equation
Next, modify the pressures so that the new velocities satisfy
the momentum equations
The velocity updates are calculated by determining 6 and P such that
(ue+/35)(up/J<5) (vN + <5)(vp 5) 0
Ax + Ay S"
8 =
S 
'"ueUp Vn~Vp^
Ax
Ay
(
2Â£ + A
Ax Ay
(31)
and
35
I
Figure 10.1
Distributive Relaxation
Now that S arid fi have been determined, let
UE
up
VN <j VN +5
vp vp 6. .
(32)
The pressure P and the neighboring pressures are updated according to
+ 5Pp
Ue ' )b8 + AB8+SPw 5Pp
2Ax, Ax Ax
V. ' +Bp<5+ **. = 0.
2 Ay, 1 Ay Ay
(33)
1 36
I '
At this point, the velocities are the old values and 8 Is the correction from (31).
I
Choosing
3
2At
i r 2 2
H + 
Re^Ax Ay
P&x8
(34)
allows an easy determination of <5PW and 5PS from (33):
5PW 5PP Ap/?AxÂ£ 
1
Re Ax 2 Axj
Ax8
5PS = 5Pp BppAy8 
1
u
Re Ay 2 Ay
Ay 8.
(35)
The determination of SPE and <5Pn Is analogous.
Due to the selection of a small At, the pressure correction term <5Pp Is much
I
larger than the other pressure corrections <5PE 5PW <5Pn and <5PS for this
problem. Therefore, for this algorithm, the smaller pressure correction terms are
I
assumed to be zero, while the main pressure correction remains unaffected:
<5PF =, 5PW = 5Pn = <5PS = 0
8P =
3 1
2 At Re
Ax2 Ay
j8Ax8
(36)
I
37
I
Utilizing this form of relaxation has allowed the solution of the NavlerStokes
equations for the flow In a grooved channel to proceed without the usual
convergence problems or diversions from the flow physics usually associated with
this class of problem.
I
38
I
11. Computational Results
The differences In convergence between the PAFAC algorithm and the
i
standard FAG algorithm are demonstrated through a comparison of the convergence
factors between Iterations of the two methods. Basic multlgrld cycles, used in both
FAC and PAFAC, have been shown [91 to converge at a linear rate that Is
Independent of mesh size, with typical convergence factors of 0.2 or less. The
parallel nature of PAFAC degredates the convergence characteristics from the 0.2
value. In fact; FAC asymptotic convergence factors [9] are roughly the square of
,i
those obtained for AFAC. Since the PAFAC code used here emulates an AFAC code,
similar results are expected.
A disturbed Polseullle flow, discussed previously, will be used for
i '
Initialization. ; A variation study will be performed for a range of Reynolds numbers.
The Reynolds number parameter Re was chosen as a critical parameter due to Its
l
substantial effect upon the overall nature of fluid flow. One may Interpret the
definition of Re In (6) on a slightly more heuristic level as
Up = .Inertial Forces , (37)
Viscous Forces
39
I
As the inertial forces increase and dominate the viscous forces, the fluid flow
i
becomes more turbulent. A range of Reynolds numbers will then present the
opportunity to test fully the algorithm.
In order to compare the FAC and PAFAC codes, some basic ground rules
I
were established. The grid structure for each method was identical, with the same
initial conditions. The global mesh in each case was 10 X 18. with a patch
i
refinement of 18 X 18. Since the multilevel schemes have been based upon the
i
residual error to estimate corrections to the solution, rather than by estimating the
solution itself, the convergence factors will be calculated using the residuals
determined when exiting the fine patch grid at the end of a cycle:
I
r.+i
Convergence factor =. (38)
' ri
.1
In both codes the norm used to calculate the residuals is the (2 norm, given by:
!
I
This norm was chosen due to the ease of calculation given the grid structure
employed.
I
40
Each method has been limited to one cycle per Iteration, except for the
iteration chosen to make the convergence factor comparison. Iteration number 10.
In iteration 10. the FAC version of the code performs 5 cycles, while the PAFAC
i
version of the code performs 8 cycles. These parameters were chosen with some
i
malice of forethought, expecting the FAC version to perform in a superior fashion.
The residuals obtained from the u and v momentum equations are presented in
i .
Tables 11.1 11.4, for FAC u and v and PAFAC u and v respectively.
Reynolds number 225 1000 7500 75000
Cycle number
1 9.10E05 1.54E04 2.19E04 2.33E04
2 5.74E06 6.55E06 9.54E06 1.04E05
3 ! 1.15E06 7.47E07 1.25E06 1.55E06
4 ! 3.9 IE07 1.78E07 2.55E07 3.09E07
5 7.74E08 4.51E08 1.55E07 2.01 E07
Table 11.1
FAC u Residuals
l
Reynolds number 225 1000 7500 75000
:
Cycle numberi
1 .! 1.72E05 2.29E05 3.35E05 3.60E05
2 3.26E06 3.28E06 3.62E06 3.77E06
3 , 6.92E07 7.35E07 9.32E07 1.00E06
4 ! 1.92E07 1.75E07 2.26E07 2.43E07
5 : 4.63E08 4.4 IE08 7.89E08 9.11E08
Table 11.2
FAC v Residuals
41
i
I
Reynolds number 225 1000 7500 75000
Cycle number
1 1 1.15E04 2.31E04 3.79E04 4.10E04
2 3.03E05 1.29E04 2.60E04 2.92E04
3 : 1.50E05 4.34E05 8.57E05 9.62E05
4 8.49E06 1.82E05 2.69E05 2.86E05
5 1 1.86E06 7.62E06 1.4 IE05 1.55E05
6 ; 6.7 IE07 2.9 IE06 5.32E06 5.81E06
7 1.34E07 3.49E07 5.52E07 5.92E07
8 7.94E08 1.73E07 2.74E07 3.08E07
Table 11.3
AFAC u Residuals
Reynolds number 225 1000 7500 75000
Cycle number 
1 ; 9.22E05 5.71 E05 8.18E05 8.84E05
2 9.76E06 1.65E05 3.02E05 3.36E05
3. ! 6.59E06 8.72E06 1.31 E05 1.42E05
4 ! 4.16E06 5.77E06 7.19E06 7.38E06
5 1.04E06 2.73E06 4.08E06 4.31E06
6 ! 2.58E07 7.00E07 1.16E06 1.25E06
7 ; 8.44E08 1.42E07 2.02E07 2.14E07
8 ! 3.88E08 7.31E08 9.94E08 1.05E07
Table 11.4
AFAC v Residuals
I
42
1
I
The convergence factors calculated for the FAC and PAFAC methods are presented
In Tables 11.511.8, for FAC u and v and PAFAC u and v respectively.
Reynolds number 225 1000 7500 75000
Convergence factor cycle numbers
12 0.063 0.043 0.044 0.045
23 0.201 0.1 14 0.131 0.149
34 0.339 0.239 0.203 0.199
45 0.198 0.253 0.610 0.651
,
Average convergence factor by root method 0.171 0.131 0.163 0.171
Table 11.5
FAC u Convergence Factors
Reynolds number 225 1000 7500 75000
i
Convergence factor cycle numbers
12 0.190 0.143 0.108 0.105
23: 0.212 0.224 0.257 0.265
34; 0.278 0.238 0.242 0.243
45i 0.240 0.252 0.350 0.375
1
Average convergence factor by root method 0.228 0.210 0.220 0.224
Table 11.6
FAC v Convergence Factors
43
I
Reynolds number 225 1000 7500 75000
Convergence factor cycle numbers
12 0.264 0.557 0.687 0.71 1
23i 0.495 0.337 0.329 0.330
34 0.567 0.418 0.313 0.297
45' 0.219 0.419 0.525 0.541
56; 0.360 0.382 0.377 0.375
67 0.200 0.120 0.104 0.102
78' 0.592 0.495 0.497 0.521
Average convjergence factor by root method 0.354 0.358 0.356 0.358
Table 11.7
AFAC u Convergence Factors
Reynolds number, 225 1000 7500 75000
Convergence, factor cycle numbers
12i 0.106 0.289 0.369 0.380
23:' 0.675 0.527 0.433 0.422
34! 0.631 0.662 0.550 0.520
451 0.250 0.474 0.567 0.584
56; . 0.249 0.256 0.286 0.291
67! 0.327 0.203 0.173 0.171
781 0.460 0.516 0.493 0.491
Average convergence factor by root method 0.329 0.386 0.383 0.382
Table 11.8
AFAC v Convergence Factors
I
I
44
The approximate averages given in Tables 11.511.8 are computed using a root
where n Is the total number of cycles performed In the iteration and k=(nl).
For both FAC and PAFAC, convergence was obtained for all Reynolds
numbers used. The FAC convergence results show good agreement with the 0.2
convergence factor achieved by McCormick [9] for simpler model problems.
Moreover, the PAFAC convergence results follow the trend reported by McCormick
[9], being roughly equal to the square root of the FAC convergence factors.
average method, as
(39)
45
12. Conclusions
Using, the flow through a grooved channel as a model, the overall
i
convergence characteristics of two versions of a multilevel Iterative scheme, FAC and
PAFAC. have 'been compared. It has been found that the FAC method had
convergence factors of 0.23 or better, while the PAFAC method had convergence
factors of 0.39 or better. These convergence factors appear to support the results
j
drawn by McCormick (4], [9], for both FAC and AFAC type codes on simpler model
problems. The results presented here are especially significant. The behavior of
FAC and AFAC for model problems has been shown to extend to more complex
applications that include systems, nonlinearity, and changing character (i.e. high
Reynolds number).
i
Examining the numbers of cycles performed for the two methods. It may be
concluded that approximately 8 PAFAC cycles are equivalent to 5 FAC cycles. This
Is a marginal jncrease In the total amount of work required, and Is more than
compensated for by the full parallellzablllty of AFAC.
i
1
46
References
[1] Llu.C;. ,Llu,Z., and McCormlck,S.F. 1991. Multilevel Adaptive Methods for
i
Incompressible Flow In Grooved Channels. J. of Computational and Applied Math,
i
38:283295. j
[2] Liu,C.,Llu,Z., and McCormlck,S.F. 1991. Multigrid Methods for Flow
Transition ina Planar Channel, Comput.Phys.Comm. 65:188200.
[3] Llu.Ci 1989. Multilevel Adaptive Methods in Computational Fluid Dynamics,
i
Ph.D thesis. University of Colorado at Denver.
[4] McCormick,S.F. 1987. Multigrid Methods, SIAM, Phil.
i
[5] Llu.Ci 1991. Multilevel Adaptive Grids for Flow Transition in a Planar
'i
Channel, Science Report, Ecodynamlcs Research Associates, Inc., Denver.
[61 ChorinAJ. 1967. A Numerical Method Jar Solving Incompressible Viscous
i
Flow Problems, J.Comput.Phys. 2:1226.
[7] Patankar.S.V., and Spalding,D.B. 1972. A Calculation Procedure for Heat
Mass and Momentum Transfer in 3Dimensional Parabolic Flows, Int. J.Heat Mass
Transfer 15:17871806.
[8] Brandt A 1984. Multigrid Techniques: Guide with Application to Fluid
Dynamics, GMD Studie, GMD, St. Augustine.
[9] McCormick,S.F. 1989. Multilevel Adaptive Methods for Partial Differential
Equations, SIAM, Phil.
I
47

PAGE 1
I :A COMPARISON OF MULTILEVELAD.API"IVE METHODS FOR INCOMPRESSIBLE FLOW IN GROOVED CHANNELS by Michael Patrick O'Brien B.S., Boston University, 1988 A thesis submitted to the Faculty of the Graduate School of the University of Colorado at Denver in partial fulfillment of the requirements for the degree of Master of Science Applied Mathematics 1992
PAGE 2
This thesis for the Master of Science degree by _Michael Patrick O'Brien has been approved for the Department of Mathematics Date
PAGE 3
O'Brien, Michael .Patrick (M.S., Applied Mathematics) A Compartsoh of Multilevel Adaptive Methods for Incompressible Flow in Grooved I I Channels Thesis directed by Professor Stephen F. McCormick ,, ABSTRACT i . Incompressible flow in a periodically grooved channel is investigated by I cUrect numerical using the finite volume method on a staggered grid. A I I secondorder., fullyimplicit time marching scheme is used together with the multigrid full approximation scheme (FAS). A comparison of convergence rates has I been made between the Fast Adaptive Composite Grid Method (FAC) and the PseudoAsynchronous Fast Adaptive Composite Grid Method (PAFAC). A local fine grid is placed about the grooved cavity to achieve increased accuracy, without the I I need for a fine grid. PAFAC, and therefore AFAC, have been shown to be a I viable solution algorithm for the general geometry fluid flow problem. I I This abstract: accurately represents the content of the candidate's thesis. I I recommend its p1,1blicat1on. I ' Signed Stephen F. McCormick iii
PAGE 4
ACKNOWLEDGMENTS I would life to thank Dr. Stephen McCormick for his guidance, support and I humor. I ani also grateful to Dr. Chaoqun Liu for his guidance and helpful I I .I discussions. :Without the patience of both of these men many things would have I I gone you is also in order to Dr. Thomas Manteuffel for sitting in I I on my coriunittee. Appreciation is given to Martin Marietta for providing the I I I funding necesscuy to complete the degree. I would also like to thank everyone at home who p* up with me during the tJme I was engaged in this project. I I !. iv
PAGE 5
CONTENTS 1. Introduction ...................................................................... 1 2. Geomet:Iy and Computational Domain .............................. 3 2.1 Staggered Grids ................................................................. 5 3. Initial Flow Conditions ...................................................... 6 4. Governtng Equations and Boundary Conditions ................ 8 5. Discretization of Governtng Equations ............................... 10 6. Intergrid 1'ransfers ............................................................ 18 7. Multilevel Adaptive Methods .............................................. 22 8. Fast Adaptive Composite Grid Method: .............................. 26 9. PseudoParallel Asynchronous Fast Adaptive Cdmposite Grid Method ..................................................... 29 10. Distributive Relaxation ...................................................... 32 11. Computational Results ...................................................... 38 12. Conclusions ...................................................................... 45 References ......................................................................................... 46 v
PAGE 6
1. Introduction I The detenn1nation of fluid properties in dynamic flow situations is a topic that has received increasing attention from both the private and the academic I i : sectors. Arguably the most popular form of solution for a fluid dynamics problem is some form nmperical flow simulation. To accurately simulate a flow using a numerical method, the full nature of I the flow must be 'captured in the mathematical model. Regions of the flow contaJntng large gradients in the variables of interest should naturally receive more attention than the regions contaJntng relatively benign flow. One technique that I addresses this purpose is a redistribution of the grid points into the regions of interest, keeJ?ing the total number of points constant. While this method works for 1 'I many cases, it often results in grid elements with inappropriate aspect ratios (element length /element width), leading to inadequate resolution. Another numerical technique, which will be considered in this thesis, that fully captures the ' I regions of interest is the multilevel adaptive method. In this scheme, multiple I coarse grids combined with an imbedded fine grid patch to increase the resolution to ia desired level in the regions of interest, with a relatively small increase I I in the total n:umber of points. Further, utilizing a uniform grid spacing for both coarse grid and patch grid levels simplifies the construction of an algorithm which incorporates multtlevel solvers. Taking into account the simple uniform grid 1
PAGE 7
structure implemented, a composite grid may be constructed using points from both the global grtd and the patch. The discrete flow equation may then be posed on the composite grid aiid equations on the constituent uniform grids may then be used to accelerate the solution process. . One possible iterative multilevel scheme solves each level of a nested sequence of grtds:in series during a socalled V cycle. While this method works well and yields excellent results, there is stlll room for improvement given the options and power available in today's computing environments. One such improvement is based upon a simultaneous processing of all grid levels, yielding a parallel version of the multllevei method. A pseudoparallelized version is used here. The code does I not actually solve the grids in parallel, but rather simulates parallel computation by reserving using previous iteration values to allow independent solution of the current grid This thesis utilizes standard and pseudoparallel versions of a multilevel I adaptive metl;lod to obtain comparative solutions for flow transition in a periodically I grooved chanp.el. 'The residuals of the solutions are then used to compare the convergence factors of the two methods, to verify the efficiency of the pseudo' parallel multilevel method. 2
PAGE 8
l : I 2. Geometry and Computational Domain The model used to illustrate the multllevel methods presented Is flow in a I channel that'is periodically grooved in the streamwlse flow direction (Figure 2.1). a l L FJgure 2.1 An Infinite Periodically Grooved Channel Here, after n6ndimensionaltzing with respect to the channel halfwidth h. a Is the I I I groove depth; Cls:the groove length, and Lis the separation distance between i grooves. L also be used as the wavelength of the spectfted mode of disturbance I added to the how: The channel is assumed to be inflnite in extent in the streamwtse I (x) and spanWise (z) directions. This corresponds to a situation in which the flow is I I fully developed, inlplying an infinite number of grooves. I 3
PAGE 9
For the tests conducted here, one period of the inflnite length grooved channel has been used as the computational domain (Figure 2.2). Since the regions of interest 1n this .flow field lie along the south wall and 1n the groove, a fine grid patch has :b:icluded (Figure 2.3) to obtain a greater resolution 1n these areas. The overall is structured to contain 5 levels of coarse grids coupled with a fine grid patch and a global grtd. All of the grid levels use a uniform structure to facilitate disc'rettzation and multlleveltmplementatlon. y 2(h) X a L F.lgure 2.2 Computational Domain Within the Grooved Channel 4
PAGE 10
I I I I I I I I I I I 1 Figure 2.3 Channel Groove With Fine Grid Patch 2 .I Grids A staggered grid scheme has been used on all grid levels. The staggered grid cell (Figure 2.4) is defined such that the flow parameters are calculated at the I midpoints of:the sides of the volume boundaries. rather than at the grid points. with the pressure being calculated at the center of the volume. The use of this type of I element alleviates the need for pressure boundary conditions at the walls of the channel. Ut,j v f,j+1 0 pi. ,J u f+ l,j Figure 2.4 Staggered Grid Element 5
PAGE 11
3. Initial FlowConditions To begin the simulation, an assumption must be made about the state of the I flow. Here, t)Ie :lil.itial flow conditions assumed to be present are defined as uo, vo. and Po. the solution to the steady state Poiseuille flow equations: U0(x,y, t) = u(x,y,O) = 1y2 v0(x,y,t) = v(x,y,O) = 0 2 P0(x,y,t) = P(x,y,O) = x +const Re which represent the flow situation depicted in (Figure 3.1). y : 0 i .. 1 : ........... a.\ ...J 7 ..,.. Figure 3.1 Planar Channel Flow 6 (1) L .. X
PAGE 12
' i The initial conditions may then be assumed as If on the walls of a smooth channel: I u(x, ...L 1, t) = v(x, 1, t) = 0 I u(x, 1i, t)::::, v(x, 1, t) = 0 I (2) The etrolution of the flow in the grooved channel then progresses from the I :I initial Polseu:Uie flow conditions. The flow description previously discussed defines a I ,I I I steady state solution, not a terribly interesting problem to study. However, If a I I I disturbance; introduced into the stationary Poiseuille flow, the flow field will then I generally bec:ome. tlme dependent. The perturbation of the pr1mitlve variables that ,, will be used here ,to develop the flow Is the model of surface roughness based upon ' periodic grooves. .I I. 7
PAGE 13
4. Govern1Ili Equations and Boundary Conditions The splution of the flow field in a periodically grooved channel (Figure 2.1) is of interest since is one possible model of the effect of surface roughness on flow transition. The tlme dependent incompressible NavterStokes equations haye been utillzed as the mathematical flow field model: (3) (4) (5) where u and Y. are the velocity components in the (x) and (y) component directions, P I is the pressure and Re 1s the Reynolds number based on the mean centerline velocity, u, the channel halfwidth h, and the dynamic viscosity of the fluid, v. I l 'uh Re=. 'v. 8 (6)
PAGE 14
I The boundazy augmenting equations (3), (4), and (5) are: u = 0 an4 v = 0 (on solid walls) u(x f L, y, t) = u(x, y, t) v(x i L,y, t) = v(x,y, t) P(xtL,y,t)=P(x,y,t)c(t)L. I I where c(t) is defined by the volume flowrate relation: I ' 4 I I Q = Ju(O,y, t)dy =.,., 3 9
PAGE 15
I 5. Discretization of Governing Equations I I I An Inipliclt secondorder finite difference scheme is used to develop the I discrete equations on the interior of the global grid and the patch. Since both I I I regions are o'r a uniform grid structure, they will be handled in a similar manner. I The finite volume method (Figure 5.1) given by Liu, Liu, and McCormick in [1] has I I been used to1 the discretization of the governing equations (3), (4), and (5). Figure 5.1 ;Finite Volume Method for xMomentum Equation I The x directi<;>h momentum equation (3) may be rewritten 1n conservative form as 'I I : au au +P)+.E_(uv_!_ au)=o. at :ax Re ax (Jy Re (Jy 10
PAGE 16
. I I Integrating 1fe x equation over the control volume of (Figure 5.1) and applying the theorem yields the d1scretized form: '! (7) 'I The terms at .the volume boundaries may be found by averaging or central differencing as I I 1 un+l ::J + un:l) e : 2 1 l+l,J I,J un+t + w 2 1l,J l,j ,I '1 un+1 n i 2 ,. I,J+l IJ un+l + Un+l) 5 I 2 .! l,j1 i,j : n+l 1 (': n+l n+l ) V n 2 Y i,j+l + V 1l,j+l I Vn+l .!.(tn+l + : 2 l,j 1l,J 11
PAGE 17
(du)n+l : un+l. l,j 1l,J dx L\x W ( :\.. )n+l n+l n+l ou u i,j+l u l,j dy n t,_y ( :l.. )'.n+l = U U ou l,J l,j1 dy ;s t,_y In a more general form (Figure 5.2) the xmomentum equation may then be written as 12
PAGE 18
where UN 0 0 0 0 Uw u 0 0 0 E 0 Pw Pp 0 0 0 Figure 5.2 Generalized Grid Neighbor Points of up 13
PAGE 19
Geneializ(ng the neighboring grid points ofvp in a like fashion (Figure 5.3) I \J the ymomentum' equation may be reWritten as I 0 0 0 0 0 0 Figure 5.3 0 0 0 Generalized Grid Neighbor Points ofvp I For south wap points (Figure 5.4) equation (8) must be modified by substituting in I the expressiops un+l = 0 8 I ''I vn+l::::: 0 'I 8 I I (au)' ., ay. I 1 3ua+l un+l I,J 3 i,j+l fly which are polynomial :Interpolation ofu1 ,1 u1 ,1 1 and u,. 14 (10)
PAGE 20
I I '. U. 1 1 1,J+ u=O FJgure 5.4 South Wall Point Neighbors The general.fonn. (10) for the south wall point equation is derived through a substitution :bf the terms 15 (11)
PAGE 21
' I, I Follo'o/ing; the same method, the required continuity equation may be I discretized I (12) where ' 1 F =F =B ;p 1& G =G N S L\ ,. y s =0. m Due to the geometiy around the cavity, a special method, given in (1] is .I required for the C1avity lip points (Figure 5.5}, resulting in a discretization of the form 3uuJ4u,;+u ( 1 (du)u+l ) i.i .i ;;'i i,J llxliy + Uu+luu+l __ + pa+l liy .!\t e e Re iJx .J 1 0 ( I 1 (du)D+l ) ( 1 (du)D+lJ U u+l.luu+I .. _ + p,u+l. .!\y + Uu+lvu+l _ L\x w I w R :l... l.J D D R ;:)., 1 I e oA "' e v1 ' D I _...!...([au. llx __ 1 ( au)ll+ lix = o, Re dy il 2 Re dy 2 2 i I .I I' I I 16 (13)
PAGE 22
where Figure 5.5 Cavity Lip Point Geometry To obtain the general form In (10) for the cavity lip points, we substitute 1 1 N N Re 6.dyz A new .i_ A old + _1 ___ 1_. P P Re (14) 17
PAGE 23
6. Intergrtd Transfers I I The desire to use nested grids and fine grid patches to fmplement the multllevel methods complicates the solution process due to the handling of the grid interface .points. Combining the interface adjustments into this algorithm requires step$ to be taken to ensure that the interfaces between the patch and global grid are handled correctly. The construction of the grid scheme and I interfaces between the grids is essential to the success of the method. The algorithm used to solve the flow field equations employs the intergrid I transfers def!f1ed !0 [2] everywhere except at the cavity lip points. Relaxation requires a of the approximations and residuals from fine grid levels to coarse grid levels. For restriction of the approx:fmattons, the following stencU I operators used: 1 2 l' 2 I 1.. 1 4 4 1 1 4. 4 (15) 18
PAGE 24
For the restriction of residuals, a full weighting stencil is used, determined from the socalled area law developed by Liu [3]: [I 1 il ciu): 'g 4 ,1 1 I : 8 4 1 1 8 8 1 1 4 4 1 1 8 8 [I il .:4 b I 1 ,_ : (16) 19
PAGE 25
' I When the corrections from a coarse grid to a fine grid, bllinear interpolation is used, which is gtven by the following stenclls: I [i] Ih ( ); [3'; !] 2h v. 4 4 Ih (P) [ :6.. 2h :' .1 _!_ 16 16 (17) The cavity lip: pofi?.ts are handled using separate stencils, gtven in [1], due to the I peculiar geometry surrounding the points (Figure 6.1): (18) 20
PAGE 26
0 ph t v2h 8 p2h Figure 6.1 Interpolation Operator at Cavity Lip Points 21
PAGE 27
7. Multilevel .!Adaptive Methods '. : I I I The purpose in using a numerical iterative scheme is to begin with an initial I .I I I guess and, through a relatively simple process, reduce the error in that guess to i .: I within some acceptable tolerance of the actual solution. A natural way to augment I I :1 such an iter8:tive 'method is to formulate an initial guess that is closer to the actual I : solution thari the1or1ginal guess was, at a low relative cost. This can be done :' naturally when using a multilevel adaptive method by first approximating the I solution on coarse grid levels and then, through successive interpolation and I ' I standard multllev'el cycles, develop approxfmations on increasingly finer levels. I ' Using the solution to the coarsest grid as an approximation to the next finer grid I I and repeating up to the finest grid level provides a good estlmate to the solution ;: on the finest grid tl" the solution function is smooth, from which it may be concluded I ' I I that the erroli a few standard multilevel cycles must then be comparable to the i ., discretization! error itself (which is all that may be expected of a practical scheme). I At the! heart of the algorithm, this form of embedded iteration employs a II multlgrid scheme. This scheme is heuristically explained, by McCormick (4], as follows: I .: 1) Petform relaxations on the finest grid level. 'i 2) Pefrorrli relaxations on progressively coarser grid levels, down to the I .I I I coarsest level, to reduce smooth error components. I I I 22
PAGE 28
3) After sufficiently good approximations are achieved, interpolate appraxtmations to the finer grid levels. Many relaxafion schemes have a "smoothing property" in that they selectively ellminate high frequency modes in the solution By selecting a relaxation scheme that incorporates, such a smoothing property, an algorithm which employs the relaxation coupled with the multlgrid basic cycling scheme then addresses both the I I smooth and escillatory components of the solution. Whilei this method works well, so far it addresses only linear models. Most , problems require an algorithm capable of addressing nonlinear models. The Full Approximation Scheme (FAS) is such an algorithm. FAS uses a nonlinear adaptation of the 'imultlgrid basic cycling scheme. A nonlinear relaxation is used, which maintains a good approximation to the equations. Moreover, the coarse grid corrections fqrmulated by FAS are also nonlinear, further enhancing the method. I When: attempting a solution of the original equation with an uncertain initial guess, whether the model is linear or nonlinear, embedded iteration is a useful approach. rn the worst case, the real solution would be purely oscillatory in nature and no improvement in the trial solution would occur. But whether or not an improvement in the accuracy of the trial solution occurs after application of the embedded iteration, the error eh = uhuh should then be oscillatory. If any smooth modes had existed in the solution, they would have been well approximated on the coarser grids. A properly designed relaxation scheme should be able to resolve the modes of the solution efficiently, so applying such a 23
PAGE 29
i relaxation scheme is the next logical step. Using relaxation here may yield quick I initial advances, but the problem of slow asymptotic convergence will manifest itself, I often before the convergence criterion is met. But slow convergence indicates that the error is implying that the time is right to reapply the embedded iteration scheme, that is, multllevel coarse grid correction. An important point in the development is that, at this point, the embedded iteration scheme should be used on the residual error r = fAii to estimate the correction to;the rather than continuing to estimate the solution itself. By estimating the error, it is then possible to transfer the results between grids and h update the current approximation to U Additional relaxation sweeps may then be performed to the error further. Repetition of this standard multilevel cycle I suggests that apJ)lying embedded iteration to the residual equation to reduce the smooth errorcomponents, and then applying the relaxation scheme on the original equations to the oscillatory error components, would be an excellent solution : technique. This is of course the basis of the multilevel methods. I Often slow convergence rates plague eachlevel of the grid refinement, up to I and the fine grid. Coarsening the grid repeatedly untll a level is reached with acceptable relaxation convergence factors leads to the next level of complexity i in the multllevel Each coarse grid would then provide a good to pte next finer grid, reducing the solution effort on the One grids. 24
PAGE 30
The multilevel scheme utlltzed for this study Is based upon an iteration of the VI Cycle, where I a single V Cycle is given by I I and Is defined in (4] as: I Step 1: Relax V1 times on A huh= fh starting with initial guess uh. Step 2: If h = h1 go to Step 4. Otherwise, set rh = fb A huh, eh =l!brb,u2b +O,u2b +MV2b(u2b,f2b). Step 3: Set uh +uh + (19) Step 4.: Relax V 2 times on A huh = fh starting with initial guess ub. (2 0) This suggests an iterative scheme in which successively coarser grtd levels are I employed until a_r,1 acceptable convergence rate for relaxation is found. Then, from I the coarsest grtd level, the algorithm would proceed back through the fine grtd levels, making corrections to the solution at each level. 25
PAGE 31
I I 8. Fast A9apttye Composite Grid Method The above, solution recipe will be augmented by the addition of a composite I grid to further enJlance the procedure. The formation of a composite grid that : I ; 'I utillzes the coarse grtd, the local fine grid patch, and the interface points I I 11 provides a more accurate discretization to guide the solution process. The definition ' and incorporation of the composite grtd into the calculations will be accomplished I via the Fast Adaptive Composite Grid (FAC) method. FAC is an iterative mesh I ,, I refinement applied to the composite grid residual equations, which are used to guide improvements to the classical process. This method, adapted from [5], may I I be described :as follows: I 1) the problem beginning with the coarsest grid down to the level of I grid, obtaining the boundary values and initial guesses for ,: the patches by interpolation from the global grid solution (21) 2) distributive relaxation on the patch = f!! (22) 2 2 2 26
PAGE 32
3) Form the composite grid defect. which consists of residuals from three I regions: the coarsegrid interior points. the patchgrid interior points, ' I and coarsepatchgrid interface points: fh0 on global grid = /1 L1 U1 on interface (23) /hPLh u: in patch. :z :z :z 4) Restrtct the composite grid residuals to the global grid. I 5) $olve the defect equation on global and patch grid levels: I in patch (24) elsewhere. 6) the globalgrid solution to define patch values: for patch interior points for patch boundary points. (25) I The FAC has been coupled with a recursive multilevel VCycle method, I I resulting in a ivecy efficient overall process. 27
PAGE 33
I I Given the .current availability and relatively low cost of multiprocessor computers, a1 natural extension to FAC w.ould allow for the parallel computation of I the various grid levels. The parallel version of FAC has been termed the I ', Asynchronous Fast Adaptive Composite Grid Method (AFAC). AFAC performs each level of refinement simultaneously. The .AFAC method inserts a new grid level between eacll paJt of adjacent grid levels that has the resolution of the coarser grid and the domain of the finer grid. These new grids are then used to elJrninate correction components duplicated by the various grid levels, avoiding the need for a I regressive underrelaxation scheme. The FAC method and an adaptation of the .AFAC method have been used in this thesis to compare results obtained from the periodically grooved channel flow problem. 28
PAGE 34
i ' 9. Fast Ada,ptlve ; Composite Grid Method : I. At in its genesis, the development of the FAC multllevel procedures I I described so (ar could be used to solve many flow problems effectively. However, : I gtven today's and the expectations of future computing environments. a version of coqe that would solve all grtd levels in parallel, with no appreciable i loss of accuracy or convergence rate, could dramatically reduce throughput time for the solution apd v{ould be highly effective in large scale simulation. For instance, I on a hypercube computer, an upperlevel host processor to an AFAC code could be I used to opttniJze the solution time required. The processor would be used to I evaluate the number of points on a grid, thereby obtaining the work required to ' solve each level. Then, based on the CPU availability of the hypercubes constituent processors, the processor could fann out the grtd levels to the machines having the necessary of available computational power. This type of method would I i : help to insure that the most is made of the available computing resources, while finishing the sjmuiatton in a mtnJmum of time. While this seems compelling, is \ I, there a way to such a multilevel scheme and implement it in code so that it l I will run this manner and continue to exhibit fast convergence? A useful.step to be taken before answering this question, which requires ,. I substantial effort to build a truly parallel version of the code, is to create artificially i 29
PAGE 35
., I a parallel version :of the FAC code which would run in single processor mode, a AFAC version (PAFAC). This code could then be used to detennine in th;e convergence factors between the FAC and the PAFAC 'I Implementations. Any deficiencies in the thought process or theoretical I development would then be forced out with a lesser undertaking than that of l a new code. A drawback to creating such a pseudoparallel version stems from being rorced to base' PAFAC upon a scheme which uses previous iteration values for the I I various grids in order to allow for simultaneous processing. The first V cycle lteration through the algorithm is made using the initial trial solution. All :! successive V cycle passes are made using the values obtained from the previous cycle. This prbcedure does not affect the implicit nature of the algorithm, since as updated values produced on the current cycle become available they are included in the The values in current use are obtained simply from the previous lteration, the' newly generated values are held in storage for application in the ; next V cycle. ; I Durtng;the conversion to the PAFAC code, many elements of the FAC algorithm updated to reflect the holdover of values. After the initlalizaUon and residual calcuhltion, the values used in the iterative V Cycle scheme are delayed by one cycle and to begin the next iteration. Following the relaxation phase, the delayed re8fduals are then passed to the restriction routines, where the are: made from the finer grid to the global grid, via the composite grid. 30
PAGE 36
: It is important to note that. by ustng the previous iteration values to generate the PAFAC algorithm. a degradation of the true convergence rate occurs. By forcing the delayed values to be used for calculation 1n the current iteration. some of the I I effectiveness is lost since the corrections being made in the current iteration were generated for. the last iteration. This algorithmic introduction of a sort of retardation Into the solution process serves to slow the overall convergence. However. this really is the basis for our ability to perform AFAC. It enables the simultaneous processing of all refinement levels. The loss of speed is not significant as the following results show. 31
PAGE 37
10. .Relaxation Throughout the previous development, many references have been made to "relaxation Many different types of relaxation techniques have been developed witf varying degrees of success, depending on the application. During I the process of formulating the dlscretlzed model equations, care must be exercised to ensure that the relaxation scheme chosen will be an effective smoother for the discrete equations. Once the incompressible NavlerStokes equations have been put into a dlscretlzed matrix form, It Is often noted that many conventional iterative solvers fail to converge d:Ue to a lack of diagonal dominance in the resulting matrix. Several attempts been made to improve upon this system of equations to allow an ' Iterative solution. Chortn (6) and Patankar and Spalding (7] detail two such I examples. uortunately, such methods often yield inconsistent physics or very slow A multtgrid approach based on distributive relaxation that does not suffer these troubles was introduced by Brandt (8) and Liu (3). This method in effect uses a preconditloner to improve the. characteristics of the discretlzed Implementation of distributive relaxation may be accomplished as follows. 32
PAGE 38
'. I i I 'I First, define tihe new variables wl, w2, and w3 as ,I (26) where I I i ;I 1 Qrh = u. +1vuY L\h, Re 8; is the backwd Euler operator, : :: 01 and oy central difference operators, and L\h is a Laplacian. Next, define the discretlzed matrix form of the NavterStokes equations as I (27) Applying (26) ,to (27) yields I : 0 [;]=[! 0 8,] [I 0 8] [ wl] :I 8: w2 ' Qrh oy o 1 (28) IQ' rh' 'I 8 .. I 8 I 0 0 0 0 Qrh w3 : y' y 33
PAGE 39
Equation (27) is then related to the linear system .. ] [wl] [fl] w2 = f2 &1 w3 f3 (29) Even though the system in (29) only applies in special cases, it is easily solved by conventional methods. and may be used as a basis to develop an I applicable method. Adapting this method to the flow problem results in the I following Begin1with the system I I UE ;Up + _V.:.:,N__V_.:...P =S. m (30) 34
PAGE 40
Then perform the distributive relaxation according to the following recipe: : 1) Freeze P Perform point GaussSeidel relaxation on the momentum equations (30) I and to obtain a new u and v 2) I For each grid cell (Figure 10.1) modify the velocities to \ satisfyithe continuity equation Next, modify the pressures so that the new velocities satisfy the momentum equations I The velocity are calculated by detenn1n1ng 8 and f3 such that (uE + /35)/35) + (vN + o)(v PD)= S Ax: !:J.y m I O=[sm ( :xu,+ vN + :J f3 =fly. llx 35 (31)
PAGE 41
I 0 Figure 10.1 Distributive Relaxation Now that 8 arid f3 have been detennined, let UE f.t UE + /38 up f.: uP VN f.i VN ,+ 8 I Vp f..1 Vp8. The pressure f1 and the neighboring pressures are updated according to I I, 36 (32) (33)
PAGE 42
At this point, the velocities are the old values and 6 1s the correction from (31). Choosing 5P ++:tuo '[ 3 1 ( 2 2 )] p : 2at Re axz ayz f3. (34) allows an eaSy determination of c5Pw and 8!>8 from (33): OP OP A f31lx8( 1 w P P Rellx2 2/lx &s ::r&pBp/3llyD( 1 2 Reily 2/ly (35) The detennination of 5P6 and 5P N 1s analogous. Due tc;> the selection of a small at. the pressure correction term 5PP is much I larger than the other pressure corrections c5P 6 c5P w c5P N and 8!>8 for this problem. for this algorithm, the smaller pressure correction terms are I assumed to zero, while the main pressure correction remains unaffected: aPE=, aPw = OPN = OPs = 0 '[ 3 1 ( 2 2 )] BP = ++f3tu.o P 2/lt Re Ax2 lly2 (36) 37
PAGE 43
titil.tzt;ng this form of relaxation has allowed the solution of the N avierStokes equations for the flow In a grooved channel to proceed without the usual convergence problems or diversions from the flow physics usually associated with this class of problem. 38
PAGE 44
' 11. Computational Results The differences in convergence between the PAFAC algorithm and the I standard FAC algorithm are demonstrated through a comparison of the convergence ' factors between iterations of the two methods. Basic multlgrid cycles, used in both I FAC and PAFAC, have been shown [9) to converge at a linear rate that is I independent of size, with typical convergence factors of 0.2 or less. The parallel nature ofPAFAC degredates the convergence characteristics from the 0.2 value. In fad, asymptotic convergence factors [9) are roughly the square of ,i those obtainetl for:AFAC. Since the PAFAC code used here emulates anAFAC code, i I similar results are expected. A Poiseuille flow, discussed previously, will be used for I I I initialization. :A variation study will be performed for a range of Reynolds numbers. The Reynolds number parameter Re was chosen as a critical parameter due to its I substantial effect the overall nature of fluid flow. One may interpret the I definition of Re in '(6) on a slightly more heuristic level as ' Re = : IneJ;iial Forces !Viscous Forces 39 (37)
PAGE 45
I As the inertial forces increase and dominate the visCous forces. the. fluid flow I becomes more turbulent. A range of Reynolds numbers will then present the opportunity to test fully the algorithm. In orderto compare the FAC and PAFAC codes, some basic ground rules I were establisped. The grid structure for each method was identical, with the same initial .. The global mesh in each case was 10 X 18, with a patch I refinement of 18 X 18. Since the multilevel schemes have been based upon the I I residual error to estimate corrections to the solution, rather than by estfmatlng the solution itself, the convergence factors will be calculated using the residuals I I I ' determined when exiting the fine patch grid at the end of a cycle: r Convergence factor = l:!:!.. r I In both codeS the nann used to calculate the residuals is the C2 norm, given by: This norm was chosen due to the ease of calculation given the grid structure I employed. 40 (38)
PAGE 46
' Each :method has been Umited to one cycle per iteration, except, for the iteration to make the convergence factor comparison, iteration number 10. In iteration 10, the FAC version of the code performs 5 cycles, while the PAFAC I version of th
PAGE 47
I Reynolds number Cycle number 1 I 2 3: 4 5' 6: 7 81 Reynolds number Cycle number 1 : 2 3.! 4! 5 6! 7 : 8 I 225 1000 1.15E04 2.31E04 3.03E05 1.29E04 1.50E05 4.34E05 8.49E06 1.82E05 1.86E06 7.62E06 6.71 E07 2.91E06 1.34E07 3.49E07 7.94E08 1.73E07 Table 11.3 AFAC u Residuals 225 1000 9.22E05 5.71E05 9.76E06 1.65E05 6.59E06 8.72E06 4.16E06 5.77E06 1.04E06 2.73E06 2.58E07 7.00E07 8.44E08 1.42E07 3.88E08 7.31 E08 Table 11.4 AFAC v Residuals 42 7500 75000 3.79E04 4.1 OE04 2.60E04 2.92E04 8.57E05 9.62E05 2.69E05 2.86E05 1.41 E05 1.55E05 5.32E06 5.81 E06 5.52E07 5.92E07 2.74E07 3.08E07 7500 75000 8.18E05 8.84E05 3.02E05 3.36E05 1.31 E05 1.42E05 7.19E06 7.38E06 4.08E06 4.31 E06 1.16E06 1.25E06 2.02E07 2.14E07 9.94E08 1.05E07
PAGE 48
The convergence factors calculated for the FAC and PAFAC methods are presented tn Tables 11.'511.8, for FAC u and v and PAFAC u and v respectively. Reynolds number 225 1000 7500 75000 factor cycle 12 0.063 0.043 0.044 0.045 23 0.201 0.114 0.131 0.149 34 0.339 0.239 0.203 0.199 45 0.198 0.253 0.610 0.651 I AveraQe conv.erqeflce 0.171 0.131 0.163 0.171 factor by root method Table 11.5 FAC u Convergence Factors Reynolds number, 225 1000 7500 75000 I Converqence factor cycle 12 0.190 0.143 0.108 0.105 23, 0.212 0.224 0.257 0.265 34: 0.278 0.238 0.242 0.243 4Si 0.240 0.252 0.350 0.375 I Averaqe converqehce 0.228 0.210 0.220 0.224 factor by root method Table 11.6 FAC v Convergence Factors i 43
PAGE 49
Reynolds number 225 1000 7500 75000 Converqence factor cycle numbers 12: 0.264 0.557 0.687 0.711 23; 0.495 0.337 0.329 0.330 34 0.567 0.418 0.313 0.297 451 0.219 0.419 0.525 0.541 56: 0.360 0.382 0.377 0.375 67 0.200 0.120 0.104 0.102 78' 0.592 0.495 0.497 0.521 Averaqe 0.354 0.358 0.356 0.358 factor by root method Table 11.7 AFAC u Convergence Factors Reynolds number. 225 1000 7500 75000 Converqence. factor cycle numbers 12: 0.106 0.289 0.369 0.380 23:' 0.675 0.527 0.433 0.422 34: : 0.631 0.662 0.550 0.520 45; 0.250 0.474 0.567 0.584 56: 0.249 0.256 0.286 0.291 67: 0.327 0.203 0.173 0.171 781 0.460 0.516 0.493 0.491 Averaqe convrrqef)Ce 0.329 0.386 0.383 0.382 factor by root method Table 11.8 AFAC v Convergence Factors 44
PAGE 50
The approximate averages given In Tables 11.511.8 are computed using a root I average method, as I : kf5:. average = Tl where n is the total number of cycles performed In the iteration and k=(n1). For both FAC and PAFAC, convergence was obtained for all Reynolds (39) numbers used. The FAC convergence results show good agreement with the 0.2 convergence factor achieved by McCormick [9] for simpler model problems. Moreover, the PAFAC convergence results follow the trend reported by McCormick [9], being roughly equal to the square root of the FAC convergence factors. 45
PAGE 51
12. Conclusions Ustng. the flow through a grooved charmel as a model, the overall I convergence of two versions of a multilevel iterative scheme, FAC and PAFAC, have;been compared. It has been found that the FAC method had convergence factors of 0.23 or better, while the PAFAC method had convergence factors of 0.39 or better. These convergence factors appear to support the results i drawn by McConl)ick (4), (9), for both FAC and AFAC type codes on stmpler model problems. The results presented here are especially sfgnificant. The behavior of FAC and AFAC for model problems has been shown to extend to more complex applications that include systems, nonlinearity, and changing character (i.e. high Reynolds number). i ExamfnJng the numbers of cycles performed for the two methods, it may be concluded approximately 8 PAFAC cycles are equivalent to 5 FAC cycles. This is a marginal increase in the total amount of work required, and is more than I compensated for by the full parallellzability of AFAC. I 46
PAGE 52
References I (1] Llu,c: . Llu,z . and McConntck,S.F. 1991. MultaevelAdapttve Methods for I Flow tn Grooved Channels, J. of Computational and Applied Math, i 38:283295. : [2] Llu,Cl,Llu.z . and McConnick,S.F. 1991. Multtgrtd.Methodsfor Plow Transition tn.:a Planar Comput.Phys.Comm. 65: 188200. [3) Liu,C: 1989. Multilevel Adaptive Methods tn Computational Fluid Dynamics. I I Ph.D thesis, University of Colorado at Denver. I i [4) McCo,anck.S.F. 1987. Multtgrid.Methods, SIAM. Phil. (5] Liu,c; 1991. Multilevel Adaptive Gridsfor flow Transition tn a Planar Sctence.Report, Ecodynamtcs Research Associates. Inc., Denver. [6] Chorth.A.J. 1967. A Numeriail. Methtxlfor Solving Incompressible Viscous i Plow Problenis, J.Comput.Phys. 2:1226. [7] . and Spaldtng,D.B. 1972. A Calculation Procedure for Heat, Mass and MJnentum Transfer tn 3Dimensional Parabolic Plows, Int.J .Heat Mass _I I Transfer 15:1!7871806. I (8] Brandt.A '.1984. Multtgrtd Techniques: GuidewtthAppltcation to F1.uid i Dynamics, GMD Studie, GMD. St. Augustine. I .. 'I [9] Mccorrruck.S.F. 1989. MulttlevelAdapttve Methodsfor Partial Differential 'I Equations. Sl.f\M. Phil. 47
