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
23 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 fullyjimplicit time steps on a staggered grid, we develop a secondorder
scheme, a secondorder in the streamwise direction and fourthorder in the normal
i
direction scheme:, and a full fourthorder scheme. A very efficient timemarching
procedure usingj a multigrid solver is shown to produce very accurate results.
The computational results show that the full fourthorder scheme gives very good
agreement in both amplitude and phase with the linear theory obtained by solving
the OrrSom!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
SecondOrder FullyImplicit Finite Difference Scheme
i
CrossChannel FourthOrder Finite Differences ....
Fourth'Order FullyImplicit 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 xmomentum 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 uequation...........................
Fourthorder approximation for vatau point (interior points).
Neighbor points of vp for vequation...........................
Neighbor points of Pp for continuity equation..................
!
Approximation of v on solid wall boundary points...............
Distributive relaxation.............
Twolevel 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) OrrSommerfeld eigenfunction
(b) OrrSommerfeld eigenfunction
(c) OrrSommerfeld 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
(secondorder fullyimplicit, 3 periods)
Disturbance energy increase
(crosschannel fourthorder, 3 periods).
Disturbance energy increase
(fourthorder fullyimplicit, 3 periods).
(a) Re = 2000, (b) Re = 7500, (c) Re = 20000
. 
Contours of disturbance stream function
(secondorder fullyimplicit, one period). . .
Contours of disturbance stream function
(crosschannel fourthorder, one period).
Contours of disturbance stream function
(fourthorder fullyimplicit, one period).
Cpmparison of contours of disturbance stream function
(Re = 7500, time=20 periods)
(a) linear theoretical result, (b) fourthorder fullyimplicit,
(c) crosschannel fourthorder, (d) secondorder fullyimplicit. .
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
OrrSommerfeld 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 vortexshaped streamlines. The vortex moves 3dimensionally 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 socalled 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 timemarching and non
staggered grids (c.f., Zang and Hussaini [10]). The difficulties with numerical
simulation for transition result from the longtime 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 timemarching 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 timemarching. 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 2D 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 fourthorder fullyimplicit scheme.
I
CHAPTER 2
I
GOVERNING EQUATIONS FOR 2D 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 twodimensional, timedependent, incompressible NavierStokes 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 halfwidth, 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. (25)
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 xdirection:
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 timedependent. 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)
SecondOrder FullyImplicit 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 secondorder central differ
I Figure 2.2. Staggered grid.
First consider the xmomentum 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 fullyimplicit finitedifference
scheme:
3
4.  2v,v"+ +
l 2At 0 2As Vl x'3
+(
+(
1l7l+1 _L 7n+1 7
?+l,j!+ UiJ Ui+l,j + Uij + ui,j Uihi +
)/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,jl + 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 xmomentum 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 = PiU
' 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 + )
*)>
05(."+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,j1>
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(<+ +
05i1 +
oss,+<;>,
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;{vV) + 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
= unil ,j)
vu = vnil,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 xdirection (for one wave length).
On the solid wall boundary, we use a modification of the xmomentum 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 +Vjij+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 ydirection, we use
i A^.ew = A}d _________
,! As 5 ' Re 3 Ay2
i new A old , ^ ^
lp ~Ap +Tuw'
13
(2.22)
CrossChannel FourthOrder Finite Differences
i
i
Since th'e solution tends to change more rapidly in the ydirection, either re
i '
1
fined gridsjor higher order differences in the ydirection can be used to balance
overall accuracy. One of the methods we study in this paper uses fourthorder
central differences in the ydirection for interior points only, while maintaining
secondorder central differences for boundary points and for interior points in
the xdirection. We use a fullyimplicit scheme, since it is generally more accu
i '
i
rate than semiimplicit or explicit methods. For simplicity, we continue to use
 !
j
secondorder differences for the continuity equation and the pressure terms in
' l
both momjentum equations. The general form for the fullyimplicit 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,i2
In (2.23), if we define
then
i i
ue 0.5t + <.),
uw 0.5> + ),
Un 0.5t> +
V. = 0.5t + ?+_>,),
vn = 05Ka+t,?_YJ+1),
^nn OSMJi + '%),
 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(rs + tit;),
