Citation

## Material Information

Title:
Multigrid and finite difference methods for flow transition in a planar channel
Creator:
Liu, Zhining
Place of Publication:
Denver, CO
Publisher:
Publication Date:
Language:
English
Physical Description:
x, 58 leaves : illustrations ; 29 cm

## Subjects

Subjects / Keywords:
Transition flow ( lcsh )
Multigrid methods (Numerical analysis) ( lcsh )
Multigrid methods (Numerical analysis) ( fast )
Transition flow ( fast )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

## Notes

Bibliography:
Includes bibliographical references.
Thesis:
Submitted in partial fulfillment of the requirements for the degree, Master of Science, Department of Mathematical and Statistical Sciences
Statement of Responsibility:
by Zhining Liu.

## Record Information

Source Institution:
Holding Location:
|Auraria Library
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
25983332 ( OCLC )
ocm25983332
Classification:
LD1190.L622 1991m .L58 ( lcc )

Full Text
MU;LTILGRID AND FINITE DIFFERENCE METHODS
1 j
! FOR FLOW TRANSITION
IN A PLANAR CHANNEL
1 by
Zhining Liu
;B.S., University of Science and Technology of China, 1986
: ^ '
j. A thesis submitted to the
! Faculty of the Graduate School of the
' ' I
| | University of Colorado in partial fulfillment
I
of the requirements for the degree of
I
Master of Science
Applied Mathematics

This thesis for the Master of Science
degree by
Zhining Liu
has been approved for the
Department of
Mathematics
by
Steve McCormick
Thomas Russell
2-3 w V
Date

Zhining Liu (M.S., Applied Mathematics)
i
Multigrid anid Finite Difference Methods for Flow Transition in a Planar Channel
Thesis direc
ed by Professor Steve McCormick
The flow
transition problem has been central to the theory of fluid motion for
over a century. Developing an understanding of this procedure will be very helpful
to the research of both theoretical and applied fluid dynamics. A computational
study of the temporal stability of planar Poiseuille flow in a channel is presented.
Using fully-jimplicit time steps on a staggered grid, we develop a second-order
scheme, a second-order in the streamwise direction and fourth-order in the normal
i
direction scheme:, and a full fourth-order scheme. A very efficient time-marching
procedure usingj a multigrid solver is shown to produce very accurate results.
The computational results show that the full fourth-order scheme gives very good
agreement in both amplitude and phase with the linear theory obtained by solving
the Orr-Som!merfeld equations.
The form and content of this abstract are approved* I recommend its publication.
11 /// / III //
Signed
Steve McCormick

To my parents

CONTENTS
List of Figures................................................... vii
List of Tables..................................................... ix
Acknowledgements.................................................... x
1 INTRODUCTION
2 GOVERNING EQUATIONS
FOR 2!
D PLANAR CHANNEL FLOW
Governing Equations and Boundary Conditions
Governing Equations in Perturbation Form
Second-Order Fully-Implicit Finite Difference Scheme
i
Cross-Channel Fourth-Order Finite Differences ....
Fourth-'Order Fully-Implicit Scheme
Interior Points
Boundary Modification
3 ALGORITHM
Distributive' Relaxation
Multigrid Algorithm .
1
4
4
5
7
14
17
17
25
28
28
31
Fine Grid Relaxation and Residual Evaluation
32

Coarse!Grid Relaxation
Interpolation of Correction to Fine Grid
4 TEST CASES................
I
I.
. !
Initial Small1 Disturbances
' J
Computational Results .
Concluding: Remarks
5 CONG
i
.USIONS AND OPEN QUESTIONS
Synopsis of the Thesis
'I' :!
Open Questions . .
BIBLIOGRAPHY
vi
33
35
37
37
43
52
53
53
54
56

2
5
7
8
9
10
12
13
19
20
21
23
26
30
32
34
36
LIST OF FIGURES
Reynolds transition experiment................................
Planar channel flow............................................
I
Staggered grid.................................................
i
, I
Grid neighbors for x-momentum equation.........................
Neighbor points of up..........................................
I
'I
Neighbor points of vp..........................................
Grid neighbors for continuity equation.........................
. I
Solid wall boundary condition. . ............................
I
Neighbor points of up for u-equation...........................
Fourth-order approximation for vatau point (interior points).
Neighbor points of vp for v-equation...........................
Neighbor points of Pp for continuity equation..................
!
Approximation of v on solid wall boundary points...............
Distributive relaxation.............
Two-level staggered grid............
l
Restriction for R,..................
Bilinear interpolation for u, v and P.

4.1.
4.1.
4.1.
4.2
4.3
4.4
4.5
4.6
4.7
4.8
4.9
(a) Orr-Sommerfeld eigenfunction
(b) Orr-Sommerfeld eigenfunction
(c) Orr-Sommerfeld eigenfunction
11
Convergence history for
(Re = 7500, mesh: 34 x 130 )
(a) u, (b) v, (c) mass........
with Re =2000.
with Re =7500.
with Re =20000
Disturbance energy increase
(second-order fully-implicit, 3 periods)
Disturbance energy increase
(cross-channel fourth-order, 3 periods).
Disturbance energy increase
(fourth-order fully-implicit, 3 periods).
(a) Re- = 2000, (b) Re = 7500, (c) Re = 20000
. |
Contours of disturbance stream function
(second-order fully-implicit, one period). . .
Contours of disturbance stream function
(cross-channel fourth-order, one period).
Contours of disturbance stream function
(fourth-order fully-implicit, one period).
Cpmparison of contours of disturbance stream function
(Re = 7500, time=20 periods)
(a|) linear theoretical result, (b) fourth-order fully-implicit,
(c) cross-channel fourth-order, (d) second-order fully-implicit. .
4.10 Comparison of different schemes with disturbance energy.

LIST OF TABLES
4.1 Temporal Frequency for Different Re.....................
j 1
4.2 Comparison of the Increasing Rate of Disturbance Energy
(in) 20 Periods)........................................
I
4.3 Comparison of the Maximum Stream Function Values
arid Phase Lagging (after 20 Periods)...................
! i '

ACKNOWLEDGEMENTS
I would jlike to take this opportunity to thank my advisor, Professor Steve
i
McCormick^ for the guidance and encouragement that he provided during all
'
phases of the preparation of this thesis. It has been a privilege to work with him,
I
|
and I hope; that I have acquired some of his insight in the process. Thanks go to
1 j . |
Professor Ghaoqun Liu, who provided me with close supervision and who guided
me step by'step through algorithm development. Thanks also go to the other two
! i
members oif my1 committee, John Ruge and Tom Russell, for their contribution
I'!
to my development.
: i
I would I also like to thank my parents and my friends for their unflagging
i
support. Their encouragement kept me going, and I am in their debt. It is
impossible| jto list them all, but a few deserve special mention: Debbie Beltz,
: i
Gaoming Yang and Jonathan Hazen.
'!
I am also indebted to NASA for their support of this work under Research
, I
Contract N!ASlrl9016.
I

CHAPTER 1
INTRODUCTION
The transition process from laminar to turbulent flow in a shear layer has been
i
the subject of intense study for over a century, but it still remains a challenging
and important area for research. The process of transition is very complicated.
Early work was initiated by Reynolds [1], Rayleigh [2] and others. Orr [3] and
Sommerfeld [4] later reformulated Rayleighs theory for viscous incompressible
fluids. Toll'mien [5] and Schlichting [6] successfully predicted the structure of
Orr-Sommerfeld modes in noninflectional flows. However, current technology
has yet to, provide an understanding of the first classic transition experiment
. i
!
completed by Reynolds more than one hundred years ago. He showed that a thin
filament of dye in water flowing along a transparent glass tube would remain well
defined and straight if the Reynolds number of the flow, fZj, was less than about
2,500; but for higher Reynolds numbers, the filament rapidly became convoluted,
thickened and broke up, the dye dispersing throughout the tube (see Figure 1 (a)
' I
and (b) ). Here, R^ = ud/u, where d is the diameter of the tube and u is the
centerline velocity of the mean flow. But by using spark photography to capture
' I
l
details of the flow, we can see that in the fully mixed region of the tube there
are many vortex-shaped streamlines. The vortex moves 3-dimensionally in the

water tank
transparent glass tube
(a)
(b)

| (c)
Figure 1.1: Reynolds transition experiment.
i
spatial direction with a frequency in the several thousands (see Figure 1 (c)).
i
As currently understood, the transition to turbulence is, in fact, a very com-
i
plicated 3-:D, nonlinear and multistage process, which begins with a sequence
of instabilities on a succession of more and more complicated basic flows (c.f.,
Bayly and lOrszag [7]). Loosely speaking, the fluid flow will be laminar when the
Reynolds number is small, and become turbulent when it is high, which means
that as the .Reynolds number increases to some value the laminar flow will begin
to become, turbulent. People have established a so-called small disturbance the-
ory to determine what happens with respect to time when a small disturbance is
added to the steady laminar flow. This theory can tell us what kind of velocity
profile is unstable, which mode of the oscillation will be amplified the fastest,
and how to change the flow parameters to delay transition. Unfortunately, it
1 2

tells us almost nothing about the nonlinear instability and the vortex breaking
stage, which are both very important to the study of transition.
' i
Numerical simulation of flow transition in a channel has been successfully done
I
by severalresearchers (c.f., Zahn, Toomre, Spiegel and Gough [8] and Orszag [9]),
most of whom use spectral methods with explicit time-marching and non-
staggered grids (c.f., Zang and Hussaini [10]). The difficulties with numerical
simulation for transition result from the long-time evolution and the sensitivities
to small changes, requiring high accuracy and efficiency from the numerical meth-
j
ods. Although it is recognized that the staggered grid (c.f., Harlow and Welsh
[11]) and implicit time-marching approaches have advantages, there has really
been little:attempt to apply them, especially in the field of transition.
The pujpose of this thesis is to develop an accurate finite difference scheme
on a staggered grid with implicit time-marching. Since we must solve a spatial
, I
system at each time step, a full approximation multigrid algorithm (Brandt [12])
is applied tjo accelerate the convergence process.
I
The test problem is a 2-D planar Poiseuille channel flow. Different Reynolds
numbers were chosen to test the accuracy of our schemes, and good results were
; I
obtained by our final fourth-order fully-implicit scheme.

I
CHAPTER 2
I
GOVERNING EQUATIONS FOR 2-D PLANAR CHANNEL FLOW
I
I
I
I
In this chapter, the governing equations for planar Poiseuille flow in analytic
and discretized forms are presented. Several types of schemes are developed for
the temporal stability problem.
I Governing Equations and Boundary Conditions
I
L . .
A two-dimensional, time-dependent, incompressible Navier-Stokes equations
can be used as the governing equations for planar channel flow:
du duu duv 1 d2u d2u. dP
dt ^ dx ^ dy Re dx2 dy2 dx
dv duv dvv 1 d2v d2v dP
dt^ dx dy Re dx2 dy2 dy
du dv
dx ^ dy
(2.1)
(2.2)
(2.3)
where u and v are velocity components in the x- and y- directions, respectively,
. I
P is the pressure and Re is the Reynolds number based on the centerline velocity,
U, of mean, flow, the channel half-width, h, and the viscosity parameter, v.
" ' * = 25. (2.4)
i v
The boundary conditions on the solid wall are
u(x, 1, t) = v(x, 1, t) = 0,

u(x, 1, t) = v(x, 1, t) = 0. (2-5)
Because we 'choose to use a staggered grid, boundary conditions for P are not
1 i
needed on jthe solid wall. For simplicity, we only consider the temporal instability
]' I '
for channel] flow, and assume periodic boundary conditions in the x-direction:
j u(x,y,t) = u(x + L,y,t),
l
j v(x,y,t) = v(x + L,y,t),
1 P(x,y,t) = P(x + L,y,t) + -^-eL, (2.6)
where L is the wavelength of a specified mode of disturbance.
\
: ]
i j
i j Governing Equations in Perturbation Form
I I
Poiseuille flow is determined as the solution functions uq,vq and Pq of (2.1)-
(2.3), with jboundary conditions (2.5)-(2.6), that satisfy the initial conditions:
i
o(*,y.o) = l j/2,
II
. i
ii ' to(*.!fi0) =0,

2
Po(x,y,0) = -X + C.
(2.7)
i
Here, C is constant. This will be called an undisturbed solution. Since u0,v0 and
Po are independent of time, we have
u0(x,y,t) = u<,(x,ji,0) = 1 ~ y3,
(2.8)
and similarly for v0 and P0.
I
If we impose some disturbance at t = 0 on Poiseuille flow, it will generally
i ,
, j i
become time-dependent. For such a disturbed flow, in order to reduce round-
off error in the computational solution, we write the primitive variables in the
I
following perturbation form:
u = u0 + u1,
v = v0 + v\
P = Po + P'. (2.9)
Substituting (2.9) into (2.1), (2.2) and (2.3), and using the fact that v0 = 0, we
may then rewrite the governing equations in the following perturbation form:
du du0u du0 duu duv 1 d2u d2u dP
dt "^"i dx dy ^ dx dy Re dx2 ^ dy2 dx
dv du0v duv dvv 1 d2v d2v dP
dt dx ^ dx ^ dy Re dx2 ^ dy2 dy
du dv
dx ^ dy
(2.10)
(2.11)
(2.12)
Here, for simplicity, we drop the prime notation, so that u, v and P actually now
l !
stand for the perturbation variables u,v' and P'. The boundary conditions may
6

be rewritten as
u(s,1,Â£) = v(x,l,t) = 0,
u(x,l,t) = v(x,l,t) = 0,
u(x,y,t) = u(x + L,y,t),
v(x,y,t) = v(x + L,y,t),
P(x,y,t) = P(x + L,y,t).
(2.13)
Second-Order Fully-Implicit Finite Difference Scheme
i'
; | j
We use; a uniform staggered grid for our problem (see Figure 2.2). A second-
order backward; Euler finite difference in time and a second-order central differ-
I Figure 2.2. Staggered grid.
First consider the x-momentum equation (2.10). The related spatial com-
1
ponents for jare shown in Figure 2.3. With superscript n used to denote u
; 7