Cw, = \ Cp =
2 At 2Ax
1
48Ay
ReK Ax2 2 A y2
Ax'
4 id* + u? 1
Similarly,
define
2At
Ve = oSMj1 +
Vw = 0.5+ + v&'j),
Vnn = + &!.),
vn = oW+O
V. = 0.5f +
= 0.5K/2 + "A).
= 05m + tfisi)
= osKi' + 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.) + (,tt + rrr),
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 secondorder scheme is maintained at the solid boundary points, the
boundary modifications are the same as the previous section.
i :
I i
. i FourthOrder FullyImplicit 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 crosschannel
' !
i ;
fourthorder fschjeme, they are inadequate for simulating transition because of the
sensitivity of suph flows to small changes and the need for longtime evolution.
For these reasons, we developed a socalled fourthorder fullyimplicit 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 xdirection.
! i
The general form of the fourthorder fullyimplicit finitedifference 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;  BunVnn + 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 ~ ui2,j,
Vi+2J,
Pi2,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
(Vt2,jX + Vx2,j+2 + U+l,jl + 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
^r2y,Vj
2AÂ£ ^
o  1 O
o 'U'WW  1 PP Uw 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 uequation.
21
Similarly, for the vequation (2.27) (see Figure 2.10), let
v.e =
Uiiillt
05(u+i j + i++2j),
,n+l
u.
Wllltrt 
os(nfij+
0.5(u{*2j f i+2,ji)>
0.5(uJt11J + "J^jx)
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 fourthorder 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 Vll,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 ' iji iiJ ij 1 iji ' *42 iji \
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) 
W2J} + ^*"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,jl Vi+l,j + Vt+2,j + M.+ljl 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,j1 ViJ + VilJ Uil,j + Vil,j Vi2,j \
. 21 2 2 2 ;
f "12T+ 2722
i^l1 + ytfi. 1 +
1 &
_L 7in+1 7in+1 4 7in+1 7j+1 4 ItV't1
+ ViJl ViJ~ 1 + V.J2 ^t.jl + V,j2 v
!i 2 2 2 '
&+i8a*j 3Qv?r+tana a*j
! 12A2
 ie
12Ay2 ]
H27P"/1 27
1 24 Ay
r^t1
't.j2
= 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 :  Vnil,j,
Plj  P 1 ml,j j
uo j  Uni2 ,j )
Vo J = Vni2,j j
Poj Pni2,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 xdirection (for one wave length).
On the 'solid; wall boundary, we change the ydirection difference to second
order, and maintain fourth order in the xdirection. Again, special attention
25
! .
i
should be j>aid to the approximate value of v. Letting (see Figure 2.12)
j+i)
Wi + + 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 ydirection become
I
Vi2,j+l t. vil,j+l \ Vi,j+1  vi+l,j+l \
Ii 1 f  . 1 1 
,1 1' ^2,i ViU 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 fourthorder 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 righthand side of the above system is written as
' I '
 ;