at the nth time increment, ito the undisturbed solution of the channel flow in
I
xdirection, v the yvelocity component defined at a u point, then we discretize
equation (2'JlO) on the above staggered grids using a fully-implicit finite-difference
scheme:
3 4. - 2v,-v"+ +
l 2At 0 2As Vl x'3
+(
+(
1l7l+1 _L 7|n+1 7 ?+l,j!+ UiJ Ui+l,j + Uij + ui,j Ui-hi +

)/A*
7,^+1.i _L 7,?71 7/?1"1 _L 7/TlT1 7I"-1"1 _L ll"?1 _L 7/r
Vl,j+l + vi,j+1 ^,j+l + _ ^-i,j + vi,j m Ui,j-l + Ui,j WA
o o o <\ H y
n+1 ,n+l
n+1 n+1 I n+1 n+1
,"+i
-k
<}j Kt1 + | A2 Ay2
p^-y i _

As
= 0.
(2.14)
i Figure 2.3. Grid neighbors for x-momentum equation.
|
' |
For simplicity, we may rewrite (2.14) in the general form (see Figure 2.4)
Apup dr Awuw + Anun + AsUg Apup CpPp + CwPw = Suy
8
(2.15)

where
Up
uw
1lN = ui,j+l,
US =
Up = Ui,j ,
Pp ; - Pi,:,
Pw = Pi-U-
' l
I ' i i i i , UN
i i ?w 0 Pw Up O Pp Up
i | i Us
Figure 2.4. Neighbor points of up.
Define
i
i
i'
i
uB =
U, =
=
0.5(Sij + )
*)>
0-5(."-+U+i + &M.
9

I
v, = 0.5(vr_Vj+ This yields
Ae |
I
I
Aw ;
Re Ax2 2Ax
(^0i,j A ^e)j
1 1 /
+ 7m\uoi,j + )i
-^AT I =
Re Ax2 2Ax
1 1
! Re Ay2 2A y
v
As^r
Ap =
1 1
------1- TTTVs,
Cp
Su
Re Ay2 2 Ay
3 1 1 . , 1 , . 2,1 1 ,
>2 At + 2Ax,U + 2Ay Vn + Re Ax' + Ay^
^-L
!; Wl + 3i + + Â£j+.)(-2w) + 2toU'F-
j Figure 2.5. Neighbor points of Vp.
I
Similarly, we write the difference equation for (2.11) in the general form (Fig-
j j
ure 2.5) j
. I
BeVe\-\- Bw^w + BffVn + Bsvs Bpvp DpPp + DsPs = Sv, (2.16)
'I 10

where
Ve
Vw =
Vff vi,j+lt
Vs = Vi,j-1>
Vp Vi,h
Pp = Pi,:,
Ps = Pi,:-1-
With
we get
Be V
Bw !=
Bn ^
d LL
Bs =p
; i
!
ue =
uw =
V =
V. =
0.5(<+ +
0-5i1 +
o-ss,+<;>,
0.5ft + Â£>),
1
1
Re Ax2 2Ax
(Oi,j "f~ ^e)j
1 1 /
"I" ^ "f" ^uijj
ReAx2 2Ax
1 1
ReAy2 2Ay
1 1
+ r~^v
ReAy2 2Ay
.,-31 1 , , 1 , * 2 1 1
Br 2Tt + ^U'-U)+2A;{v--V-) + Te(A*i+W
:! li
),

_ I 1
Dp i = Dw = ,
j = *rlz*?L.
1 2A t
Figure 2.6. Grid neighbors for continuity equation.
The continuity equation may be discretized as (Figure 2.6)
Feue Fpup + GpfVpf GpVp = Sm,
where
fE = irp= aJ
Gk = Gs = SJ
Sm = 0.
(2.17)
The boundary is modified as follows. At the periodic boundary, we write
= uni-l ,j)
vu = vni-l,ji
Plj = Pni-\,ii
uni,j =
12

Vni,j = V2,j,
Pni,j = Pl,ji
(2.18)
where ni is the number of grid points in the x-direction (for one wave length).
On the solid wall boundary, we use a modification of the x-momentum equa-
|
tion (Figure 2.7). For example, for j = 2 we use the approximations
d2u
4 Uff up up 0 4 1 4
dy2 3Ay Ay \Ay 3 Ay2** AyzUP'
! duv 1 ^Vij+i +Vj-ij+i uN + up ^
(2.19)
(2.20)
dy Ay 2 2
i i
We simply, modify the appropriate coefficients in (2.14) at the boundary points.
For example, for j = 2 we make the changes
I
1 1
Anew a old
N ~An+ ReZAtf'
Apew = Afd +
l__2_
Re Ay2
(2.21)
UN
Up

I u = 0
I Figure 2.7. Solid wall boundary condition.
For j = nj 1, where nj is the number of grid points in the y-direction, we use
i A^.ew = A}d ________-_
,! As 5 ' Re 3 Ay2
i new A old , ^ ^
lp ~Ap +Tuw'
13
(2.22)

Cross-Channel Fourth-Order Finite Differences
i
i
Since th'e solution tends to change more rapidly in the y-direction, either re-
i '
1
fined gridsjor higher order differences in the y-direction can be used to balance
overall accuracy. One of the methods we study in this paper uses fourth-order
central differences in the y-direction for interior points only, while maintaining
second-order central differences for boundary points and for interior points in
the x-direction. We use a fully-implicit scheme, since it is generally more accu-
i '
i
rate than semi-implicit or explicit methods. For simplicity, we continue to use
| !
j
second-order differences for the continuity equation and the pressure terms in
' l
both momjentum equations. The general form for the fully-implicit scheme in
'll
this case can be written as follows:
Aeue -f Awiiw + AffjfUifff + Anun + As us + Assuss
i
I
! Apup + C\nrP\v CpPp = Su,
I Beve + BwVw + BifffVpfif + BpfVff + Bsvs + Bss^ss
; i
Bpvp -f DsPs DpPp = Svy
i
Feue FpUp + GffVff Gpup = Stti.
Here, similar to. the previous section,
uNN = uij+2,
Uss = U{j _2y
14
(2.23)
(2.24)
(2.25)

1
VNN = Vij+2,
Vss = Vi,i-2-
In (2.23), if we define
then
i i
ue 0.5t- + <.),
uw 0.5|> + ),
Un 0.5t> +
V-. = 0.5t + ?+_>,),
vn = 0-5Ka+t,?_YJ+1),
^nn O-SMJi + '%),
- 0.5t* + ?.),
V-
- o.25(rt+<;+!+;
n+1
,.n+l
Ae
AW
Ann
An
As
Ass
1
Re Ax2 2Ax
(uOi,j + Ue),
! 1
1
+ 0 A (UOi,j + Uw),
Re Ax2 2Ax
i ^ 1 1 J_J_
j Re 12Ay2 48 Ay
| i 4 ii. .
Re 3Ay2 48 Ay 27Un)
I4 1 1 /nn
+ -7 (27v, v),
Re 3 Ay2 48 Ay
' 1 1 J_J_
Re 12Ay2 48 Ay
15
i

Ap =
3 1 . , 27 , . 1 2 5 .
+ - u") + - Vs) + s-(-r-s + tit-;),
Cw, = \ Cp =
2 At 2Ax
1
48Ay
ReK Ax2 2 A y2
Ax'
4 id* + u? 1
Similarly,
define
2At
Ve = o-SMj1 +
Vw = 0.5+ + v&'j),
Vnn = + &!.),
vn = o-W+O
V. = 0.5f +
= 0.5K/-2 + "A).
= 0-5m + tfis-i)-
= o-sKi' + C)-
Then, for equation (2.24), we have
Bp 1= i
Re Ax2 2Ax
(Uoij + ue),
Bw \ =
1 1 ,
+ 7rr~\uoi,j + u)i
; Re Ax2 2 Ax
Bnn != I
1 1
A t n ^nn,
i?el2Ay2 48 Ay
Bn |^e3Ay2 + 48Ay(Vnn 27Vn)'
; ; 1 4 1 1 ,
Bs '{= r> o a + --r-(2Tv, -u,,),
-#S5 =
Re 2 Ay2 48 Ay
1 1 J_J_
ile 12Ay2 48 Ay
16

'! ! 3 1 . , 27 , , 1.2 5 ,
Bp ~ + 7m(v" ~ v) + TcHTtK v.) + (,-t-t + rr-r),
2A t 2 Ax
Ds .r Pp = Si'
48 A y
n1
Sv '=
:-4v- 4- v?~
2Ai
ReKAx2 2 Ay2'
There is noj change to the continuity equation, so
I- I
= Fp=Ai
Gk = Gf = 7sj
Sm = 0.
Since the second-order scheme is maintained at the solid boundary points, the
boundary modifications are the same as the previous section.
i :
I i
. i Fourth-Order Fully-Implicit Scheme
' [. ~ :
II :
, i i
Since this scheme is quite complicated, we treat the interior points and boundary
I; i
points separately.
' |
' I i
' i !
i, 1
Interior pjoints
__ i
Though ithe above schemes are fairly accurate, especially the cross-channel
' !
i ;
fourth-order fschjeme, they are inadequate for simulating transition because of the
sensitivity of suph flows to small changes and the need for long-time evolution.
For these reasons, we developed a so-called fourth-order fully-implicit scheme on
ill. 1 . m
a staggered grid1, with the aim of resolving the problems of lagging of the phase
! 17
1' 1
: I :

I i
and damping of the rate of disturbance energy increase that we observed in the
above schemes. iSince we only consider the temporal problem at this time, which
i I
l ;
means that 5 we only consider the instability with respect to time, there is no
I i
II
difficulty iniapplying periodic boundary conditions in the x-direction.
! i
The general form of the fourth-order fully-implicit finite-difference scheme can
:
be written as:
11
i :
' | r
11 ; Aeeuee + AeV>e + Ayfuw + Awwuww +
!
| ^ AjfifUfjfj + AprUff + Asus + Assuss Apup +
! ; CwwPwW + CwPw + CnPn CpPp = Su,
: I
I
: i
' j i Bee^ee + Beve -f Bw^w + Bwwvww +
I |
i; | B-un-Vnn + BprVff + Bsvs + Bssvss Bpvp +
DssPss + DsPs + BnPn DpPp = Sv,
Feeuee + Feue + Fwuw FpUp +
Gnnvnn + Gnvn + Gsvs ~ Gpup = Sm.
(2.26)
(2.27)
(2.28)
Here,
, i i
, i ;
UEE = ut+2,jj
UWW ~ ui-2,j,
Vi+2J,
Pi-2,31
VEE =
Vww -
Pww =
18

order approximation for v at the midpoints of vertical cell edges. The fourth-
i I
I
order approximation we use here involves the points shown in Figure 2.9. The
l
formula is i j

= +ajn) -
(, .n+l I ,,Tl+l 1 n+l I n+l
(Vt-2,j-X + Vx-2,j+2 + U+l,j-l + Vi+l,j+2))-
n+l
n+1
,.n+l
According to the above approximation, we obtain the following coefficients
i
i
for the uequation:
Aee ;?=
|
Ae
Aw ,=F
;i
i
i
i4lvw
1 1 , A \
+ + 4u0i.i),
. I 12ReAx2 48 Ax
+ 48^(u "27- ~32uoij)
+ 48^(27 + 32uijh
i 1 1
12i2eAx2 48Ax
(Uww t 41X01,7 )i
20

Ann ==
i
An =
* t
-4^5 =f
+
1 1
Re 12Ay2 48 Ay
I1 4 1 1 / V
iZe3Ay2 + 48 A~ 27u")>
14 1 1
i*e3Ay2 + 48 Ay 27v* V"
11 11
iZe 12Ay2 48 Ay

Ce ==
24 As
Cw = Cp =
CVw =
27
24 As
24 As
Su
-4 <, + 7^71
-----^-r--------2y,Vj
2AÂ£ ^
o - 1 O
o 'U'WW - 1 PP U-w j Vp o ue
o . f Ps ^t>s o
o 1 PSS^V3S t VSS o
w,

Figure 2.10. Neighbor points of up for u-equation.
21

Similarly, for the v-equation (2.27) (see Figure 2.10), let
v.e =
Uiiillt
0-5(u+i j + i++2j),
,n+l
u.
Wllltrt ------
o-s(nfij+
0.5(u{*2j -f i+2,j-i)>
0.5(uJt11J- + "J^j-x)-
This yields ;
Bee =f
Be :4=
Bw !=j=
Bww =
Bnn =
Bn 4=
Bs =
Bss =
BP =
I
Dn
Ds
Bss
Sv
1
12ReAx2 48Ax
4 1
+ 7TXZ(U + 4uov)>
3ReAx2 48A
4 1
+
ZReAx2 48A
1 1
, 12i2eA2 48 A
: 1 1
(uee 27ue 32uoiJ),
(27 u
W 'U'ww f~ SSlioijj)?
(^uitu I- 4uof,i);
1 1
I n a ^nn j
! Re 12Ay2 48 Ay
i 4 ii. ,
+ _ 27u),
Re 3 Ay2 48 Ay
i 4 1 1 ,0, .
Re3Ay2 + 48Ay(27v>~v^
; 1. 1 11
~Re 12Ay2 ~ 48 Ay*"
SS + 4^ ~ ^ + 5SJ+ +
1
2fZe'A2 A y2
24 Ay
= DP =
27
24 A
' 24 Ay
Ki + <7*
2At
22

VNN
Uw
Up
Pp
\Vp
UE
I Vs

figure 2.11. Neighbor points of Pp for continuity equation.
The continuity equation, (2.28), becomes fourth-order with the following coeffi-
cients (see Figure 2.11):
Fee
Fe
Fw
Gnn
Gn
Gs
Srn
1
24 A
= FP =
27
P 24 A
1
24 A
1
~ 24Ay
= Gp
27
24 Ay
24 Ay
0.
The complete discrete system at an interior point can now be written as follows:
Ki1 f + <71
2A t
+ u0i,j'
12A
23

i
i
' i
i, i
. <&+<&++<+i

24 A!'
g^rj1 + Â£ J j 2| 2 + 2 2
1 .1 "+!_ j. r*+1. ?+i. j_ n+i ,.ri+1. _L ."+i
7#??1 X l/??1 7)r"1'1. 7I?+1 J- 7/?T* 7J?"'-1. 4. 7)?T1
______t ^.J+2 Mt,J + l Vl-l,J + 2 + Vl,J + 2 , gyUt,J + l ' Ul,J Vl~l,J + l V1,J + 1
' 24Ay V i 2 2 2 2
t/}1 + un+17- + vpf1 u?t1 + u"}1 - + v?t\
__g7 *41 ' ij-i -iiJ ij 1 ij-i ' *4-2 ij-i \
I 2'; ' 2 2 ' 2 J
1 f + 16 Re{ i
12A2
+ :16<Â£1 301#1 + _ ?+j2
I I
12A y2
)
1
~2Vi ^KT + + &1j + S+1) -
W-2J-} + ^*"-24+2 + + "+lJ+l)]
. -pdj +!27i7ihl wfEtt +
24A
= 0,
(2.29)
-i'4u^- 4- v- 1
11 4 ' t,J
*4
2AI + ....... 12A1
jn+l , n+1 n+1 , n+1 n+1 , n+1 n+1 , n+1
'/ ^+2,j + Ut+2,j-l Vi+l,j + Vt+2,j + M.+lj-l V,tJ + Vi+l.j
j:1 ; 2 2 2 2
^1 1 i.n+l ^..n+l 1 ,,w+l ,.n+l 1 n+1 -.**+1 I ,..n+l
JI +Ui,j-1 ViJ + Vi-lJ Ui-l,j + Vil,j Vi2,j \
. 21 2 2 2 ;
f "1------2---------------T-------+ 27--------2---------------2
i^l1 + ytf-i. 1 + 1 &
_L 7in+1 7in+1 4- 7in+1 7j+1 4- ItV't1
+ ViJ-l ViJ-~ 1 + V.J-2 ^t.j-l + V,j-2 v
!i 2 2 2 '
&+i8a*j 3Qv?r+tana a*j
! 12A2
- ie 12Ay2 ]
H27P"/1 27
1 24 Ay
r^t1
't.j-2
= 0,
(2.30)
24
I

- 24Ax
~v?j+2 jb 27^ 27V?}1 + vtf_\
, 24Ay
+
= 0.
(2.31)
l
Boundary [Modification
I
At the periodic boundary, we write
*1.}
1J : - Vni-l,j,
Plj - P 1 m-l,j j
uo j - Uni2 ,j )
Vo J = Vni2,j j
Poj Pni-2,ji
UniJ = U2,j,
VniJ = V2J,
1 ru,j = P2J,
Uni+lJ
Vni+l'j 7 U3,j,
Pni+lj = ftj,
(2.32)
where ni is the number of grid points in the x-direction (for one wave length).
On the 'solid; wall boundary, we change the y-direction difference to second
order, and maintain fourth order in the x-direction. Again, special attention
25
! .

i
should be j>aid to the approximate value of v. Letting (see Figure 2.12)
j+i)
-W-i + + tf.j+i)).
we obtain an approximation error for v of 0(Ax4 + Ay2). Thus, for j = 3 and
nj 2, the
coefficients of the discrete equations in the y-direction become
I
Vi-2,j+l t. vi-l,j+l \ Vi,j+1 | vi+l,j+l \
Ii 1 f | . 1 1 |
,1 1' ^-2,i Vi-U vi,j 1 Vi+l,j
Fijgure 2.12. Approximation of v on solid wall boundary points.
Ann =
An |=
a. i|
ij=
A, j|-
i
Bnn |j=
Q,
_1___1_ Vn
Ay2 Re 2Ay2
11 v,
AtfRe + 2 Ay2
0,
3 1 1 5 2.0 27
2 At + Re ^ 2Ax2 + Ay2' + 48 A
0,
26
^ + 2Ay(v Vg)
I

Bn =
Bs =
Bss =
Re Ay2 2 Ay
1
+
Re Ay2 2 Ay
0,
0 . 3 1 1.5 2.0 , 27 . . 1 ,
2 At ReK2Ax2 Ay2 48Ax
Bss = 0>
2Ay
.4=
i
DP =
Â£)j\r = 0,
Ay
J_
Ay
Gnn
Gn
Gp
Gs
= o,
1
Ay
J_
Ay
= 0.
For j = 2 and nj 1, we must redefine some of these coefficients. For j = 2, we
make the changes
1
Re 3Ay2
Anew Aold ,
aN ~ AN +
Anew Aold
Ap ~ p ^ Re Ay2
1 2
and, for j = nj 1, we set
A new a old 1
~A* + Re 3Ay~2
Anew Aold , ^ ^
Ap -Ap +TeAy*'
27