= (Aeeu{+2j + AEUi+i'j + Awiiiij + Aww^ii,j
: i i
j,j  +ApfpfUij+2 f ApfUiij+1 ( A$Ui j_l + AssUi,j2
i i 1
: j \CwwPi2,j + CwPii,j + CeP{+ij) + {Su)p,
(Sv)p = (BEEVi+2,j + BEVi+l,j + B\yVil,j + BwwVi2,j
\ ;
i i
j j +BffNUij+2 + BpiVi,J+i + BsV{,ji f BSsVi,j2
I 1 +DssPi,j2 + DsPi,j\ + DffPi,j+i) ( (5v)p,
L
(Sjn)p = (FEEui+2,j + FEui+u + FwUii,j
h. !
1 : !
, I
i I
i i 1
+Gjvj\fViiJ+2 + GrtUi'j+i + GsVijx) + (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 Gau6sSeidel
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
GaussSeidel 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 socalled 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
ible NavierStokes equations. The basic idea is to start with traditional point
29
i
i
GaussSeidel 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,
 (310)
Up
Multigrid Algorithm
For the large grids and longtime 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 twole;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. Twolevel 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
handsides,, 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 fullweighting stencils, which
come from: the socalled 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?wUi1,5 +
I +^5wJ+2 + + A2shsUitj_2
i i
! '
^4plU;ij + Cw\V Pi2,i + CftPiU 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,j1 Dp1 Pij
+D2NhPitj+l) + Phh(v)Rv + (Sv)2ph,
nh<5 O =
~{^EhE^i+2,j + P^^i+l ,3 + Fwil,j ^Pht,j
+ ^MVt,j+2 + OiJ+l + G2ShVi,j1 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)(P2hPhh(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 Vcycle 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 fullyimplicit 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)eiaxiu,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 OrrSommerfeld eigenvalue problem. The OrrSommerfeld
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 50thorder
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)SaxWRt).
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), (44)
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) OrrSommerfeld eigenfunction with /Ze=2000.
40
eigenfunction of u
I
Figure 4.1.(b) OrrSommerfeld 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 semiimplicit secondorder and semiimplicit partial
i i
i
fourthorder schemes, were also used, they did not yield satisfactory results.
For the first two schemes, secondorder and ydirection fourthorder, we took
; i
Re = 7500 and tested for accuracy. All three Reynolds numbers are used for the
I '
fourthorder 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.34.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
(secondorder fullyimplicit, 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
(crosschannel fourthorder, one period).
48
T/To=0.25
7T
T/To=0.5
T/To=0.75
Figure 4.8. Contours of disturbance stream function
(fourthorder fullyimplicit, one period).
49
I
Theoretical result, T/To=20, ^MAX=2.46E04
fourthorder in ydirection, T/T0=20, ^LaX=2.13E04
! I
0 ^ 7T 2n
1 ] secondorder, T/T0=20, ^MAX=1.42E04
Figure 4.9. Comparison of contours of disturbance stream function
( Re=7500, time=20 periods).
(a) linear theoretical result, (b) fourthorder fullyimplicit,
(c) crosschannel fourthorder, (d) secondorder fullyimplicit.
I 50
From the above results, we see clearly that the fourthorder fullyimplicit
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 fourthorder crosschannel secondorder
, ' results fullyimplicit fourthorder fullyimplicit
average value 0.2235 xlO3 0.2245 xlO3 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 fourthorder crosschannel secondorder
1 i i i ; results fullyimplicit fourthorder fullyimplicit
; 1 observed value 2.46 x 10"4 2.48 x 10~4 2.13 x 104 1.42 x 104
:  ; 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
Highorder 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 fullyimplicit schemes com
I
' i ;
parable in cost to the explicit schemes at each time step.
The simple twodimensional 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 finitedifference method? Can higherorder finitedifference
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 shorttime 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 highorder 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 socalled 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 fourthorder fullyimplicit scheme
I
over a range of Reynolds numbers.
I
We summarize our results as follows:
1. Staggered grids and fullyimplicit timemarching schemes can obtain high
accuracy.
2. Multigrid methods substantially accelerate relaxation processes, making the
i
more accurate and stable fullyimplicit and semiimplicit timemarching
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 2dimensional problems, so
a new relaxation method must be developed for three dimensions. Like line re
laxation for the 2dimensional 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 xdirection would be changed to a Dirichlet condition
upstream and a buffer domain for the nonreflecting 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 474487, 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 969, 1907.
[4] A. Sommerfeld. Ein Beitrag zur hydrodynamicschen Erklarung der turbu
lenten Flu'ssigkeitsbewegung. In Atti. Congr. Int. Math., ^th Rome, pages
116124, 1908.
[5] W. Tollmein. f/ber die Entstehung der Turbulenz. Nachr. Ges. lYiss.
Gottingen \Math.Phys. Kl., pages 2144, 1929.
[6] H. Schlichting. Zur Enstehung der Turbulenz bei der Plattenstromung. In
Nachr. Ges. Wiss. Gottingen Math.Phys. Kl., page 182, 1933.
I
[7] B.J. Bayly and S.A. Orszag. Instability mechanisms in shearflow transition.
In Ann. Rev. Fluid Mech. 20, pages 359391, 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 319345, 1974.
l
[9] S.A. Orszag and A.T. Patera. Subcritical transition to turbulence in plane
channel flows. Phys. Rev. Lett. 45, pages 989993, 1980.
[10] T.A. Zang and M.Y. Hussaini. On spectral multigrid methods for the time
dependent NavierStokes equations. Appl. Math. Comput. 19, pages 359372,
1986.
i

[11] F.H. Harlow and J.E. Welsh. Numerical calculation of timedependent vis
cous incompressible flow with free surface. Phys. Fluids 8, pages 21822189,
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 1226, 1967.
i
[14] S.V. Patankar and D.B. Spalding. A calculation procedure for heat, mass
and momentum transfer in 3dimensional parabolic flows. Int. J. Heat Mass
Transfer 15, pages 17871806, 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 188200, 1991.
[17] L. Kleiser and T.A. Zang. Numerical simulation of transition in wallbounded
shear flows. Ann. Rev. Fluid Mech., Vol. 23, 1991.
i
j,
[18] G.Danabasoglu, S. Biringen and C.L. Streett. Numerical simulation of
spatiallyevolving instability. ICASE Transition and Stability Workshop Pro
ceedings, 1989.
i
58

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 secondorder scheme, a in the streamwise direction and fourthorder in the normal direction and a full fourthorder scheme. A very efficient timemarching procedure a multigrid solver is shown to produce very accurate results. The results show that the full fourthorder 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
:II ,, 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 FullyImplicit Finite Difference Scheme I 'I I CrossChannel FourthOrder Finite Differences . I .. FullyImplicit 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 xmomentum 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 uequation. 7 8 9 10 12 13 19 2.9 Fourthorder approximation for v at a u point (interior points). 20 :I 2.10 N;eighbor points of Vp for vequa.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 Tlolevel 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) OrrSommerfeld eigenfunction with Re =7500. 4.1. 4.2 4.3 4.4 4.5 4.6 4.7 4.8 (c): OrrSommerfeld 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 (secondorder fullyimplicit, 3 periods) ............... i energy increase ( cJosschannel fourthorder, 3 periods). I '' energy increase (f9urthorder fullyimplicit, 3 periods). (a) Re.= 2000, (b) Re = 7500, (c) Re = 20000 .. '! Contours of disturba.nce stream function (secondorder fullyimplicit, one period). :I cJntours of disturbance stream function I (crosschannel fourthorder' one period). clntours of disturbance stream function 'I I (fourthorder fullyimplicit, 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) fourthorder fullyimplicit, 40 41 42 44 45 45 46 47 48 49 ( t) crosschannel fourthorder, (d) secondorder fullyimplicit. 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 OrrSommerfeld 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,ortexshaped streamlines. The vortex moves 3dimensionally 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 socalled 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 timemarching and non sta.ggered grids (c.., Zang and Hussaini (10]). The difficulties with numerical I I' simulation for transition result from the longtime 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 timemarching 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 timemarching. 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 2D 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 fourthorder fullyimplicit scheme. :I II I o I o 1 [ I I 3
PAGE 14
CHAPTER 2 I GOVERNING EQUATIONS FOR 2D 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 twodimensional, timedependent, incompressible Na.vierStokes 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 Xand ydirections, respectively, I Pis the pressure a.nd Re is the Reynolds number based on the centerline velocity, '' '' U, of flow, the channel halfwidth, 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 xdirection: 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) = 1y2 (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 timedependent. 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). FullyImplicit 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 secondorder 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 xmomentum 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 theyvelocity component defined at au point, then we discretize I I equation on the above staggered grids using a fullyimplicit finitedifference 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 Uil,j ui+l,j 0 0 0 pil,j p,.. t,] Vil,j v. t,] 0 Uij1 0 0 I 2.3. Grid neighbors for xmomentum equation. 'i ; I (2.14) For simppcity, we may rewrite (2.14) in the general form (see Figure 2.4) ; I : AEuE Awuw + ANuN + AsusApupCpPp + 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 Ui1,j! UN Ui,j+l! us Ui,jt, Up Ui,j, Pp Pi,;, Pw pil,j UN 1uw Up UE 0 10 11Pvv 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 Vil,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. . vil,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 )( ) Up1 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 Vil,jl VN Vi,j+b vs Vi,jb Vp Vi,j, Pp Pi,;, Ps Pi,j1 0 5( n+l + n+l ) ui+1,j ui+l,j1 0 5( n+ 1 + n+l ) ui,i ui,i1 0 5( n+l + n+l) vi,j+t vi,; ; 0 5( n+l + n+l) vi,j1 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_p1 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 xdirection (for one wave length). On the wall boundary, we use a. modification of the xmomentum equal tion (Figure 2.7). For example, for j = 2 we use the approximations 82u::::::: _4_('1LNUp_ 0 = _ 4_Up, 8y2 3L\y L\y lAy ) 3 Ay2 Ay2 (2.19) 8uv ::::::: _1_( + t1il,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 ydirection, 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
CrossChannel FourthOrder Finite Differences 'I Since tends to change more rapidly in the ydirection, either re I I fined grids[ higher order differences in the ydirection can be used to balance overall a.ccJracy. One of the methods we study in this paper uses fourthorder central difflences in the ydirection for interior points only, while maintaining secondorder central differences for boundary points and for interior points in :I I the We use a fullyimplicit scheme, since it is generally more accu1 : I rate than semHmplicit or explicit methods. For simplicity, we continue to use secondordl dikerences for the continuity equation and the pressure terms in .I I I both equations .. The general form for the fullyimplicit scheme in i I this case be written as follows: I I I '' Apup + CwPwCpPp = Su, iBEvE + Bwvw + BNNVNN + BNvN + Bsvs + BssVss I :I Bpvp + DsPsDpPp = 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&,j2 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,,; ui1,; 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 vit,i+1 0 5( n+1 + n+l ) Vnn Vi,j+2 Vil,j+2 v. 0 5( n+ 1 + n+ 1 ) vi,j vi1,j 0 5( n+ 1 + n+ 1 ) vi,j1 vit,j1 ii .J. 0 25( n+1 + n+1 + n+1 + n+l. ) vi,; vi,j+1 vi1,j vi1,j+t 1 1 ilef1x2 2t1x (uoi,j + ue), 1 1 lfet1x2 + 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 vit,; 0 5( n+l+ n+l ) V V "+t 7 .,, t,J v, = 0 5( n+l + n+l ) vi,i vi,j1 Vu 0 5( n+l + n+l ) V. 2 V. 1 ':t,Ju., 0 5( n+l + n+l ) ui+1,j ui+l,j1 0 5( n.+l + n+1 ) Uw ui,j ui,j1 .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 n1 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. FourthOrder FullyImplicit 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 crosschannel I : I fourthordelschfme, they are inadequate for simulating transition because of the I sensitivity su1=h flows to small changes and the need for longtime evolution. I 1 For these we developed a socalled fourthorder fullyimplicit 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 xdirection. ; I I The form of the fullyimplicit finitedifference 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 Ui2,j, VEE Vi+2,j, vww Vi2,j, Pww pi2,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.. 1u 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 uequation '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 + . + 1l,J I2,J 1_l,J I2,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 yvelocity component defined at a u point), we must ai.lrolmate them with greater accuracy. Standard fourthorder central differences:Layibe used to obtain fourthorder 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 'Pi2,j1 ', 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' ,, Vil,j+l t I t I Vil,j v. l,] Vi,j+l J I t I v. l,] t I Vi+l,j+2 Vi+l,j1 t I I Figure Fpurthorder 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 {vnn27vn), .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 vequation. 21
PAGE 32
Similarly, for the vequation (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 Vil,j vi2,j 0 5( n+l + n+l ) Uee Ui+2,j Ui+2,j1 __:_ + (uee + 4uoij), 4 1 + (uee27ue32uoi,;), 4 1 + (27uwUww + 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 f0 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 fourthorder 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 + n1 n+ 1 + 8 n+ 1 8 n+ 1 + n+ 1 ui,i : ui,j ui,j 7J.i+2,j ui+1,j uit,j ui2,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 : t1,J '' tl,J + t1,J t2,J tl,J I2,J) ': 21 2 2 2 I 1 I + v!'+lt'+2 + v!'t+t2 + v!'+lt'+t + v!'t+tl +' .( ,t,J t1J 11J . 1 1J + 27 t,J I,J 1tJ t,J 246yl : 2 2 2 2 + n+l + n+1 nf1 + n+l n+1 + n+1 _27 ''! : a,J1 vit,i vi,j + ... ui,j2 vi1,j1 1'i.il) 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 ui2,j Re I +: + + 11J. 11J 1 1J 11J11J) :: 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 )] vi2,i1 'l?i2,i+2 vi+1,j1 vi+1,i+2 +i27 P!'71 27 p!'+t. + + 1 11J .t11 ) t2,J Q I = l I 24Ll:z: (2.29) : i :I 3 n+l '14 + n1 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 '26t 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,j1 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,J1 vi,j vi1,j + uit,; uil,j1 V&1,j vi2,;) 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,Jl t,J 11)1 + 11]1 11J2 1 1Jl t1J2) f 2; 2 2 2 1 + 16v!'++ 1 1 30v!'T1 + 16v!'+1 1 v!'+2 1 ( ,: I t' I 1J . 11J 11J Re 126.:z:2 + [6v!'T+1 1 30v'!"t1 + l61/'T1 1 v!'t1 2 + 1,,,, '' t,JloJ.) : 126.y2 I I I + Pi7/_;.,\ + :27 27 Pi7f\ + ui,j!2 = 0, '! 246y (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,j1 _.....::;_:......._; __ = 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 'UnilJ, Vnil,j, Pnil,j, '1Lni2,j, Vni2,j, Pni2J, u2,j, v2,j, P2,;, u3,j, v3,j, P3,;, (2.31) (2.32) where ni is ihe number of grid points in the xdirection (for one wave length). On the solid: wall bounda.ry, we change the ydirectlon difference to second order, and mai11tain fourth order in the xdirection. 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 nj2, of the discrete equations in theydirection become 'I i I 'I I I :I Vi2 j+1 Vi1 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 Vi1,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 nj1, 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: nj1, 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 fourthorder 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 righthand 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 GaussSeidel 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 socalled 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 coupled with the multigrid technique, is much more efficient for solving incompressible equations. The basic idea is to start with traditional point ': 29
PAGE 40
GaussSeidel 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 =Upou, new c Vp = VpoV, 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 )6uCp6P 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 longtime 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(u2hI 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. Twolevel 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 righthandsides,: respectively. For restriction of approximations, we use the following stencils: (3.14) For restriction of residuals, we use the following fullweighting stencils, which come from the socalled 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 WUil,j WWUi2,j +A2h + A2h+ A2h+ A2hNNUi,j+2 N Ui,j+l S ssUi,j2 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 wVi1,j WWVi2,j +B2h + B2h+ B2h. + B2hNNVi,j+2 N Vi,j+l S Vi,j1 SSVi,j2 34 (3.15) (3.16)
PAGE 45
2h2h 2h 2h Bp vi,j + D55 Pi,j2 + D5 Pi,j1Dp Pi,j (3.17) (F2h + p2h+ r;o2hp2hEEUi+2,j E Ui+1,j rw Ui1,j p Ui,j +G2h + G2h+ G2hG2h) NNVi,j+2 N Vi,j+1 S Vi,j1 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 Vcycle 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 fullyimplicit 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 OrrSommerfeld eigenvalue problem. The OrrSommerfeld equation that governs the linear stability of parallel shear flow is i 2 2 a?)2 iRe{(a:uoa:)= 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 50thorder 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) OrrSommerfeld 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) OrrSommerfeld 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) OrrSommerfeld 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 semiimplicit secondorder and semiimplicit partial i I i schemes, were also used, they did not yield satisfactory results. I For flrst' two schemes, secondorder and ydirection fourthorder, 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.34.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, secondorder fullyimplicit T/To i :I I I I 'I :I 0.4 I' I. :I Figure 4.3. Disturbance energy increase (secondorder fullyimplicit, 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, crosschannel fourthorder T/To Figure 4.4. Disturbance energy increase (crosschannel fourthorder, 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 (fourthorder fuJlyimplicit, 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 (secondorder fullyimplicit, 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 (crosschannel fourthorder, 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 (fourthorder fullyimplicit, one period). 49 211" 211" 211" 211" 211"
PAGE 60
0.5 1T Theoretical result, T/T0=20, llj01MAX=2.46E04 fourthorder, T/T0=20, llj01MAX=2.46E04 fourthorder in ydirection, T/T0=20, llj0111.u=2.13E04 secondorder, T/T0=20, lljOIMAx=1.42E04 I Figure 4.9. Comparison of contours of disturbance stream function 'I ( Re=7500, time=20 periods). linear theoretical result, (b) fourthorder fullyimplicit, (c) crosschannel fourthorder, (d) secondorder fullyimplicit. 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 fourthorder fullyimplicit crosschannel fourth order second order fullyimplicit 4.10. Comparison of different schemes with disturbance energy . I From the above results, we see clearly that the fourthorder fullyimplicit 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 fourthorder crosschannel secondorder I results fullyimplicit fourthorder fullyimplicit I ... ... ...... average 0.2235 x 103 0.2245 X 103 0.1931 X 103 0.1223 x 103 .. 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 fourthorder crosschannel secondorder I ; I results fullyimplicit fourthorder fullyimplicit I I I ; I observeo value 2.46 X 104 2.48 x 104 2.1a x 104 1.42 x 104 :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 fullyimplicit schemes com1 1 i I para,ble cost to the explicit schemes at each time step. I : I I ; The ;simple twodimensional 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 finitedifference method? Can higherorder finitedifference 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 shorttime 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 highorder 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 socalled 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 fourthorder fullyimplicit scheme I over a range of Reynolds numbers. I We summarize our results as follows: I 1. Staggered grids and fullyimplicit timemarching schemes can obtain high accuracy; 2. Multigrid methods substantially accelerate relaxation processes, making the more accurate and stable fullyimplicit and semiimplicit timema.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 2dimensional problems, so a new relaxation method must be developed for three dimensions. Like line relaxation for 2dimensional 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 tothe 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 xdirection would be changed to a Dirichlet condition upstream and. a buffer domain for the nonreflecting 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 474487, 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 969, 1907. [4] A. Sommerfeld. Ein Beitrag zur hydrodynamicschen Erklarung der turbulenten Flussigkeitsbewegung. In Atti. Gongr. Int. Math., ,4th Rome, pages 116124, 1:908. [5) W. Tollmein. Uber die Entstehung der Turbulenz. Nachr. Ges. Wiss. Gottingen ;Math.Phys. Kl., pages 2144, 1929. [6) H. Schlichting. Zur Enstehung der Turbulenz bei der Plattenstromung. In Nachr. Wiss. Gottingen Math.Phys. Kl., page 182, 1933.
PAGE 67
(7] B.J. Bayly and S.A. Orszag. Instability mechanisms in shearflow tra.nsition. In Ann. Rev. Fluid Mech. pages 359391, 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 319345, 1974 . [9] S.A. Orszag and A.T. Patera. Subcritical transition to turbulence in plane channel Phys. Rev. Lett. 45, pages 989993, 1980. [10] T.A. Zang and M.Y. Hussaini. On spectral multigrid methods for the timedependent NavierStokes equations. Appl. Math. Comput. 19, pa.ges 359372, 1986. ! (11] F.H. Harlow and J.E. Welsh. Numerical calculation of timedependent viscous incompressible flow with free surface. Phys. Fluids B, pages 21822189, 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 1226, 1967. (14] S.V. Patankar and D.B. Spalding. A calculation procedure for heat, mass and momentum transfer in 3dimensional parabolic flows. Int. J. Heat Mass I Transfer 15, pages 17871806, 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 188200, 1991. I [17] L. Kleiser T.A. Zang. Numerical simulation oftransition in wallbounded I shear flows. Ann. Rev. Fluid Mech., Vol. 23, 1991. i [18] S. Biringen and C.L. Streett. Numerical simulation of I spatiallyevolving instability. !CASE Transition and Stability Workshop Pro' ceeding.s, i989. ; 58