, ! ' i
' I "
I'
! I
'' i I
I
CHAPTER 3
! ALGORITHM
i i
i
i
! Distributive Relaxation
'I :
|i ;
At eachltimje step, we must solve a spatial system, which we write as
! ApUij CpPij = (Su)p, | ] (3.1)
BpVij DPPi,j = (Sv)P, (3.2)
! -Fpmj GpVij = (Sm)P, (3.3)
for i = 2,3 1, j = 3,4,...,nj 1. We take the full fourth-order scheme
'! 1
as the example1, since the other cases can be obtained by setting some of the
' I I
I I
1 I I
coefficients' to zero. The right-hand side of the above system is written as
' I '
| ;
= (Aeeu{+2j + AEUi+i'j + Awiii-ij + Aww^i-i,j
: i i
j,j | +ApfpfUij+2 -f ApfUiij+1 -(- A\$Ui j_l + AssUi,j-2
i i 1
: j -\-CwwPi-2,j + CwPi-i,j + CeP{+ij) + {Su)p,
(Sv)p = (BEEVi+2,j + BEVi+l,j + B\yVi-l,j + BwwVi-2,j
\ ;
i i
j j +BffNUij+2 + BpiVi,J+i + BsV{,j-i -f BSsVi,j-2
I 1 +DssPi,j-2 + DsPi,j-\ + DffPi,j+i) -(- (5v)p,
L
(Sjn)p = (FEEui+2,j + FEui+u + FwUi-i,j
h. !
1 : !
, I
i I
i i 1

+Gjvj\fViiJ+2 + GrtUi'j+i + GsVij-x) + (Sm)p.
Now (3.1)|(3.3) can be expressed in the following matrix form:
Ap 0 Cp
0 Bp Dp
Fp GP 0
(3.4)
Up Su
vp - Sv
Pp Sm
(3.5)
It is obvious that (3.5) is not diagonally dominant, so the traditional point Gauss-
I
Seidel relaxation scheme may not converge. Also, although block Gau6s-Seidel
can be used to solve this kind of problem, it is fairly expensive.
' I
Some researchers (c.f., Chorin [13] and Patankar [14]) have developed methods
I
for solving system (3.5) that involve the usual type of relaxation on u and v,
, 'i
followed either by adding a ^ term to the continuity equation to convert the
j
third equation to the same form as the other two equations so that the traditional
Gauss-Seidel relaxation can be used, or by using a pressure correction iteration
to make the continuity equation satisfied. The problem with these methods is
i i
that, not only do they tend to converge very slowly, but they may not always be
consistent ,with ,the real physical problem in some cases.
Here, -vye introduce a so-called distributive relaxation scheme, which is similar
: i
. i
to but somewhat different from the distributive scheme developed by Brandt^12!.
Numerical experiments show that our distributive relaxation method, when cou-
pled with the multigrid technique, is much more efficient for solving incompress-
29
i

i
Gauss-Seidel relaxation for the first and second equations in (3.5), which we
write as follows:
Up
vP
(iSu)p + CpPp
Ap
(Sv)p + DpPp
Bp
(3.6)
(3.7)
We then adjust up,vp and Pp simultaneously to satisfy the continuity equation
I
and both momentum equations, using the corrections given in the form
ttp = Up Su,
U*Â£w = uE + 8u,
vÂ¥w = vP- 8v,
vw = vjf + 8v,
P2ew = PP + 6P. (3.8)
The corrections of other components in the system are set to zero.
'VN + 6
uE + (36
Figure 3.1. Distributive relaxation.
I
The defining equations for these corrections are
j
i {Ap -)- Ap)8u Cp6P 0,
: 30

(jE?at + Dp)Sv DpSP = 0,
(Fe + Fp)6u + (Gn + Gp)8v = (Sm)p, (3*9)
where
(\$m)p = (Feeuee + Feue Fpup + Fwuw
+C?jv jv "^jv jv + GrfVpf Gpvp + Gsvs) + (Sm)p.
This yields
8u Cp Bff + Bp A y
8v Dp Ap + Ap Ax
(Sm)P
(Fe + Fp)P + (Gn 4- Gp)'
= (ISv,
- (3-10)
Up
Multigrid Algorithm
For the large grids and long-time evolution models considered here, the usual
relaxation methods are much too inefficient by themselves. To obtain optimal
efficiency, we use a multigrid method based on the relaxation schemes developed
j
in the previous section. For simplicity of discussion, we consider only the two-
grid case as depicted in Figure 3.2. We use a full approximation scheme(FAS) to
accommodate nonlinearities.
I
A two-le;vel FAS algorithm for an equation of the form
P
6v
8 u
8P
LhV'h. = fh
31
(3.11)

may be described loosely as follows:
i) relax on LhUh = fh,
ii) solye L2hu2h = L2hllhuh + i^{fh Lhuh),
' i
iii) replace uh <- uh + I2h{u2h ~ (3.12)
i
1
The notation we have introduced includes the difference operators Lh and L2h,
the restriction operators lÂ£h (for the approximation) and 7^ (for the residual),
Figure 3.2. Two-level staggered grid.
The following subsections describe this algorithm in more detail.
Fine Grid'Relaxation and Residual Evaluation
We first perform distributive relaxation for (3.5), then calculate the residuals on
the fine grid level:
32
I

a; = (S)r (Apvij CrPb),
nS = - (Bpvfj DPpK),
St = (&*){, (FpuSj GfvSj). (3.13)
Coarse Grid Relaxation
Relaxation on the coarse grid requires restriction of the approximations and resid-
i
uals from the fine grid to the coarse grid to provide initial guesses and right-
hand-sides,, respectively. For restriction of approximations, we use the following
l
stencils:
I2Hh(u)
ll\v)
llh{P) :
1 l
2 2
(3.14)
For restriction of residuals, we use the following full-weighting stencils, which
come from: the so-called area law developed by Liu[15](Figure 3.3):

ill
8 4 8
111
8 4 8
33

jfW :
' I
I
I I
I

1 1
8 8
1 1
4 4
1 1
8 8 5
1 1
4 4
1 1
4 4
I
I
Figure 3.3. Restriction for Ru-
Tliis leads to the following coarse grid system:
i A2h ,2h /ry2h ri2h
i -AP utj CP Pij =
\ 1
I (-^Â£Â£t+2,j + -^Â£**+1,7 + A?wUi-1,5 +
I +-^5wJ+2 + + A2shsUitj_2
i i
! '
^4plU;ij + Cw\V Pi-2,i + CftPi-U Cp1 Pitj
^ +cfi+w) + i?()fl + (Su)f,
j ;
! n2hn2h n2h p2h
i
I {^EE^i+2,3 + BlsVi+lj + + BlyWVi_2,j
I
i ;
I 1 3A
(3.15)
(3.16)

(3.17)
1 Bph V{j + DfsP^2 + Dgh Pi,j-1 Dp1 Pij
+D2NhPitj+l) + Phh(v)Rv + (Sv)2ph,
-nh<5 O =
~{^EhE^i+2,j + P^^i+l ,3 + Fwi-l,j ^Pht,j
+ ^MVt,j+2 + OiJ+l + G2ShVi,j-1 G2pVij)
+Phh(m)Rm + (Sm)2p , (3.18)
where ii, i; and P denote restrictions of the respective u, u and P from the fine
grid approximations. We then perform one or more relaxations on (3.16)-(3.18),
using these restricted values as the initial guesses. (In the practical multilevel
I
scheme, we would recursively solve (3.16)-(3.18) by way of corrections from yet
j
coarser grids.)
Interpolation of Correction to Fine Grid
I
The next step is to interpolate the corrections from the coarse grid to the fine
I
grid: j
I
I
Uh - Uh + I^h{u){u2h Ilh(u)uh),
Vh r vh + I2h{v){v2h Ilh{v)vh),
PH + lZh(P)(P2h-Phh(P)Ph).
35
Ph
(3.19)

Here we use bilinear interpolation, given by the following stencils:

&() :
&(P)
3 1
4 4
_9_ _3_
16 16
iL JL
16 16
(3.20)
Figure 3.4 depicts the relations between u2h and uh, v2h and vh, P2h and Ph
corresponding to the above stencils.
O
u
2 h
o-
u
o
o-
v
2 h
V
o o o o
>2 h
O ph O
o o
Figure 3.4. Bilinear interpolation for u, v and P.
The FAS V-cycle based on the ingredients described here performed well in
all of the cases we have considered. Typical convergence factors per V(l,l) cycle
(c.f., Brandt [12]) were about 0.3 for the fully-implicit scheme.
36

CHAPTER 4
1
TEST CASES
In this: chapter, we first describe the initial conditions used for our experi-
' I
ments, then show well how our schemes work for the cases of low, medium and
high Reyn bids numbers.
Initial Small Disturbances
i
I 1 .
To study transition, we considered the least stable eigenmodes of the linear
! i
theory. For such modes, the linear theory provides an approximate analytic
1!
solution (which is accurate for small disturbances) to the governing equation
(2.9)-(2.11) given by
! i
u = eReal(u(y)eiax-iu,t),
1 I '
!
'! v = eReal{v(y)eiax~iu,t), (4.1)
|
where e is Ja small positive number and Real stands for the real part of. The
complex function = R + ii and the temporal frequency u; = wr + iujj come
from solution to the Orr-Sommerfeld eigenvalue problem. The Orr-Sommerfeld
I.
equation that governs the linear stability of parallel shear flow is
((^2 ~ 2)2 ^{(-o - ~ a) ~ Q^o}]^ = 0,

(x, -1) = (4.2)
The particular problems chosen here for study used Re = 2000, 7500 and 20000,
and a = 1;. Equation (4.2) was solved by a spectral method with 50th-order
Chebyshev polynomials. The only observed growing mode for each Reynolds
number is shown in Table 4.1. The eigenfunctions for those modes are shown in
Figure 4.1 (a), (b) and (c), respectively.
Re Wr <*>i
2000 0.31210050 -0.01979860
7500 0.24989154 0.00223497
20000 0.20966415 0.00330604
Table 4.1 Temporal Frequency for Different Re.
Note that the stream function is given by
i i> =
I
so, the streamwise and normal disturbance velocities can be obtained as
i u =
v = -i{y)Sax-WRt).
Setting a = 1, we obtain the eigenfunctions for both u and v in the following
r(y) = + =
38
form:

P(y) = r + <4) = -i'Kv)-
(4.3)
The initial jdisturbance is chosen from (4.1) at t = 0:
' i
u = e(^cosi ffisinx),
v = e(Vjsinx), (4-4)
where e is set to 0.0001. We compared the results obtained from our numerical
scheme against (4.1), which is the analytic solution for a small disturbance; as a
performance measure, we used the disturbance kinetic energy
E = i J ^iy J (u + v,)da
(4.5)
For sufficiently low amplitude waves, this should exhibit exponential behavior:
or
E(t) = E0e2wit,
= 2 w/Â£,
C'O
(4.6)
where Eq is the disturbance energy at t = 0.
39

eigenfunction of u
eigenfunction of v
width of channel
I
Figure 4.1.(a) Orr-Sommerfeld eigenfunction with /Ze=2000.
40

eigenfunction of u
I
Figure 4.1.(b) Orr-Sommerfeld eigenfunction with i2e=7500.
41

: I
eigenfunction of u
l

I
Computational Results
We chose Re, ujr and o/j to be the values shown in Table 4.1, which define
the least stable, modes for the corresponding Reynolds numbers. The parameters
used in our
experiments were
A< = ^i = -/500
500 u;,.
L = 27r,
nj = 34,
ni = 66 and 130.
Though other schemes, like semi-implicit second-order and semi-implicit partial
i i
i
fourth-order schemes, were also used, they did not yield satisfactory results.
For the first two schemes, second-order and y-direction fourth-order, we took
; i
Re = 7500 and tested for accuracy. All three Reynolds numbers are used for the
I '
fourth-order scheme. Distributive relaxation coupled with FAS multigrid proved
efficient for each of these schemes, with a convergence rate ranging between 0.3
J
and 0.4 for the Simple V(l,l) cycle (see Figure 4.2). Figures 4.3-4.10 display the
computational results.
i 1
43

LoglO(reaidm) LoglO(reaidv) LogiO(reaidu)
i
i
(c)
Figure 4.2. Convergence history for
(Re = 7500, mesh: 34 x 130)
(a) u, (b) v, (c) mass.
44

log(E(t)/E(0)) log(E(t)/E(0))

log(E(t)/E(0)) log(E(t)/E(0))

7T
T/To=0.5
T/To=0.75
Figure 4.6. Contours of disturbance stream function
(second-order fully-implicit, one period).
, 47

T/To=0
T/To=0.25
7T
T/To=0.5
T/To=0.75
0
7T
T/To=l
27T
Figure 4.7. Contours of disturbance stream function
(cross-channel fourth-order, one period).
48

T/To=0.25
7T
T/To=0.5
T/To=0.75
Figure 4.8. Contours of disturbance stream function
(fourth-order fully-implicit, one period).
49

I
Theoretical result, T/To=20, |^|MAX=2.46E-04
fourth-order in y-direction, T/T0=20, |^LaX=2.13E04
! I
0 ^ 7T 2n
1 ] second-order, T/T0=20, |^|MAX=1.42E-04
Figure 4.9. Comparison of contours of disturbance stream function
( Re=7500, time=20 periods).
(a) linear theoretical result, (b) fourth-order fully-implicit,
(c) cross-channel fourth-order, (d) second-order fully-implicit.
I 50

From the above results, we see clearly that the fourth-order fully-implicit
scheme is very accurate and efficient for simulating the process of transition, at
least at the small disturbance stage. We compare our results with the linear
theory in Tables 4.2 and 4.3.
. 1 theoretical fourth-order cross-channel second-order
, ' results fully-implicit fourth-order fully-implicit
average value 0.2235 xlO-3 0.2245 xlO-3 0.1931 xlO3 0.1223 xlO3
relative error. : i 0.0% 0.45% 13.60% 45.28%
Table 4.2. Comparison of the Increasing Rate of Disturbance Energy
: | (in 20 Periods).
I 51

1 1 I theoretical fourth-order cross-channel second-order
1 i i i ; results fully-implicit fourth-order fully-implicit
; 1 observed value 2.46 x 10"4 2.48 x 10~4 2.13 x 10-4 1.42 x 10-4
: | ; relative error 0.0% 0.81% 13.41% 42.27%
phase lagging 0 0.025tt 0.278tt 0.266tt
Table 4.3. Comparison of the Maximum Stream Function Values
j and Phase Lagging (after 20 Periods).
, |
l Concluding Remarks
!
i
.1 ;
From the above test cases, we conclude that
High-order finite differences can be very accurate for simulating transition.
Distributive relaxation coupled with multigrid can accelerate the iteration
process and make the more accurate and stable fully-implicit schemes com-
I
' i ;
parable in cost to the explicit schemes at each time step.
The simple two-dimensional problem studied here shows promise as a very
j I
efficient and feasible approach for more practical problems, including three
, i ,
dimensions and complicated regions, especially in the context of adaptive
refinement.
I I
I
52

CHAPTER 5
CONCLUSIONS AND OPEN QUESTIONS
Synopsis of the Thesis
I
The fundamental questions motivating this work are: Is the finite difference
method accurate enough to simulate transition? Can multigrid efficiently increase
the power of the finite-difference method? Can higher-order finite-difference
schemes obtain acceptable accuracy, enough to compare with spectral methods?
In the work of Liu, Liu and McCormick [16], the first two questions were
answered for the case of short-time evolution with Reynolds number Re = 7500.
But the results were not fully satisfying, especially for the problem of phase
lagging. The results presented in this thesis represent a significant advance over
I
their results.
In Chapter 2, we established three numerical schemes based on the staggered
grid. Care must be taken in developing high-order schemes because u,v and P
/ '
are not represented at the same points. This task was somewhat simplified by
the fact that we only considered the small disturbance problem.
Chapter 3 described the tools for solving the equations developed in Chapter 2,
including the so-called distributive relaxation method and the full approximation
I
scheme (FAS) of multigrid method. We then showed that these methods are

easily applied to simulate flow transition.
r
Chapter 4 documents a first step in using the developed method to simulate
I
i
the process of transition. The temporal problem was chosen for illustration.
Excellent performance was observed for the fourth-order fully-implicit scheme
I
over a range of Reynolds numbers.
I
We summarize our results as follows:
1. Staggered grids and fully-implicit time-marching schemes can obtain high
accuracy.
2. Multigrid methods substantially accelerate relaxation processes, making the
i
more accurate and stable fully-implicit and semi-implicit time-marching
schemes comparable in cost to explicit techniques at each time step.
3. For small disturbances, finite difference schemes with high order accuracy
achieve excellent agreement with the linear theory in both phase and am-
I
plitude. i
Open Questions
Though numerical simulation of transition has experienced tremendous progress
;l
over the past decade (c.f., Kleiser and Zang [17]), there are many remaining un-
solved problems. In the special direction of our research, a number of avenues for
further investigation appear to be open.

The schemes developed in Chapter 2 are very effective for the simple model
problems we have treated. They seem to promise a very efficient and feasible ap-
proach for more practical problems, including three dimensions and complicated
regions, especially in the context of adaptive refinement. However, the relaxation
method introduced in Chapter 3 can only be used for 2-dimensional problems, so
a new relaxation method must be developed for three dimensions. Like line re-
laxation for the 2-dimensional problem, perhaps a plane relaxation method may
be used to obtain better smoothing. For multigrid, although there seems to be no
apparent difficulty in extending it to more complicated problems, we must pay
special attention to the boundary points.
In Chapter,3, we only used the temporal problem in our tests. For more
practical cases, we should also consider the spatial problem, where the periodic
boundary condition in the x-direction would be changed to a Dirichlet condition
up-stream and a buffer domain for the non-reflecting outflow boundary (c.f.,
Danabasoglu, Biringen and Streett [18]). The main difficulty here will be that
the computational domain will be 20 or 30 times as large as the one we have
taken in this thesis, and local refinement should be used downstream to simulate
the details of transition.
These are but a few of the directions for further work that appear to be
promising. The study of multilevel adaptive methods for transition flow with
i
roughness elements, an area which is almostly unexplored, appears to be a rich
area for continuing research.
55

BIBLIOGRAPHY
[1] 0. Reynolds. An experimental investigation of the circumstances which
determine whether the motion of water should be direct or sinuous, and of
the law of resistance in channels. In Phil. Trans. Roy. Soc., A 174, page 935,
1883.
[2] Lord Rayleigh. On the stability, or instability, of certain fluid motions. In
Scientific Papers of Lord Rayleigh 1, pages 474-487, New York: Dover, 1964.
i
[3] W.M.F. Orr. The stability or instability of the steady motions of a perfect
liquid and; of a viscous liquid. Proc. Roy. Irish. Acad. 27, pages 9-69, 1907.
[4] A. Sommerfeld. Ein Beitrag zur hydro-dynamicschen Erklarung der turbu-
lenten Flu'ssigkeitsbewegung. In Atti. Congr. Int. Math., ^th Rome, pages
116-124, 1908.
[5] W. Tollmein. f/ber die Entstehung der Turbulenz. Nachr. Ges. lYiss.
Gottingen \Math.-Phys. Kl., pages 21-44, 1929.
[6] H. Schlichting. Zur Enstehung der Turbulenz bei der Platten-stromung. In
Nachr. Ges. Wiss. Gottingen Math.-Phys. Kl., page 182, 1933.
I

[7] B.J. Bayly and S.A. Orszag. Instability mechanisms in shear-flow transition.
In Ann. Rev. Fluid Mech. 20, pages 359-391, 1988.
[8] J.-P. Zahn, J. Toomre, E.A. Spiegel and D.O. Gough. Nonlinear cellular
motions in Poiseuille channel flow. J. Fluid Mech. 64, pages 319-345, 1974.
l
[9] S.A. Orszag and A.T. Patera. Subcritical transition to turbulence in plane
channel flows. Phys. Rev. Lett. 45, pages 989-993, 1980.
[10] T.A. Zang and M.Y. Hussaini. On spectral multigrid methods for the time-
dependent Navier-Stokes equations. Appl. Math. Comput. 19, pages 359-372,
1986.
i
|
[11] F.H. Harlow and J.E. Welsh. Numerical calculation of time-dependent vis-
cous incompressible flow with free surface. Phys. Fluids 8, pages 2182-2189,
1965. ;
[12] A. Brandt. Multigrid techniques: Guide with application to fluid dynamics.
GMD Studie, GMD, St. Augustine, 1984.
[13] A.J. Chorin. A numerical method for solving incompressible viscous flow
problems. J. Comput. Phys. 2, pages 12-26, 1967.
i
[14] 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, pages 1787-1806, 1972.
: 57

[15] C. Liu. Multilevel adaptive methods in computational fluid dynamics. Pli.D.
thesis, University of Colorado at Denver, i989.
i
[16] C.Liu, Z.Liu and S.McCormick Multigrid methods for flow transition in
planar channel. Comput. Phy. Comm. 05, pages 188-200, 1991.
[17] L. Kleiser and T.A. Zang. Numerical simulation of transition in wall-bounded
shear flows. Ann. Rev. Fluid Mech., Vol. 23, 1991.
i
j,
[18] G.Danabasoglu, S. Biringen and C.L. Streett. Numerical simulation of
spatially-evolving instability. ICASE Transition and Stability Workshop Pro-
ceedings, 1989.
i
58

Full Text

PAGE 1

I I ' : ' ' I AND FINITE DIFFERENCE METHODS I FOR FLOW TRANSITION IN A PLANAR CHANNEL by Zhining Liu University of Science and Technology of China, 1986 ; I 'i :I "' ', ; I I i ; I I I A thesis submitted to the of the requirements for the degree of Master of Science Applied Mathematics 1991

PAGE 2

' I This thesis for the Master of Science degree by Zhining Liu has been approved for the Department of Mathematics by ------Steve McCormick Thomas Russell

PAGE 3

' 'I 'I Zhining Liu (M.S., Applied Mathematics) 'I I I Multigrid Fipite Difference Methods for Flow Transition in a Planar Channel Thesis Professor Steve McCormick I! I The flow probleiD: has been central to the theory of fluid motion for over a century. Developing an understanding of this procedure will be very helpful 'I I to the research both theoretical and applied fluid dynamics. A computatiomtl study of th; stability of planar Poiseuille fiow in a channel is presented. Using time steps on a staggered grid, we develop a second-order scheme, a in the streamwise direction and fourth-order in the normal direction and a full fourth-order scheme. A very efficient time-marching procedure a multigrid solver is shown to produce very accurate results. The results show that the full fourth-order scheme gives very good ; I agreement both amplitude and phase with the linear theory obtained by solving I the equations. :I I The form and tent of this abstract are 'I I Steve McCormick 'I 'I I I

PAGE 4

I I i i I I I '' i I i :' I. '' '! To my parents

PAGE 5

:I-I ,, CONTENTS II List of Bigures . . . I I List of 'IDables l ............................ i I I : 1 INTRDUCTION I ; : 2 EQUATIONS FOR 2:jD PLANAR CHANNEL FLOW Equations and Boundary Conditions . . . . . I I Governiig in Perturbation Form . . . T I Fully-Implicit Finite Difference Scheme I 'I I Cross-Channel Fourth-Order Finite Differences . I .. Fully-Implicit Scheme . I Interior Points . . . . . I Modification ....... 3 ALGORITHM n t h lr ; R 1 r Is n live' e axa ton . . . . . . . . . . . . . Multigrid . . . . . . . . . . . . File Grid Relaxation and Residual Evaluation . . . . . I I :I I II I i I Vll ix X 1 4 4 5 7 14 17 17 25 28 28 31 32

PAGE 6

,i , '' , ; : I I I Grid Relaxa.tion . . . . . i l i of Correction to Fine Grid : ., .............. 33 .............. 35 I 4 TEST CASES 37 II : 'I Initial Small: Disturbances . . . . . . . . . . . . 37 C :1. :aiR 1 esu ts . . . . . . . . . . . . . 43 I' I I 1temarks . . . . . . . . . . . . . . 52 I '' I 5 CONCLUSIONS AND OPEN QUESTIONS 53 i .: ' ' Synopsis of the Thesis . . . . . . . . . . . . . . 53 II 0 I 0 ,, Open . . . . . . . . . . . . . . . 54 I I o BIBL]OG]ltAPHY i 0 I I I :.1 I 'I 0! I I , I I I oO . . . . . . . . . . . . . . 56 Vl

PAGE 7

I I LIST OF FIGURES 1.1 Reynolds' transition experiment.. . . . . . . . . . 2 I : 2.1 :Planar channel fl. ow. . . . . . . . . . . . . 5 2.2 2.3 2.4 2.5 2.6 2.7 2.8 I 'staggered grid. ............. i I Grid 'neighbors for x-momentum equation. I Neighbor points of up. I I I Neighbor points of vp. :I Orid neighbors for continuity equation. I wall boundary condition ..... Neighbor points of up for u-equation. 7 8 9 10 12 13 19 2.9 Fourth-order approximation for v at a u point (interior points). 20 :I 2.10 N;eighbor points of Vp for v-equa.tion. . . . . . . . . 21 2.11 points of Pp for continuity equation. 23 I 2.12 of v on solid wall boundary points .. 26 'I I 3.1 Distributive relaxation. 30 3.2 I Tlo-level staggered grid. 32 I 3.3 Restriction for Ru. .. 34 3.4 Bilinear interpolation for u, v and P. II 36 ; I

PAGE 8

' I 4.1. ( a:)l eigenfunction with Re =2000. I 4.1. (b) Orr-Sommerfeld eigenfunction with Re =7500. 4.1. 4.2 4.3 4.4 4.5 4.6 4.7 4.8 (c): Orr-Sommerfeld eigenfunction with Re =20000 ........ I co:nvergence history for (Re = 7500, mesh: 34 x 130 ) (a) u, (b) v, (c) mass.. . . .................. t b. ur ance energy mcrease (second-order fully-implicit, 3 periods) ............... i energy increase ( cJoss-channel fourth-order, 3 periods). I '' energy increase (f9urth-order fully-implicit, 3 periods). (a) Re.= 2000, (b) Re = 7500, (c) Re = 20000 .. '! Contours of disturba.nce stream function (second-order fully-implicit, one period). :I cJntours of disturbance stream function I (cross-channel fourth-order' one period). clntours of disturbance stream function 'I I (fourth-order fully-implicit, one period). I 4.9 C<;>mpa.rison of contours of disturbance stream function I (Re = 7500, time=20 periods) I ( al) linear theoretical result, (b) fourth-order fully-implicit, 40 41 42 44 45 45 46 47 48 49 ( t) cross-channel fourth-order, (d) second-order fully-implicit. 50 4.10 Comparison of different schemes with disturbance energy. I 51 I I ... Vlll I

PAGE 9

, I 'i : i LIST OF TABLES 4.1 Frequency for Different Re. . . . . . . 38 'I I I 4.2 Compa;rison of the Increasing Rate of Disturbance Energy Periods). . . . . . . . . . . . . 51 4.3 I Cdmparison of the Maximum Stream Function Values PKase Lagging (after 20 Periods) .......... i. I I ,, I i. i 'I II ['' I I. 'I : I I I i 52

PAGE 10

ACKNOWLEDGEMENTS :I I would 1like. to take this opportunity to thank my advisor, Professor Steve McCormick; for the guidance and encouragement that he provided during a.ll I I phases of thf of this thesis. It has been a privilege to work with him, I and I hopeJhaf I have acquired some of his insight in the process. Thanks go to 'I. I Professor Ghaoqun Liu, who provided me with close supervision and who guided I I I me step step through algorithm development. Thanks also go to the other two : i I members q my: committee, John Ruge and Tom Russell, for their contribution I 1: to my development. 'I : i I wouldl also like to thank my parents and my friends for their unflagging i' support. Wheir encouragement kept me going, and I am in their debt. It is I : impossible1 ro list them all, but a few deser:ve special mention: Debbie Beltz, :I Gaoming Yang and Jonathan Hazen. I I am to NASA for their support of this work under Research II 'I Contract 'I ,I I I I X I I i

PAGE 11

I I CHAPTER 1 INTRODUCTION The process from laminar to turbulent flow in a sheax layer has been I the subject of intense study for over a century, but it still remains a challenging and impor\rt for reseaxch. The process of transition is very complicated. Early work 'was initiated by Reynolds [1 J, Rayleigh [2) and others. Orr [3) and Sommerfeld! [4) later reformulated Rayleigh's theory for viscous incompressible lluids. [5) and Schlichting [6) successfully predicted the structure of Orr-Sommerfeld modes in noninflectional flows. I has yet to: an understanding of the first 'i However, current technology classic transition experiment completed by Reynolds more than one hundred years ago. He showed that a thin I filament of: fYe in water flowing along a transparent glass tube would remain well defined and/ stra.ight if the Reynolds number of the flow, Rd, was less than about 2,500; but for higher Reynolds numbers, the filament rapidly became convoluted, thiclrened up, the dye dispersing throughout the tube (see Figure 1 (a) I and (b) ). Here, Rd = udjv, where dis the diameter of the tube and u is the I centerline of the mean flow. But by using spark photography to ca.pture I details of the flow, we can see that in the fully mixed region of the tube there are many 1,ortex-shaped streamlines. The vortex moves 3-dimensionally in the :I

PAGE 12

.;rater tank transparent glass tube (a) (b) (c) Figure 1.1: Reynolds' transition experiment. spatial direftio!). with a frequency in the several thousands (see Figure 1 (c)). I I As currJntly understood, the transition to turbulence is, in fact, a very comi plicated 3-::0, nonlinear and multistage process, which begins with a sequence I of instabilities on a succession of more and more complicated basic flows ( c.f., I I Bayly and ;Orszag [7]). Loosely speaking, the fluid flow will be laminar when the '. Reynolds is small, and become turbulent when it is high, which that as 1Reynolds number increases to some value the laminar flow will begin I to become: People have established a so-called small disturbance theory to what happens with respect to time when a small disturbance is '' added to the steady laminar flow. This theory can tell us what kind of profile is unstable, which mode of the oscillation will be amplified the fastest, and how to change the flow parameters to delay transition. Unfortunately, it 2 :I

PAGE 13

. I 0 I I '' tells us :nothing about the nonlinear instaobility and the vortex breaking 0 0 stage, are both very important to the study of transition. I I Numerical simulation of flow transition in a channel has been successfully done I I I 0 by several ( c.f., Zahn, Toomre, Spiegel and Gough (8] and Orszag (9]), I most of spectral methods with explicit time-marching and non sta.ggered grids (c.., Zang and Hussaini (10]). The difficulties with numerical I I' simulation for transition result from the long-time evolution and the sensitivities to small requiring high accuracy and efficiency from the numerical meth1 ods. Althc;>hgh )t is recognized the staggered grid ( c.f., Ha.rlow and Welsh (11]) and time-marching approaches ha.ve advantages, there has really been to apply them, especially in the field of transition. The pjposl of this thesis is to develop an accurate finite difference scheme I on a .grid with implicit time-marching. Since we must solve a spatial 'I I I system at e1ach time step, a full approximation multigrid algorithm (Brandt (12]) '' is applied tf accelerate the convergence process. I The test problem is a 2-D planar Poiseuille channel flow. Different Reynolds I ; i . numbers were chosen to test the accuracy of our schemes, and good results were 'I :I obtained our final fourth-order fully-implicit scheme. :I II I o I o 1 [ I I 3

PAGE 14

CHAPTER 2 I GOVERNING EQUATIONS FOR 2-D PLANAR CHANNEL FLOW I In this chapter, the governing equations for pla.nar Poiseuille flow in analytic I and discretized 'forms are presented. Several types of schemes are developed for I the temporal stability problem. I Governing Equations and Boundary Conditions I I A two-dimensional, time-dependent, incompressible Na.vier-Stokes equations can be useq as the governing equations for planar channel flow: i 'I : i + 8uu + 8uv _.!:_( 82u + 82u) + 8P = O, 8t 8:v 8y Re 8:v2 8y2 8:v Bt' 8uv 8vv 1 ( 82v 82v) 8P -+-+----++-=0, (Jt 8:v 8y Re 8:v2 8y2 8y 8u 8v -+-=0, 8:v 8y (2.1) (2.2) (2.3) where U V are velocity components in the X-and ydirections, respectively, I Pis the pressure a.nd Re is the Reynolds number based on the centerline velocity, '' '' U, of flow, the channel half-width, h, and the viscosity parameter, v: I l j 'I Uh Re=-. v The conditions on the solid wall are I u(x, -1, t) = v(x, -1, t) = 0, 'I 'I (2.4)

PAGE 15

'' u(x,l,t) = v(x,l,t) = 0. (2.5) y I 1 i! :! ::::--.. '' 'i I Q_ I l :I .. 1 I / '' I' /' X 1 :' : --' i I Figure 2:1. Planar channel flow . I I we :Choose to use a staggered grid, boundary conditions for P are not : I j, needed on jt!he sblid wall. For simplicity, weo,nly consider the temporal insta.bility :1 : for channel I flow, and assume periodic boundary conditions in the x-direction: I I I i I u(x, y, t) = u(x + L, y, t), v(x,y,t) = v(x + L,y,t), 2 P(x,y,t) = P(x + L,y,t) + ReL, where Lis ,the wavelength of a specified mode of disturbance. : i I :' I! i i Governing Equations in Perturbation Form 0 'I 'I (2.6) flow is determined as the solution functions uo, Vo and Po of (2.1 )(2.3), [boundary conditions (2.5)-(2.6), tha.t satisfy the initial conditions: I :I 'I II '! i I 'I i II 0' v0(x, y, 0) = 0, 5

PAGE 16

i 2 P0(x,y,O) =-Rex+ C. (2.7) Here, C constant. This will be called an undisturbed solution. Since u0 v0 and Po are in:dependent of time, we have uo(x,y,t) = u0(x,y,O) = 1-y2 (2.8) and for v0 and P0 If we imp()se some disturbance at t = 0 on Poiseuille flow, it will generally I I i I I become time-dependent. For such a disturbed flow, in order to reduce roundoff error in the computational solution, we write the primitive variables in the I following :perturbation form: I I I I I I u = u0 + u, I v = v0 + v, P =Po+ P1 (2.9) SubstitutiJg (21.9) into (2.1), (2.2) and (2.3), a.nd using the fact that v0 = 0, we I may then the governing equations in the following perturbation form: (2.10) (2.11) (2.12) I: Here, for siwplicity, we drop the prime notation, so that u, v and P a.ctua.lly now stand for variables u1 v1 and P1 The boundary conditions ma.y I 6 I

PAGE 17

I :I : be rewrit tell as I I I I .. I I u(:z:, -1, t) = v(:z:, -1, t) = 0, u(:z:, 1, t) = v(x, 1, t) = 0, u(x, y, t) = u(:z: + L, y, t), v(x,y,t) = v(x + L,y,t), P(:z:,y,t) = P(:z: + L,y,t). Fully-Implicit Finite Difference Scheme 1. II (2.13) ': I We use+ uniform staggered grid for our problem (see Figure 2.2)'. A secondorder backJ..ard: Euler finite difference in time and a second-order central differ1 I : ence in space aie applied first to discretize equations (2.10) (2.12). 'I I I I 'I I : i : Figure 2.2. Staggered grid. First consider the x-momentum equa.tion (2.10). The related spatial com-, I : : I ponents 1!ui,i iare shown in Figure 2.3. With superscript n used to denote u i 7 'I i I :I

PAGE 18

' I at the nth' increment, 1Lo the undisturbed solution of the channel flow in I v they-velocity component defined at au point, then we discretize I I equation on the above staggered grids using a fully-implicit finite-difference scheme: u.'!-++1 1 + u'!-t1 u'!-+1. + u'!-t1 u':l+l. + u'!-t1 ,J I,J -&-l,J &,J &-l,J )/A:v 2 2 2 + 'J ,_ ,J = 0 A:v u '+1 '' o, 0 t,] 0 0 '; I I' I I I Vi,j+l Ui-l,j ui+l,j 0 0 0 pi-l,j p,.. t,] Vi-l,j v. t,] 0 Uij-1 0 0 I 2.3. Grid neighbors for x-momentum equation. 'i ; I (2.14) For simppcity, we may rewrite (2.14) in the general form (see Figure 2.4) ; I : AEuE Awuw + ANuN + Asus-ApupCpPp + CwPw = Su, (2.15) 8 i : i i I

PAGE 19

where 1Define I I '' 'i I' '' I I '! 'I I I I i I I I : I i :[' I I I I ; I : I I I UE Ui+t,j, uw Ui-1,j! UN Ui,j+l! us Ui,j-t, Up Ui,j, Pp Pi,;, Pw pi-l,j UN 1uw Up UE -0 1-0 11-Pvv Pp Figure 2.4. Neighbor points of up. 0 5 ( n+ 1 + n+ 1) U+l U 1 ,, .,, 0 5( n+l + n+l) U 1 U 1 ,, .,, 0 5( n+l + n+l ) Vn Vi-l,j+l vi,j+l 9

PAGE 20

'I I I I This I 'i Aw;= AP' I I Cp != 0 5( n+l + n+l) v. -. vi-l,i vi,i 1 1 : Ref!1x2 2!1x ( Uoi,i + Ue ), 1 1 Ref!1x2 + 2!1x (uoi,j + uw), 1 1 ------v 'Ret1y2 2!1y n' 1 1 I + --v., Ret1y2 2!1y 31 1 1 2 1 1 :2 tit + (u.,uw) + 2!1y (vnv.) + Re ( + !1y2 ), 1 1Cw = -, tix ; I 11 ( n+l n+l n+1 n+l )( ) Up-1 -4tt.p Su = -v. + v. '+1 + v. 1 + v. 1 .+1 -2y; + ; I" '4 t,J ,,, t-,, t-,, 2!1t I :I I ; I I I I I I t t j_ I I I 0 PN 0 0 t tVN I I I 0 Pp o 0 tVw tVP tVE I I I 0 Ps 0 0 tVs t I .. .. I I Figure 2.5. Neighbor points of Vp. write the difference equation for (2.11) in the general form (Figure 2.5) I I : i I (2.16)

PAGE 21

i I I I where ; i With I I I I : I : i 'i :I I I I I I I I : I :,I I I II : I I I I : i 'I I' I u., Uw -Vn -.v. -VE vi+l,j, vw Vi-l,jl VN Vi,j+b vs Vi,j-b Vp Vi,j, Pp Pi,;, Ps Pi,j-1 0 5( n+l + n+l ) ui+1,j ui+l,j-1 0 5( n+ 1 + n+l ) ui,i ui,i-1 0 5( n+l + n+l) vi,j+t vi,; ; 0 5( n+l + n+l) vi,j-1 v,,; I I i we get 1 BE T 1 1 Refl.z2 2/l.z ( Uoi,j + u.,), I i Bw :.: 1 1 I Refl.z2 + 2/l.x (uoi,j + uw), I 1 1 1 BN Refl.y2 -2/l.-y Vn, Bs l 1 1 :t Refl.y2 + 2/l.y v., I i 31 1 1 2 1 1 :T 2 ll.t + .2/l.z ( u., -Uw) + 2Ay ( Vn -v.) + Re ( fl.z2 + fl.y2 ), I 11 I I 'I

PAGE 22

I 'l .I 1 Dp1 = Dw = -, : !1y I Sv I= I .\ v_p-1 -4vp 2!1t tVN I _Y,p OPp Y,E tVP I Figure 2.6. Grid neighbors for continuity equation . The co*tinuity equation may be discretized as (Figure 2.6) :I where I I 'I I I :I I I I 1 FE=Fp=-, !1x 1 GN = Gs = -, !1y Sm=O. (2.17) The boundary is modified as follows. At the periodic boundary, we write :I I

PAGE 23

(2.18) where ni is the number of grid points in the x-direction (for one wave length). On the wall boundary, we use a. modification of the x-momentum equal tion (Figure 2.7). For example, for j = 2 we use the approximations 82u::::::: _4_('1LN-Up_ 0 = _ 4_Up, 8y2 3L\y L\y lAy ) 3 Ay2 Ay2 (2.19) 8uv ::::::: _1_( + t1i-l,j+l 'ILN +Up O). 8y L\y 2 2 (2.20) We simply. the appropriate coefficients in (2.14) at the boundary points. I For for j = 2 we make the changes Anew= Aold + _1 ___ 1_1 N N Re3L\y2 0 i Anew = Aold + P P Re Ay2 (2.21) UN fUp '-. ... -.... u=O Figure 2. 7. Solid wall boundary condition. For j = nj .,..1, where nj is the number of grid points in the y-direction, we use i I A new A old + _!_ _1_ s s Re3L\y2' Anew = Aold + _!__2 __ P P Re Ay2 13 (2.22)

PAGE 24

Cross-Channel Fourth-Order Finite Differences 'I Since tends to change more rapidly in the y-direction, either re I I fined grids[ higher order differences in the y-direction can be used to balance overall a.ccJracy. One of the methods we study in this paper uses fourth-order central difflences in the y-direction for interior points only, while maintaining second-order central differences for boundary points and for interior points in :I I the We use a fully-implicit scheme, since it is generally more accu1 : I rate than semHmplicit or explicit methods. For simplicity, we continue to use second-ordl dikerences for the continuity equation and the pressure terms in .I I I both equations .. The general form for the fully-implicit scheme in i I this case be written as follows: I I I '' -Apup + CwPw-CpPp = Su, iBEvE + Bwvw + BNNVNN + BNvN + Bsvs + BssVss I :I -Bpvp + DsPs-DpPp = Sv, Here, similJr the previous section, '; I :I I I I I I I I 14 (2.23) (2.24) (2.25)

PAGE 25

I I I I i i I I i I Vss V&,j-2 In (2.23), if we define then Aw As Ass '! I I '' I : I i : !] I I I .': i. I i j 'i I Ue 0 5( n+1 + n+1 ) U U + 1 .,, ,J Uw 0 5( n+ 1 + n+ 1 ) u,,; ui-1,; Un 0 5( n+l + n+1 ) ui,j ui,i+ 1 u. 0 5( n+1 + n+1 ) u.. u . 1, '' ''Vn 0 5( n+1 + n+l ) vi,i+1 vi-t,i+1 0 5( n+1 + n+l ) Vnn Vi,j+2 Vi-l,j+2 v. 0 5( n+ 1 + n+ 1 ) vi,j vi-1,j 0 5( n+ 1 + n+ 1 ) vi,j-1 vi-t,j-1 ii .J. 0 25( n+1 + n+1 + n+1 + n+l. ) vi,; vi,j+1 vi-1,j vi-1,j+t 1 1 ilef1x2 2t1x (uoi,j + ue), 1 1 lf-et1x2 + 2t1x ( Uoi,j + Uw), 1 1 1 1 _, Re 12t1y2 + 48 tiy Vnn, =l 14 11 ---+ ---(v -27v ) i fle 3y2 48 tiy nn n + 1e + :8 (27v.-v .. ), I 1 1 1 1 1 i ...., Re 12t1y2 48 tiy v.., I I 15 I I I

PAGE 26

3 1 27 1 2 5 -' + (u.,uw) + 48Ay (vnv.) + Re ( Ax2 + 2Ay2 ), 1 Cw:. = ; Cp Su1 I I ' + 1 ,,, ,,, 2y ;fii,i. i define 'I 0 5( n+1 + n+1 ) v., vi,j vi+1 ,; 0 5( n+l + n+t ) vi,i vi-t,; 0 5( n+l+ n+l ) V V "+t 7 .,, t,J v, = 0 5( n+l + n+l ) vi,i vi,j-1 Vu 0 5( n+l + n+l ) V. 2 V. 1 ':-t,J-u., 0 5( n+l + n+l ) ui+1,j ui+l,j-1 0 5( n.+l + n+1 ) Uw ui,j ui,j-1 .i Then, for.equation (2.24), we have 1 1 1 BE (=.I 2Lix (uo,,; + u.,), 1 I 1 1 Bw ':R 6. 2 + 2A (uoi,j +uw), e x ux 1 1 1 1 !-Re + 48 v., .. I I I 1 1 4 1 1 BN. l = 'Re 36.y2 + 48 Ay ( Vnn -27v .. ), I ; 1 4 1 1 Bs, := :Re 36.y2 + 48 6.y (27v,-v,.), 1 1 1 1 Bss. := --. ---.-v I Re 126.y2 48 Ay ""' I .I 16

PAGE 27

' '; I I 'I 'i Bp ,1= 3 1 27 1 2 5 + 2Jj.:z: ( Vn v.) + 48Jj.y ( Vn -v.) + Re (Jj.x2 + 2Ay2 ), 'I I Ds .i= I I Sv '' I I ; 1 pp =A I L.J.:J: 1 n n-1 -4v .. + v .. '' '' There is chahge to the continuity equation, so I i i I I 1 i i FE -Fp=-, I fj,:z: I GN 1 I Gp=-, ,, I Jj.y I ' Sm 0. I I '' Since the scheme is maintained at the solid boundary points, the boundary are the same as the previous section. I II. Fourth-Order Fully-Implicit Scheme I l I Since this is quite complicated, we treat the interior points and boundary I points I I I I '' i I 'I I Interior I I ;. I i Though ithe above schemes are fairly accurate, especially the cross-channel I : I fourth-ordelschfme, they are inadequate for simulating transition because of the I sensitivity su1=h flows to small changes and the need for long-time evolution. I 1 For these we developed a so-called fourth-order fully-implicit scheme on I I . ,, i I a with the aim of resolving the problems of lagging of the phase 17

PAGE 28

I i i and the rate of disturbance energy increase that we observed in the i above schenites. :since we only consider the temporal problem at this time, which ,, I I means that! we;only consider the instability with respect to time, there is no 'I difficulty periodic boundary conditions in the x-direction. ; I I The form of the fully-implicit finite-difference scheme can I ,I be written as: I I Here, I I I i "I I I :I I 1 1 I 'i ; I 'I :I I I I 'I i I I' : i : I I I 'I : i I I II 1: ; I I' 'I ,, I I I !" : i I I I I (2.26) (2.27) (2.28) UEE ui+2,il uww -Ui-2,j, VEE Vi+2,j, vww Vi-2,j, Pww pi-2,j, 18

PAGE 29

,I I Pss P;,,;-2 The other are defined in the previous section. Figure 2.8 depicts the neighbor points:lof for xmomentum equation. We determine the coefficients as follows. : I I nn v :I I 0 0 0 0 ,; ,. Vn .. ; vw Uww u v Uw Up Ue Uee u ---_,g.. ... _,g.. _,g.. ... _,g.. 1-u E I Pivw Pw Pp PE I Vs I .. I I' I d 0 0 0 I Vss I I I . . :.I rgure 2.8. Neighbor pomts of Up for u-equation 'I I I I that u and v a.re several orders of magnitude smaller than u0 we still use averages to approximate uu, uv and vv, i.e., we still use :t I + . + 1-l,J I-2,J 1_-l,J I-2,J Uww Uww = .. 2. .. -2 I I and so fori'h. However, since the terms u0 8 8u and are several orders larger :I I z y :.1 : . than the other (here, v is 'the y-velocity component defined at a u point), we must ai.lrolmate them with greater accuracy. Standard fourth-order central differences:Layibe used to obtain fourth-order accuracy, as long as we use a high : i 19 I I

PAGE 30

'l I order appr9xiniation for v at the midpoints of vertical cell edges. The I I order apptoxirdation we use here involves the points shown in Figure 2.9. The I formula is I I I I : 'i t I I i I : i I ,J I i l I i I I 'Pi-2,j-1 ', I I 1:, I : T. '' I 1 (9( n+1 + n+l + n+1 + n+l ) V V '+1 V 1 V 1 '+1 -32 '' '' ,_ ,, -,J ( n+ 1 + n+ 1 + n+ 1 + n+ 1 ) ) V 2 1 V 2 '+2 V+1 1 V+1 '+2 .J.J ,J-' ,, Vi-l,j+l t I t I Vi-l,j v. l,] Vi,j+l J I t I v. l,] t I Vi+l,j+2 Vi+l,j-1 t I I Figure Fpurth-order approximation for v at au point (interior points). I I. I I I to the above approximation, we obtain the following coefficients I AEE ; 1 1 + ( Uee + 4uoi,j ), ': AE -,-' 4 1 + ( Uee -21ue -32uoi,j ), ,, '' Aw : .,.. I 'I I Aww ,-;-I I I. ; 4 1 + Uww + 32Uoi,j)! ; 1 1 i (uww + 4uoi,;), 20 ' I '' I I

PAGE 31

ANN AN As Ass Ap CE Cw Cww Su I I I T I* I i ..!... T I ...:... T 1 : i I :-I I I u.. I :I I I I I I I .I ' I I I 1 1 1 1 -;' Re 12.6.y 2 + 48 Lly Vnn, ;1 4 1 1 Re 3.6.y2 + 48 Lly {vnn-27vn), .1 4 1 1 Re 3Lly 2 + 48 Lly (27v" v.,.,), 1 1 1 1 -Re 12Lly 2 48 Lly v.,.,, 3 27 27 5 1 1 2Llt + {u., -uw) + 48Liy (vnv.,) + 2Re { + Lly2 ), 1 27 Cp = 1 + ?L':':1 .,, . 2 26t YjVi,j NN ..... ---.. -r .. I PN fVnn 0 fVN ...... T 0 Pp fVn Uww Uw !VP r 0 Ps fVs fVS I 0 PssfVss 1_Vss I 0 0 Ue 0 0 Figure 2.10. Neighbor points of Vp for v-equation. 21

PAGE 32

Similarly, for the v-equation (2.27) (see Figure 2.10), let This yields : BEE -;-BE ,_!_ Bw ._;_ ' Bww '' ,_ BNN -BN l'i ...,... I Bs _:_ I' I Bss I -Bp -DN -Ds Dss I I I Sv 0 5( n+l + n+l ) Vww Vi-l,j vi-2,j 0 5( n+l + n+l ) Uee Ui+2,j Ui+2,j-1 __:_ + (uee + 4uoij), 4 1 + (uee-27ue-32uoi,;), 4 1 + (27uw-Uww + 32uoi,j), 1 1 (uww + 4uoi,;), : 1 1 1 1 :Re + 48 Vnn, 1 4 1 1 -' --+--(v -27v) 48 nn n 1 4 1 1 Re + 48 {27v.v .. ), : 1 1 1 1 --;.----v .. Re 48 27 27 5 1 1 + (vnv.) + 48Lly (vnv.) + 2Re ( + ), 1 24Lly' 27 = 1 + v!l:-1 I,J ,J 22

PAGE 33

. I ' 'I I I I I 0 Up 0 f-0 t v NN I 0 0 tVN I Pp UE UE -0 tVP I I 0 0 tVs "I .. figure 2.11. Neighbor points of Pp for continuity equation. I I The continuity equation, (2.28), becomes fourth-order with the following coeffi-: I dents (see Figure 2.11): I I I FEE ol FE :I Fw :I GNN I 0 : GN rl Gs : I Sm I -= 1 24Ll:z:, 27 Fp = 24Ll:z:' 1 24Ll:z:, 1 24.6y, 27 Gp = 24.6y' 1 24Lly, 0. The compllte discrete system at an interior point can now be written as follows: I 3 n+ 1 o I, 4 n + n-1 n+ 1 + 8 n+ 1 8 n+ 1 + n+ 1 ui,i : ui,j ui,j -7J.i+2,j ui+1,j -ui-t,j ui-2,j 0 2Llt + Uoi,i ... ---. -y2:Li:z: I 23 :I

PAGE 34

I ,, 1 ' + +---'(t+2,J 't+1,J 24Ll:z:l i 2 2 I + + + + -27 t,JI : t-1,J '' t-l,J + t-1,J t-2,J t-l,J I-2,J) ': 21 2 2 2 I 1 I + v!'+lt'+2 + v!'t+t2 + v!'+lt'+t + v!'t+tl +---' .( ,t,J t1J 1-1J . 1 1J + 27 t,J I,J 1tJ t,J 246-yl : 2 2 2 2 + n+l + n+1 nf1 + n+l n+1 + n+1 _27 ''! : a,J-1 vi-t,i vi,j + ... ui,j-2 vi-1,j-1 1'i.i-l) I 2: 2 2 2 I 1 I n+ + 16 n+1 30 n+l + 16 n+1 n+l --( -;ui+2,j ui+t,i ui-2,j Re I +: + + 11J. 11J 1 1J 11J-11J) :: 126.y 2 [ (: n+1 n+1 n+1 n+1 ) -2y3 -,9 V + V '+1 + V 1 + V 1 '+1 -:32 : '' 1 '1 I'1 1 I I I ( n+l '; + :n+1 + n+1 + n+1 )] vi-2,i-1 'l?i-2,i+2 vi+1,j-1 vi+1,i+2 +i27 P!'7-1 -27 p!'+t. + + 1 11J .t-11 ) t-2,J Q I = l I 24Ll:z: (2.29) : i :I 3 n+l '14 + n-1 n+l + 8 n+l 8 n+l + n+1 V V V -V+2 V+ 1 V 1 V 2 loJ I . I,J t,J + tL I ,J I ,J 1,J 1-,J '26-t OI,J I i 1 +1 1 +1 +1 +1 +1 +1 1 : i + _1 v"+ + v" u" + u" + +---. (1 1 1J I 1J i+1,j o+2,j + 27 i+11j ofl,j-1 1 1J I 1J 246-:z:; I 2 .. :t' .. 2 2 + u.' n+l + n+l n+l + n+1 n+1 + n+l _27 '', a,J-1 vi,j vi-1,j + ui-t,; ui-l,j-1 V&-1,j vi-2,;) 21 2 2. 2 1 I n+l + n+l n+1 + n+1 n+l + n+1 n+1 + n+1 ---, [ Vi,i+2 vi,j+t vi,j+2 vi,i+1 27 V&,j+1 vi,j vi,j+1 vi,j + I 2 2 + 2 2 v!'t: 1 + v!'t1 v!'t1 + v!'t1 v'!"t 1 + v!'t1 v!'t1 + v!'t1 -27 t,J i : I,J-l t,J 11)-1 + 11]-1 11J-2 1 1J-l t1J-2) f 2; 2 2 2 1 + 16v!'++ 1 1 30v!'T1 + 16v!'+1 1 v!'+2 1 --( ,: I t' I 1J . 1-1J 1-1J Re 126.:z:2 + [6v!'T+1 1 30v'!"t1 + l61/'T1 1 -v!'t1 2 + 1,,,, '' t,J-loJ.) : 126.y2 I I I + -Pi7/-_;.,\ + :27 -27 Pi7f-\ + ui,j!2 = 0, '! 246-y (2.30) I I 1. 24 .I

PAGE 35

427u':'t1 + u':'-+1 1 t ,J t ,J t,J t-,J + ' n+l _l 2..,. n+l 27 n+l + n+l -vi,j+2 T' '1vi,j+lvi,j vi,j-1 _.....::;_:......._; __ = 0. I I : Boundary :Modification i At the periodic boundary, we write I 'i I; I ul,j vl,j pl. ,J Uo,j Vo,; Po ,J 'Uni,j Vni,j Pni,j Uni+l,j Vni+l,j Pni+l,j 'Uni-lJ, Vni-l,j, Pni-l,j, '1Lni-2,j, Vni-2,j, Pni-2J, u2,j, v2,j, P2,;, u3,j, v3,j, P3,;, (2.31) (2.32) where ni is ihe number of grid points in the x-direction (for one wave length). On the solid: wall bounda.ry, we change the y-directlon difference to second order, and mai11tain fourth order in the x-direction. Again, special attention 25 I.

PAGE 36

should be Jaid. to the approximate value of v. Letting (see Figure 2.12) I I I J j v. .,, 1 (9( n+1 + n+1 + n+l + n+1 ) V V '+1 V 1 V 1 '+1 32 lJ -J -J i I we obtain approximation error for v of 0( + .6.y2). Thus, for j = 3 and nj-2, of the discrete equations in they-direction become 'I i I 'I I I :I Vi-2 j+1 Vi-1 j+1 v '+1 Vi+1,j+1 ( t I t'J t .I :I .II I I I I v . z,] j _t 1 t ,I I' I I I Vi-1,j v. z,] Vi+1,j :I I I :I FiJure 12.12. Approximation of v on solid wall boundary points. ll I :L ANN 0, AN 'L I i As iL r ,I Ass IL r I Ap :L :r I BNN 1! i= II 'I I 1 1 Vn .6,y2 Re 2.6.y2 1 1 -:--+lly2 Re 2lly2 0, 0, 26

PAGE 37

BN -1 Vn. Retl.y2 2tl.y' Bs 1 v, Retl.y2 + 2tl.y' '' Bss 0, Bp --;Dss 0, Ds ...!.._ .1 T tl.y Dp 1 -' D.y DN 0, GNN '-;-0, GN 1 T D.y Gp 1 -' D.y Gs j_ 0. T For j = 2 a.hd nj-1, we must redefine some of these coefficients. For j = 2, we j make the cna.nges Anew_ Aold + N N Re 3D.y2 Anew Aold + _.!.__2_ P P Re tl.y2 and, for j d: nj-1, we set i 'I A new A old + _...!:._ -2:__. 5 5 Re36y2 Anew Aold + 1 2 p p ---. Re tl.y2 27

PAGE 38

; I f, 11 'I I I, I I I I i : I CHAPTER 3 ALGORITHM Distributive Relaxation At step, we must solve a spatial system, which we write as "I i I 'i .i I I -Apu CpP. = (Su)p t,J t,J -Bpv DpP. = (Sv)p 11J 11J -Fpu -Gpv = (Sm)p (3.1) (3.2) (3.3) for i = 2, 3, ... r,.i -1, j = 3, 4, ... nj -1. We take the full fourth-order scheme I IJ I as the since the other cases can be obtained by setting some of the 1) I IJ coefficients' .to zero. The right-hand side of the above system is written as .! I 'i I : i : i I! I I,

PAGE 39

I I I I 'I Now (3.1):(3.3) can be expressed in the following matrix form: -Ap 0 -Cp Up Su 0 -Bp -Dp Vp Sv -Fp -Gp 0 Pp Sm I (3.5) It is (3.5) is not diagonally dominant, so the traditional point Gauss Seidel scheme may not converge. Also, although block Gauss-Seidel can be used to solve this kind of problem, it is fairly expensive. . I :I : Some ( c.f., Chorin [13] and Patankar [14]) have developed methods I for solving (3.5) that involve the usual type of relaxation on u and v, I followed either 'by adding a term to the continuity equation to convert the I I I third to the same form as the other two equations so that the traditional relaxation can be used, or by using a pressure correction iteration to make continuity equation satisfied. The problem with these methods is I I that, not o4ly they tend to converge very slowly, but they may not always be consistent ;with ,the real physical problem in some cases. Here, w,l introduce a so-called distributive relaxation scheme, which is similar I I to but different from the distributive scheme developed by Brandtl121. 'I : :I Numerical show that our distributive relaxation method, when cou-pled with the multigrid technique, is much more efficient for solving incompressible equations. The basic idea is to start with traditional point ': 29

PAGE 40

Gauss-Seidel relaxation for the first and second equa.tions m (3.5), which we write a.s follows: (Su)p + CpPp Upf---. Ap (3.6) (Sv)p + DrPP Vp f--.... Bp (3.7) I We then ad,.iust up, vp and Pp simultaneously to satisfy the continuity equation I a.nd both :rhomentum equations, using the corrections given in the form new c Up =Up-ou, new c Vp = Vp-oV, The corrections of other components in the system are set to zero . I ; I tl.p {36 0 Pp+6P Vp6 Figure 3.1. Distributive relaxation. The defining equations for these corrections are (AE + Ap )6u-Cp6P 0, 30 (3.8)

PAGE 41

(3.9) where This yields f3 6u Cp BN+Bp /).y --="-'-, 6v Dp AE + Ap /).z 6v (Sm)p = +F;))7+(GN + Gp )' 6u {36v, 6P BN + Bp6v. Dp (3.10) Algorithm For the large grids and long-time evolution models considered here, the usual relaxation are much too inefficient by themselves. To obtain optimal I efficiency, w;e use a multigrid method based on the relaxation schemes developed I in the prev(ous section. For simplicity of discussion, we consider only the twogrid case as depicted in Figure 3.2. We use a full approximation scheme(FAS) to accommodate nonlinearities. I A FAS algorithm for an equation of the form I I 'I (3.11)

PAGE 42

. I I i I I : may he described loosely as follows: I I' uh +---uh + 1;h(u2h-I I (3.12) The notatio.n we have introduced includes the difference o.perators Lh and L2h, ':. the operators (for the approximation) and (for the residual), and the interpolation operator Figure 3.2. Two-level staggered grid. The following subsections describe this algorithm in more detail. Fine Grid !Relaxation and Evaluation : I ....... :U. . -. ,, We first distributive relaxation for (3.5), then calculate the residuals on the fine gri'd level: 32

PAGE 43

Rh = (Sv)ph DpP!'-.), v Rh =(Sm)ph m t,J t,J (3.13) Coarse Gtid Relaxation Relaxation bn the coarse grid requires restriction of the a.pproximations and residuals from fhe fine grid to the coarse grid to provide initial guesses and righthand-sides,: respectively. For restriction of approximations, we use the following stencils: (3.14) For restriction of residuals, we use the following full-weighting stencils, which come from the so-called area law developed by Liu[15)(Figure 3.3): [ : l 33

PAGE 44

I' 1 1 8 8 1 1 4 4 1 1 I' 8 8 I I II [ : l. I I I I I I uh I I Figure 3.3. Restriction for Ru.. ; I II I This leads, to the following coarse grid system: I I I 'I :I I I :I I 'i 'I '' 'i 'I 'I '' I I :I 1 = 1,] 1,] (A2h -+ A2h-+ A2h-+ A2h -EEUi+2,j E Ui+t,j WUi-l,j WWUi-2,j +A2h -+ A2h-+ A2h-+ A2hNNUi,j+2 N Ui,j+l S ssUi,j-2 2h 2h 2h 2h -Ap u + CwvP: 2 + Cw P. 1 -Cp p .. 1,] 1-,] -,] 1,] _ Dp2h = .,, .,, (B2h -+ B2h-+ B2h-+ B2h -EEVi+2,j E Vi+t,j wVi-1,j WWVi-2,j +B2h -+ B2h-+ B2h-. + B2hNNVi,j+2 N Vi,j+l S Vi,j-1 SSVi,j-2 34 (3.15) (3.16)

PAGE 45

2h-2h -2h -2h --Bp vi,j + D55 Pi,j-2 + D5 Pi,j-1Dp Pi,j (3.17) (F2h -+ p2h+ r;o2h-p2hEEUi+2,j E Ui+1,j rw Ui-1,j -p Ui,j +G2h -+ G2h-+ G2h-G2h-) NNVi,j+2 N Vi,j+1 S Vi,j-1 -p Vi,j (3.18) where u, v :and P denote restrictions of the respective u, v and P from the fine grid We then perform one or more relaxations on (3.16)-(3.18), using these values as the initial guesses. (In the pra.ctical multilevel I scheme, would recursively solve (3.16)-(3.18) by way of corrections from yet I coarser grifls.) Interpola'tiori of Correction to Fine Grid I The next is to interpolate the corrections from the coarse grid to the fine I grid: Uh +uh + 1;h(u)(u2hVh +Vh + 1:h( v )( v2h v )vh), ph +-Ph+ 1:h(P)(P2h{3.19) 35

PAGE 46

. I I I I Here we bmnear interpolation, given by the following stencils: I I I ' I I ,I I : i :I I Figure 3.4: the relations between u 2 h and uh, v 2 h and vh, I correspondipg to the above stencils. I (3.20) p2h and ph The FA,T V-cycle based on the ingredients described here performed well in all of the cases we have considered. Typical convergence factors per V(l,l) cycle I ( c.f., Branqt [12]) were about 0.3 for the fully-implicit scheme . I I II 36

PAGE 47

, I i I CHAPTER 4 TEST CASES In this: chapter, we first describe the initial conditions used for our experi, I I I ments, sh?w well how our schemes work for the cases of low, medium and high Reyndlds numbers. Initial Small Disturbances I ; To study transition, we considered the least stable eigenmodes of the linear I theory. FJr such modes, the linear theory provides an approximate a.nalytic '! solution (which is accurate for small disturbances) to the governing equation I (2.9)-(2.q) given by I 'I I : i I (4.1) where e ja small positive number and Real stands for 'the real part of'. The complex .P R + i,Pr and the temporal frequency w WR + iw1 come from solutton the Orr-Sommerfeld eigenvalue problem. The Orr-Sommerfeld equation that governs the linear stability of parallel shear flow is i 2 2 a?)2 iRe{(a:uo-a:)= 0,

PAGE 48

l/>(x,-1) = 4>(x,1) = l/>'(x,-1) = 4>'(x,1) = 0. (4.2) The problems chosen here for study used Re = 2000, 7500 and 20000, I I and a = 1:. Equation ( 4.2) was solved by a spectral method with 50th-order I I Chebyshev I polynomials. The only observed growing mode for each Reynolds number is shown in Table 4.1. The eigenfunctions for those modes are shown in Figure 4.1 (a), (b) and (c), respectively. Re WR W[ 2000 0.31210050 -0.01979860 7500 0.24989154 0.00223497 ........ ... ......... -;. i 20000 0.20966415 0.00330604 ..... -. Table 4.1 Temporal Frequency for Different Re. Note that the stream function is given by so, the and normal disturbance velocities can be obtained as Setting a ::!: 1, 'we obtain the eigenfunctions for both u and v in the following I form: q,u(y) l/>R. + ilj>( = l/>'(y ), 38

PAGE 49

+ = -i
PAGE 50

. eigenfunction of u 6 Real( ) u 0 -2 -4 -6 -1 0 1 I eigenfunction of v 0.06 Im( ) v 0.00 -0.02 -0.04 -0.06 -1 0 1 ' eigenfunction of v 2.5 1.5 Real( ) v 1.0 0.5 0 1 width of channel I Fi1gure 4.1.(a) Orr-Sommerfeld eigenfunction with Re=2000. 40

PAGE 51

6 4 2 0 -2 -4 0.06 0.04 0.02 0.00 -0.02 -0.04 rffi( ) u eigenfunction of u Real( ) u 0 1 eigenfunction of v Im( ) v -1 0 1 2.5 2.0 1.5 1.0 0.5 I : eigenfunction of v Real( ) v 0 width of channel 4.1.(b) Orr-Sommerfeld eigenfunction with Re=7500. 41 I I I 1

PAGE 52

1.5 1.0 0.5 0.0 -1 I :I I. I i I I ; i I I ' I I 'I eigenfunction of u 0 width of channel 4.1.( c) Orr-Sommerfeld eigenfunction with Re=20000. :t 42 1

PAGE 53

1:[ r I I I ; I I I Computational Results I We choie !lre, WR and wr to be the values shown in Table 4.1, which define the least stible:modes fur the corresponding Reynolds numbers. The parameters used in ou:rl experiments were b.t To L -= -/500, 500 w.,. L 27r, nJ 34, I I :I nt -66 and 130. \ Though other s.chemes, like semi-implicit second-order and semi-implicit partial i I i schemes, were also used, they did not yield satisfactory results. I For flrst' two schemes, second-order and y-direction fourth-order, we took :I I Re = 7500 and tested for accuracy. All three Reynolds numbers are used for the scheme. Distributive relaxation coupled with FAS multigrid proved efficient each of these schemes, with a convergence rate ranging between 0.3 I a.nd 0.4 It he V ( 1 ,1) cycle (see Figure 4.2). Figures 4.3-4.10 display the computational results. II I :I II ll I :I 43

PAGE 54

> "d 'il II "" 0 .... toG 0 ..t convergence history of u (during one time step) .... ---...... --multigrid one level grid I I work units o I 200 400 600 I (a) convergence history of v (during one time step) -7.0 ........... -7.5 .. ........ .......... ..... -8.0 multigrid -8.5 one level grid -9.0 work units 0 200 400 600 (b) I :convergence history of mass (during one time step) -4.0 ---:!. __ I ---1 multigrid one level grid work units 0 200 400 (c) Figure 4_.2. Convergence history for (Re = 7500, mesh: 34 x 130) (a) u, (b) v, (c) mass. 44 600

PAGE 55

-0 -......... ..... --tiLl 0 -tiLl 0 I I I I I 'I 'I I I T I I '' I :I I I ; I I I .I 0.!4 I ; I I I I I : i :I i I I DO x x I '' I lbeorelical solution 34 X 130 34 X 66 ,1' xi x 01.0 :o.o 0.5 1.0. 1.5 2.0 2.5 3.0 Re=7500, second-order fully-implicit T/To i :I I I I 'I :I 0.4 I' I. :I Figure 4.3. Disturbance energy increase (second-order fully-implicit, 3 periods) theoretical solution Oi 0 34 X 130 X X 34 x 66 :I o:. 0 ............. _._.......,..'-::-'-----'-_.__"-="' 10.o o.5 t.o 1.5 2.0 2.5 3.o Re=7500, cross-channel fourth-order T/To Figure 4.4. Disturbance energy increase (cross-channel fourth-order, 3 periods). : 1 45 'i :I I I 'I 'I I

PAGE 56

0 -10 :::::: 0 .... -1.5 .2 -2.0 0.4 0.0 0.0 0.8 0 ........ 2: 0.4 .2 oo X 0 1 1.0 2.0 (b) 0 0 0 the1oretfcal 34 130 34 ?' 66 X X X 1.0 2.0 (c) Figure 4.5. Disturbance energy increaose (fourth-order fuJly-implicit, 3 periods). 1 (a) Re = 2000, (b) Re = 7500, (c) Re = 20000. I 46 'I I I X X 3.0 T/To 3.0 T/To

PAGE 57

7{ T/To=O 7{ T/To=0.25 0.5 0.0 I i' 7{ T/To=0.5 7{ T/To=0.75 7{ T/To=l Figure 4.6. Contours of disturbance stream function (second-order fully-implicit, one period). 47 2n 2n 2n 2n

PAGE 58

0.5 0.0 -0.5 0 0.5 7T T/To=O 7T T/To=0.25 7T T/To=0.5 7T T/To=0.75 7T T/To=l Figure 4.7. Contours of disturbance stream function (cross-channel fourth-order, one period). 48 27T 2n 27T 27T

PAGE 59

1.0 0.5 o.o 1i T/To=O 1i T/To=0.25 1i T/To=0.5 1i T/To=0.75 -0.5 0 1i T/To=l Figure 4.8. Contours of disturbance stream function (fourth-order fully-implicit, one period). 49 211" 211" 211" 211" 211"

PAGE 60

0.5 1T Theoretical result, T/T0=20, llj01MAX=2.46E-04 fourth-order, T/T0=20, llj01MAX=2.46E-04 fourth-order in y-direction, T/T0=20, llj0111.u=2.13E-04 second-order, T/T0=20, lljOIMAx=1.42E-04 I Figure 4.9. Comparison of contours of disturbance stream function 'I ( Re=7500, time=20 periods). linear theoretical result, (b) fourth-order fully-implicit, (c) cross-channel fourth-order, (d) second-order fully-implicit. 1 so I 21T

PAGE 61

; dh3turbance ener distribution Re=7500 -.-.... 0 2.5 t;r 1.5 ......... -....., ----1.0 bll 0 -0.5 theoretical solution fourth-order fully-implicit cross-channel fourth -order second -order fully-implicit -----------------4.10. Comparison of different schemes with disturbance energy . I From the above results, we see clearly that the fourth-order fully-implicit scheme is lery accurate and efficient for simulating the process of transition, at least at small disturbance stage. We compare our results with the linear theory in Tables 4.2 and 4.3. theoretical fourth-order cross-channel second-order I results fully-implicit fourth-order fully-implicit I ... --... ...... -------average 0.2235 x 10-3 0.2245 X 10-3 0.1931 X 10-3 0.1223 x 10-3 .. relative error. 0.0% 0.45% 13.60% 45.28% :I I I Table 4.2. Comparison of the Increasing Rate of Disturbance Energy I (in 20 Periods). s1 I I .I

PAGE 62

, I I 'I '1 I theoretical fourth-order cross-channel second-order I ; I results fully-implicit fourth-order fully-implicit I I I ; I observeo value 2.46 X 10-4 2.48 x 10-4 2.1a x 10-4 1.42 x 10-4 :I 0.0% 0.81% 13.41% 42.27% :I phase 11aggiin.g I 0 0.02571" 0.27871" 0.26671" 4.3. I I Comparison of the Maximum Stream Function Values and Phase Lagging (after 20 Periods). '! I I I I i I Concluding Remarks From the above test cases, we conclude that I I I finite differences can be very accura.te for simulating transition. relaxation coupled with multigrid can accelerate the iteration I make the more accurate and sta.ble fully-implicit schemes com1 1 i I para,ble cost to the explicit schemes at each time step. I : I I ; The ;simple two-dimensional problem studied here shows promise as a very I I I: and feasible approach for more practical problems, including three ,I I dimensions and complicated regions, especially in the context of adaptive I fil re n:ement. I I I I :I 'I ,, I :I i I 52

PAGE 63

CHAPTER 5 AND OPEN QUESTIONS Synopsis of the Thesis I The questio!ls motivating this work are: Is the finite difference method accurate enough to simulate transition? Can multigrid efficiently increase the power of the finite-difference method? Can higher-order finite-difference schemes o:btain acceptable accuracy, enough to compare with spectral methods? In the work of Liu, Liu and McCormick (16], the first two questions were answered for the case of short-time evolution with Reynolds number Re = 7500. But the results were not fully satisfying, especially for the problem of phase lagging. The results presented in this thesis represent a significant advance over their results. In Chapter 2, we established three numerical schemes based on the staggered grid. Care must be taken in developing high-order schemes because u,_ v and P / .. are not represented at the same points. This task was somewhat simplified by the fact that we only considered the small disturbance problem. Chapter 3 described the tools for solving the equations developed in Chapter 2, including the so-called distributive relaxation method and the full approxima.tion scheme (FAS) of multigrid method. We then showed that these methods are

PAGE 64

easily applied simulate flow transition. I' Chapter 4 documents a first step in using the developed method to simulate the process o transition. The temporal problem was chosen for illustration. I r: Excellent perf<;)rmance was observed for the fourth-order fully-implicit scheme I over a range of Reynolds numbers. I We summarize our results as follows: I 1. Staggered grids and fully-implicit time-marching schemes can obtain high accuracy; 2. Multigrid methods substantially accelerate relaxation processes, making the more accurate and stable fully-implicit and semi-implicit time-ma.rching schemes comparable in cost to explicit techniques at each time step. : I 3. For smarl disturbances, finite difference schemes with high order accuracy achieve excellent agreement with the linear theory in both phase and am I plitude. Open Questions Though numerical simulation of transition has experienced tremendous progress I over the past decade ( c.f., Kleiser and Za.ng [17]), there are many remaining unsolved problems. In the special direction of our research, a number of avenues for further investigation appear to be open. 54

PAGE 65

The schemes developed in Chapter 2 are very effective for the simple model problems we have treated. They seem to promise a very efficient and fea.sible approach for more practical problems, including three dimensions and complicated regions, especially in the context of adaptive refinement. However, the relaxation method introduced in Chapter 3 can only be used for 2-dimensional problems, so a new relaxation method must be developed for three dimensions. Like line relaxation for 2-dimensional problem, perhaps a plane relaxation method may be used to obtain better smoothing. For multigrid, although there seems to be no apparent difficulty in extending it to more complicated problems, we must pay special attention to-the boundary points. In Chapter; 3, we only used the temporal problem in our tests. For more practical cases, we should also consider the spatial problem, where the periodic boundary condition in the x-direction would be changed to a Dirichlet condition up-stream and. a buffer domain for the non-reflecting outflow boundary ( c.f., Danabasoglu, Biringen and Streett [18]). The main difficulty here will be that the computational domain will be 20 or 30 times as large as the one we have taken in this tnesis, and local refinement should be used downstream to simulate the details of tJ;"ansition. These a.re hut a few of the directions for further work that appear to be promising. study of multilevel adaptive methods for transition flow with I roughness elements, an area which is almostly unexplored, appears to be a rich area for continuing research. 55

PAGE 66

BIBLIOGRAPHY [1] 0. Reynolds. An experimental investigation of the circumstances which determine whether the motion of water should be direct or sinuous, and of the law of resistance in channels. In Phil. Trans. Roy. Soc., A 174, page 935, 1883. [2] Lord Rayleigh. On the stability, or instability, of fluid motions. In Scientific Papers of Lord Rayleigh 1, pages 474-487, New York: Dover, 1964. I [3] W.M.F. Qrr. The stability or instability of the steady motions of a perfect liquid and:'of a viscous liquid. Proc. Roy. Irish. Acad. 27, pages 9-69, 1907. [4] A. Sommerfeld. Ein Beitrag zur hydro-dynamicschen Erklarung der turbulenten Flussigkeitsbewegung. In Atti. Gongr. Int. Math., ,4th Rome, pages 116-124, 1:908. [5) W. Tollmein. Uber die Entstehung der Turbulenz. Nachr. Ges. Wiss. Gottingen ;Math.-Phys. Kl., pages 21-44, 1929. [6) H. Schlichting. Zur Enstehung der Turbulenz bei der Platten-stromung. In Nachr. Wiss. Gottingen Math.-Phys. Kl., page 182, 1933.

PAGE 67

(7] B.J. Bayly and S.A. Orszag. Instability mechanisms in shear-flow tra.nsition. In Ann. Rev. Fluid Mech. pages 359-391, 1988. [8] J.-P. Zahn, J. Toomre, E.A. Spiegel and D.O. Gough. Nonlinear cellular motions in Poiseuille channel flow. J. Fluid Mech. 64, pa.ges 319-345, 1974 . [9] S.A. Orszag and A.T. Patera. Subcritical transition to turbulence in plane channel Phys. Rev. Lett. 45, pages 989-993, 1980. [10] T.A. Zang and M.Y. Hussaini. On spectral multigrid methods for the timedependent Navier-Stokes equations. Appl. Math. Comput. 19, pa.ges 359-372, 1986. ! (11] F.H. Harlow and J.E. Welsh. Numerical calculation of time-dependent viscous incompressible flow with free surface. Phys. Fluids B, pages 2182-2189, 1965. [12] A. Multigrid techniques: Guide with application to fluid dynamics. GMD GMD, St. Augustine, 1984. (13] A.J. Chorin. A numerical method for solving incompressible viscous flow problems. J. Comput. Phys. 2, pages 12-26, 1967. (14] 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 I Transfer 15, pages 1787-1806, 1972. 57

PAGE 68

[15] C. Liu. adaptive methods in computa.tional fluid dynamics. Ph.D. I thesis, of Colorado at Denver, i989. I :: [16] C.Liu, Z.Liu and S.McCormick Multigrid methods for flow transition in planar channel. Comput. Phy. Comm. 65, pages 188-200, 1991. I [17] L. Kleiser T.A. Zang. Numerical simulation oftransition in wall-bounded I shear flows. Ann. Rev. Fluid Mech., Vol. 23, 1991. i [18] S. Biringen and C.L. Streett. Numerical simulation of I spatially-evolving instability. !CASE Transition and Stability Workshop Pro' ceeding.s, i989. ; 58