Citation

## Material Information

Title:
Multigrid methods for simulation of flow transition
Creator:
Liu, Zhining
Publication Date:
Language:
English
Physical Description:
xxix, 286 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.
General Note:
Submitted in partial fulfillment of the requirements for the degree, Doctor of Philosophy, Applied Mathematics.
General Note:
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:
28863589 ( OCLC )
ocm28863589
Classification:
LD1190.L622 1993d .L58 ( lcc )

## Auraria Membership

Aggregations:
Auraria Library
University of Colorado Denver Theses and Dissertations

Full Text
MULTILGRID METHODS FOR
SIMULATION OF FLOW TRANSITION
by
Zhining Liu
B.S., University of Science and Technology of China, 1986
M.S., University of Colorado at Denver, 1991
A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
Applied Mathematics
1993

1993 by Zhining Liu

This thesis for the Doctor of Philosophy
degree by
Zhining Liu
has been approved for the
Department of
Mathematics
Chaoqun Liu
Thomas Russell

Date

Liu, Zhining
(Ph.D
Applied Mathematics)
Multigrid Methods for Simulation of Flow Transition
Thesis directed by Professor Steve McCormick
ABSTRACT
The purpose of this study is to develop a fast and accurate multigrid solver
for simulating the flow transition process. The research identified several tasks,
including 2-D and 3-D simulations with various specific conditions.
The efficient multigrid full approximation scheme (FAS) has been successfully
applied to perform direct numerical simulation (DNS) of flow transition in 2-D
and 3-D planar channels, 2-D and 3-D flat plate boundary layers, and 3-D chan-
nel and 3-D boundary layers with 2-D and 3-D isolated and randomly-distributed
multiple roughness elements in a curvilinear coordinate system. The algorithms
developed include a fourth-order finite difference technique on stretched and
staggered grids, a fully-implicit time-marching scheme, a full-coarsening multi-
grid method based on point distributive relaxation, a semi-coarsening multigrid
method based on line distributive relaxation, and a new treatment of the outflow
boundary condition that needs only a very short buffer domain to damp all wave
reflection.
These approaches made the multigrid DNS code very accurate and efficient.
It enabled not only spatial DNS for 3-D channel and flat plate flows at low

computational costs, but also spatial DNS for transition in 3-D boundary layers
with randomly-distributed and general shaped 3-D multiple roughness elements.
This represents research never before accomplished. Numerical results show good
agreement with linear stability theory, secondary instability theory, and a number
of laboratory experiments.
This study shows that the above approaches on spatial DNS of flow transi-
tion are very successful. The discretization schemes used are highly accurate and
stable, the multigrid method is very efficient, the new version of the governing
equations in general coordinates successfully performs DNS for transitional flow
with roughness elements, and the basic approach is extendible to problems with
more general geometry. Also, the new treatment of the outflow boundary dra-
matically reduces computational cost, making realistic spatial simulation much
more feasible.
The form and content of this abstract are approved. I recommend its publication.
Signed
Steve McCormick
IV

To Mom

CONTENTS
List of Figures............................................... xi
List of Tables...............................................xxvii
Acknowledgements...........................................xxviii
1 INTRODUCTION................................................. 1
General Overview............................................. 1
Advances Represented by This Study........................... 6
2 GOVERNING EQUATIONS IN PHYSICAL SPACE....................... 10
Original Form............................................... 10
Nondimensional Forms for Planar-Channel-Type Flow ............ 12
Nondimensional Forms for Flat-Plate-Type Flow............... 15
Linear Stability Theory..................................... 16
Blasius Similarity Solution................................. 20
Initial and Boundary Conditions for Planar-Channel-Type Flow .... 23
Initial and Boundary Conditions for Flat-Plate-Type Flow...... 25
3 TRANSFORMATION OF GOVERNING EQUATIONS............. 27
Vector Form of the Navier-Stokes Equations.................. 27
General Transformation
28

General Curvilinear Coordinates System......................... 31
4 NUMERICAL METHODS FOR TWO-DIMENSIONAL
PROBLEMS....................................................... 36
Discretization on a Uniform Staggered Grid..................... 37
Discretization on a Stretched Staggered Grid................... 47
Outflow Boundary Treatment..................................... 55
Consistency and Stability ..................................... 60
Consistency............................................... 60
Stability................................................. 65
Distributive Relaxation and Line-Distributive Relaxation....... 72
Distributive Relaxation................................... 72
Line-Distributive Relaxation.............................. 75
Multigrid Methods.............................................. 80
Full-Coarsening........................................... 80
Semi-Coarsening .......................................... 84
5 NUMERICAL METHODS FOR THREE-DIMENSIONAL
PROBLEMS....................................................... 88
Discretization of the Navier-Stokes Equations on a Uniform Staggered
Grid...................................................... 88
Discretization of the Navier-Stokes Equations on a Nonuniform Stag-
gered Grid..................................................... 95
vii

Analytical Mapping...................................... 95
Governing Equations under Special Mapping............... 98
Discretization..........................................100
Line-Distributive Relaxation.................................107
Multigrid Methods............................................113
Full-Coarsening.........................................114
Semi-Coarsening ........................................119
Inflow and Outflow Boundary Points .....................123
Outflow Boundary Conditions..................................125
6 METHODS FOR ANALYSIS OF NUMERICAL RESULTS 127
Perturbation Energy..........................................127
N-Factor.....................................................129
Fourier Transformation.......................................131
7 TWO DIMENSIONAL SIMULATION.................................. 135
Two-Dimensional Spatial Channel..............................135
Two-Dimensional Spatial Flat Plate...........................139
Parallel Base Flow......................................139
Non-Parallel Base Flow..................................148
8 THREE-DIMENSIONAL PLANAR CHANNEL
SIMULATION...................................................155
viii

Temporal Approach..............................................155
Linear Instability Stage..................................155
Secondary Instability.....................................162
Spatial Approach...............................................190
Linear Instability Stage..................................190
Secondary Instability................................... 191
9 THREE-DIMENSIONAL FLAT PLATE SIMULATION.......................195
Linear Instability Stage.......................................195
Parallel Base Flow .......................................195
Non-Parallel Base Flow....................................196
Secondary Instability..........................................196
10 EFFECTS OF ROUGHNESS ELEMENTS ON
TRANSITION ................................................. 202
Two-Dimensional Roughness Element for Channel Flow.............202
Isolated Roughness....................................... 202
Multiple Roughness Elements...............................204
Two-Dimensional Isolated Roughness Element for Flat Plate Boundary
Layer Flow ...............................................211
Grid Generation...........................................211
Small-Scale Roughness.....................................214
Medium-Scale Roughness....................................226
ix

Large-Scale Roughness...................................236
Two-Dimensional Distributed Multiple Roughness..........244
Three-Dimensional Isolated Roughness Element for Flat Plate Bound-
ary Layer Flow..............................................247
Base Flow...............................................248
Disturbance.............................................250
Three-Dimensional Distributed Roughness Element for Flat Plate Bound-
ary Layer Flow..............................................261
Base Flow...............................................262
Disturbance.............................................263
11 CONCLUDING REMARKS..........................................272
Summary and Conclusions.....................................272
Open Questions..............................................274
BIBLIOGRAPHY................................................276
x

LIST OF FIGURES
1.1 Secondary instability patterns: (a) K-type, (b) H-type........... 9
2.1 Neutral curve for (a) planar channel, (b) boundary layer flow. . . 22
3.1 Transformation from the physical (x,y,z) space to computational
(Â£,77,C) space................................................... 29
4.1 Staggered grid structure for 2-D channel flow. For the primary
pressure point inside the computational domain, the index is marked
from i = 2, .ni 1, j = 2, ,nj 1......... 37
4.2 Indexing scheme for the staggered grid, illustrating primary vari-
able locations for the ijih cell............................. 40
4.3 Neighbor points for the discrete xmomentum equation................. 43
4.4 Neighbor points for the discrete ymomentum equation................. 43
4.5 Neighbor points for the discrete continuity equation................. 44
4.6 Interpolation for vatau point (a) interior point, (b) point adjacent
to the solid walls............................................... 44
4.7 y-direction stretched grid....................................... 48
4.8 Staggered grid for flat plate flow in computational (Â£,77) plane. . . 50
4.9 Neighbor points for U momentum equation......................... 51
4.10 Interpolation for and at the adjacent of upper
and lower boundaries............................................... 54
4.11 Extended computational domain....................................... 55

4.12 Buffered outflow boundary points................................. 59
4.13 Neighbor points for BTCS scheme at Uij point..................... 61
4.14 Distribution of the corrections for distributive relaxation at an
interior point.................................................... 74
4.15 Distribution of the corrections for distributive relaxation at bound-
ary points: (a) inflow; (b) far field; (c) solid wall; (d) solid wall
and inflow corner................................................. 74
4.16 Distribution of corrections in the computational plane. . . 76
4.17 Two-level staggered grid (full coarsening)....................... 81
4.18 Restriction for iZu............................................ 82
4.19 Bilinear interpolation for u,v and P............................. 83
4.20 Two-level staggered grid (normal direction semi-coarsening). ... 84
4.21 Bilinear interpolation for corrections to U, V and P............. 87
5.1 3-D planar channel flow......................................... 89
5.2 Staggered grid structure.......................................... 91
5.3 Neighbor points for x momentum equation: (a) cross-section in
(x,y)-plane, z-axis orthogonal to view, (b) cross-section in (x,z)-
plane, i/-axis orthogonal to view................................. 93
5.4 Staggered grid structure in the computational (Â£,t/,Â£) space. . 101
5.5 Neighbor points for Â£momentum equation (U are at the same
points as u and are not shown here).............................105
5.6 Neighbor points of fourth-order approximation for V^_i -k and
Pi-\ik........................................................106
5.7 Extrapolation of at grid point adjacent to solid wall.............106
xu

5.8 Distribution of corrections in the (Â£,q)-plane.....................110
5.9 Projection of two-level staggered grid on the (,y)-plane (2-axis
is orthogonal to view, and the actual u,v,w and P points in the
computational domain are located in different (x,y)-planes). . 114
5.10 Full weighting restriction.........................................116
5.11 Bilinear interpolation for Au in the (jr,z)-plane................118
5.12 Trilinear interpolation for pressure points......................118
5.13 Two-level staggered grid structure for xz direction semi-coarsening:
(a) fine grid projection on (x,y) plane, (b) coarse grid projection
on (x,y) plane, (c) fine and coarse grid projection on (x,y) plane. 120
5.14 Full-weighting restriction for (a) x-momentum equation, (b) 2-
momentum equation...................................................121
5.15 Full-weighting restriction for y-momentum and continuity equations.121
5.16 Bilinear interpolation for Au......................................122
5.17 Bilinear interpolation for Aw......................................122
5.18 Bilinear interpolation for Av and Ap, (AP and Av are not located
in the same (x, 2) plane).........................................122
5.19 Modification of the two-level staggered grid structure for x z
direction semi-coarsening: (a) fine grid projection on (x,y) plane,
(b) coarse grid projection on (x,y) plane........................124
7.1 Convergence history of semi-coarsening multigrid method............136
7.2 Comparison of the numerical and theoretical solution at t=12T.
numerical solution, 0 0 linear stability solution, computational
domain: 5 T-S wavelengths physical + 1 wavelength buffer, grids:
192 x 64, e = 0.001, scheme: fourth-order finite difference. . 137
xm

7.3 Disturbance vorticity contours for the first 3 T-S wavelengths of the
physical domain; contour intervals=10-3, Re = 5000, e = 0.001,
and flow direction is from left to right. .......................138
7.4 Disturbance streamfunction contours for the first 3 T-S wave-
lengths of the physical domain; contour intervals=5 X 105, Re =
5000, e = 0.001, and flow direction is from left to right........138
7.5 Streamwise and wall-normal velocity eigenfunctions of an unstable
mode (Re% = 900) in flat plate flow..............................140
7.6 Convergence history at a fixed time for multigrid and single-grid
relaxation. .....................................................144
7.7 Comparison of the numerical and theoretical velocity components
near the solid wall (y* = 1.3137 for u and y* = 1.2448 for w).
Re* = 900, Ft = 86, e = 0.0001\/2, parallel base flow assumption
is used grids: 362 x 50 (11 T-S wavelengths physical domain + 1
wavelength buffer domain)..........................................144
7.8 Comparison of the numerical and LST velocity profiles at x* =
439.2 and 608.4. ................................................ 145
7.9 Maximum streamwise and wall-normal direction amplitudes of fun-
damental wave Ui and Vi with iZeJJ = 900, Ft = 86 and e =
V2 x 10~4..........................................................146
7.10 Test of the efficiency of the new outflow boundary treatment with
different amplitude of disturbances (Re* = 900, t = 4T, grids:
66 x 50).......................................................... 147
7.11 Amplitude of streamwise and wall-normal disturbance velocity com-
ponents of 2-D disturbance near the solid wall (y* = 1.3137 for u
and y* = 1.2448 for v). Re5 = 900, e = 0.0001-\/2. Non-parallel
Blasius base flow assumption is used, and the effect of non-parallel
base flow is shown.................................................150
xiv

7.12 Contour plots of (a) disturbance streamfunction and (b) distur-
bance vorticity at t = 13T for Re^ = 900, e = \/2 x 10-4, flow
direction is from left to right. Non-parallel base flow is used, and
only part of the computational domain is shown............... 151
7.13 Comparison of maximum streamwise amplitudes of fundamental
wave Ui, mean-flow distortion Uo, and first harmonic Ui with those
obtained by DP-81 of Joslin et al. (1992) for Req = 688.315, Ft =
86 and e = 0.0025\/2............................................152
7.14 Comparison of streamwise disturbance profiles of fundamental wave
U\ and mean-flow distortion u0 against normal distance from wall
between DP-81 of Joslin et al. (1992) and our approach (DNS+MG)
at Re* = 1519.................................................... 153
7.15 Comparison of wall-normal disturbance profiles of fundamental
wave Vi and mean-flow distortion v0 against normal distance from
wall between DP-81 of Joslin et al. (1992) and our approach
(DNS+MG) at Re* = 1519............................................154
8.1 Velocity eigenfunctions of a stable mode in plane-parallel channel
flow (Re = 5000, a = 1.0, /? = 1.005). (a) streamwise veloc-
ity eigenfunctions, (b) real part of vertical velocity eigenfunction,
(c) imaginary part of vertical velocity eigenfunction, (d) spanwise
velocity eigenfunctions...........................................157
8.2 Velocity eigenfunctions of an unstable mode in plane-parallel chan-
nel flow (Re = 7500, a = 1.0, = 0.3). (a) streamwise velocity
eigenfunctions, (b) real part of vertical velocity eigenfunction, (c)
imaginary part of vertical velocity eigenfunction, (d) spanwise ve-
locity eigenfunctions.......................................... 158
8.3 Comparison of convergence histories for relaxation only and for
relaxation used with multigrid...................................159
8.4 Perturbation energy distribution in 3 T-S periods for Re = 5000. 159
8.5 Perturbation energy distribution in 20 T-S periods for Re = 7500. 160
xv

8.6 Comparison of the numerical and theoretical velocity at t = 3T
along y = 1.555, z = 11.127 line for u and v, along y 1.555, z =
10.472 line for w. Re = 7500, e = 0.001, grid: 34 x 130 x 18. ... 161
8.7 A typical 3-D contour plot of the relative helicity for the three-
dimensional temporal channel. Re = 1500, a = 1, 0 = 1, and
the grid size is 34 x 66 x 50............................ 165
8.8 Total disturbance vorticity contours at 0.05 intervals for t = 7.99,
Re = 1500, a = 1, and 0 = 1 on a 34 X 66 x 50 grid.............166
8.9 Total disturbance vorticity contours at 0.1 intervals for t = 15.8,
Re = 1500, a = 1, and 0 = 1 on a 34 x 66 x 50 grid.............167
8.10 Total disturbance vorticity contours at 0.3 intervals for t = 23.6,
Re = 1500, a = 1, and 0 = 1 on a 34 x 66 x 50 grid.............168
8.11 Total disturbance vorticity contours at 0.4 intervals for t = 27.56,
Re = 1500, a = 1, and /3 = 1 on a 34 x 66 x 50 grid............169
8.12 Total disturbance vorticity contours for t = 29.1, Re = 1500,
a = 1, and /3 = 1 on a 34 x 66 x 50 grid, contours from 0.05 to 45.170
8.13 Streamwise disturbance vorticity contours at 0.20 intervals for t
7.99, Re = 1500, a = 1, and (3 = 1 on a 34 X 66 X 50 grid;
dashed lines indicate negative contours....................171
8.14 Streamwise disturbance vorticity contours at 0.20 intervals for t =
15.8, Re = 1500, a = 1, and (3 = 1 on a 34 x 66 x 50 grid; dashed
lines indicate negative contours...................................172
8.15 Streamwise disturbance vorticity contours at 0.40 intervals for t =
23.6, Re 1500, a = 1, and 0 = 1 on a 34 x 66 x 50 grid; dashed
lines indicate negative contours.................................. 173
8.16 Streamwise disturbance vorticity contours at 0.50 intervals for t =
27.56, Re = 1500, a = 1, and 0 = 1 on a 34 X 66 x 50 grid;
dashed lines indicate negative contours............................174
xvi

8.17 Streamwise disturbance vorticity contours at 0.50 intervals for t =
29.1, Re = 1500, a = 1, and /3 = 1 on a 34 x 66 x 50 grid; dashed
lines indicate negative contours...............................175
8.18 Contour plots of du/dy at t = 0.15 in the (*,y)-plane. Re =
1500, a = 1, and (3 = 1. The contour interval is 0.2, flow
direction is from left to right................................176
8.19 Contour plots of du/dy: (a) one-spike stage, experimental results
by Nishioka et al.(1981); (b) computations at t = 7.99 in the (sea-
plane, contour interval is 0.2. Flow direction is from left to right. 176
8.20 Contour plots of du/dy: (a) three-spike stage, experimental results
by Nishioka et al.(1981); (b) computations at t = 15.8 in the (sea-
plane, contour interval is 0.2. Flow direction is from left to right. 177
8.21 Contour plots of du/dy: (a) five-spike stage, experimental results
by Nishioka et al.(1981); (b) computations at t = 23.6 in the (x,y)-
plane, contour interval is 0.2. Flow direction is from left to right. 178
8.22 Contour plots of du/dy at: (a) t = 27.56; (b) t = 29.1. Flow
direction is from left to right. Results show that the high-shear
layer has broken and spread to the whole domain, and the scale of
the flow has become smaller....................................179
8.23 Total disturbance vorticity contours at 0.025 intervals for t = 23.8,
Re = 7500, a = 1, and /3 = 0.5 on a 50 x 66 x 98 grid.........181
8.24 Total disturbance vorticity contours at 0.05 intervals for t = 47.3,
Re = 7500, a = 1, and /3 = 0.5 on a 50 x 66 x 98 grid.........182
8.25 Total disturbance vorticity contours for t = 47.3, Re = 7500,
a = 1, and = 0.5 on a 50 x 66 x 98 grid, contours from 0-180
(nonuniform contour intervals are used)..................183
8.26 Streamwise disturbance vorticity contours at 0.20 intervals for t
23.8, Re = 7500, a = 1, and fi = 0.5 on a 50 x 66 x 98 grid;
dashed lines indicate negative contours...........................184
xvn

8.27 Streamwise disturbance vorticity contours at 0.20 intervals for t =
47.3, Re = 7500, a = 1, and /3 = 0.5 on a 50 x 66 x 98 grid;
dashed lines indicate negative contours....................185
8.28 Streamwise disturbance vorticity contours at 0.40 intervals for t =
59.52, Re = 7500, a. = 1, and /3 = 0.5 on a 50 x 66 x 98 grid;
dashed lines indicate negative contours....................186
8.29 Contour plots of du/dy at t = 23.8, Re = 7500, a = 1, jd = 0.5.
Flow direction is from left to right..............................187
8.30 Contour plots of du/dy at t = 47.3, Re = 7500, a = 1, /? = 0.5.
Flow direction is from left to right..............................188
8.31 Contour plots of du/dy at t = 59.52, Re = 7500, a = 1, = 0.5.
Flow direction is from left to right..............................189
8.32 Comparison of the numerical (^>) and theoretical () solutions
at selected xlines at t = 7T on a 182 x 66 X 18 grid (5 T-S
wavelengths physical domain 1 wavelength buffer), with Re =
7500, a; = 0.25, e = 10-4, and (3 = 0.3. Flow direction is from left
to right............................................................192
8.33 Contour plots of du/dy at 0.30 interval for selected (s,y)-planes
on a 162 x 66 x 50 grid at t = 91.23 (7 T-S wavelengths physical
domain 1 wavelength buffer domain), with Re = 7500,u =
0.25,/d = 0.3, 62d = 0.03, 63,* = 0.01. Flow direction is from left
to right............................................................193
8.34 Total disturbance vorticity contours at 0.05 intervals for selected
(*, z)-planes on a 162 x 66 x 50 grid (7 T-S wavelengths physical
domain + 1 T-S wavelength buffer domain), with Re = 7500,
u; = 0.25, (3 = 0.3, and t = 91.23. Flow direction is from left to
right...............................................................194
xvui

9.1 Comparison of the numerical and theoretical velocity components
near the solid wall (yjj = 1.132 for u and w and yJJ = 1.027 for u).
Re^ = 900, e = 10~3, parallel wall-bounded base flow assumption
is used........................................................198
9.2 Streamwise, wall-normal, and spanwise disturbance velocity com-
ponents of a 3-D disturbance near the solid wall (yj = 1.132 for u
and w and yj = 1.0247 for v). Re\J = 900, e = 103, non-parallel
base flow is used, and the effect of non-parallel base flow is shown. 199
9.3 Spanwise vorticity contours on the central (*, y)plane of the com-
putational domain at (a) t=244.17, and (b) t=270.97................. 200
9.4 High-shear layer at one spike stage of transition in boundary layer
flow, (a) Numerical result in a selected part of the central (x,y)-
plane, (b) experimental results obtained by Kovasznay et al.(1962). 200
9.5 Total vorticity contours on the yj = 0 (x,z)plane at (a) t =
244.17, (b) t = 270.97; Re% = 900, e2d = 0.03, e3d = 0.01,
grids: 202 x 50 x 34, and contour interval is 0.035. Flow direction
is from left to right. ..........................................201
10.1 Grids with one isolated roughness for the channel flow............205
10.2 Contour plots of the streamfunctions of the base flow with one
roughness element. Re = 5000, k = 0.15. Flow direction is from
left to right...................................................205
10.3 Instantaneous contour plots of the perturbation streamfunctions
for the one roughness element case. Re = 5000, k = 0.15, e =
0.0025v/2, u> = 0.33017. Flow direction is from left to right. . 205
10.4 Maximum amplitudes of fundamental wave u\, Vi, mean-flow dis-
tortion u0, v0, and first harmonic u2, v2 for Re = 5000, / =
0.15, w = 0.33017, and e = 0.0025-\/2 with one roughness element
(grid: 170 x 50 x 4)
206

10.5 Contour plots of the perturbation streamfunctions and vorticities
for Re = 5000, k = 0.12 u = 0.33017, and e = 0.0025\/2 with two
roughness elements (grid: 402 X 66 x 4)......................... 207
10.6 Contours of total streamfunctions and vorticities for Re = 5000, K\ =
0.12, u) = 0.33017, and e = 0.0025\/2 with two roughness elements
(grid: 402 x 66 x 4)............................................ 207
10.7 Maximum amplitudes of fundamental wave U\, V\, mean-flow dis-
tortion Uq, Vo, first harmonic wave u2, v2, and second harmonic
wave ii3, u3 for Re = 5000, k -- 0.12, u> = 0.33017 and e =
0.0025\/2 with two roughness elements (grid: 402 X 66 x 4). ... 208
10.8 Contour plots of the instantaneous perturbation streamfunctions
and vorticities for Re = 5000, k = 0.12, u> = 0.33017, and e =
0.0025\/2 with seven roughness elements (grid: 402 x 66 x 4). . 209
10.9 Contour plots of total streamfunctions and vorticities for Re =
5000, k = 0.12, u> = 0.33017, and e = 0.0025'\/2 with seven rough-
ness elements (grid: 402 x 66 x 4).............................. 209
10.10 Maximum amplitudes of fundamental wave u\, Vi, mean-flow
distortion uq, Vo, first harmonic wave u2, v2, and second har-
monic wave -u3, v3 for Re = 5000, it = 0.12, u) = 0.33017, and
e = 0.0025\/2 with seven roughness elements (grid: 402 X 66 X 4). 210
10.11 Grid generation..............................................212
10.12 A typical grid used for simulation of rough flat plate boundary
flow (part of the grid is shown) ..............................213
10.13 Computational domain for the 2-D isolated roughness cases. . 215
10.14 Comparison of the normal direction distribution of the streamwise
velocity component of the base flow with small-scale roughness case
I and the Blasius profiles at (a) x = 20, (b) x = 30, (c) x = 40,
(d) x = 60. Roughness element located at x = 35.41, k = 0.27,
Re* = 500....................................................... 218
xx

215
10.15 Streamwise distribution of 6*, 6, and H for small-scale roughness
case I.........................................................
10.16 Comparison of the amplification factor of disturbance for small-
scale roughness case I and smooth flat plate. ReJ = 500, k = 0.27,
u> = 0.07, e = 0.003..................................... 219
10.17 Maximum streamwise and normal amplitudes of mean flow dis-
tortion iio, Vo, fundamental wave ulf V\, and first harmonic wave
v>2, V2, with nondimensional downstream distance for small-scale
roughness case I and smooth flat plate. Re$= 500, u) = 0.07, e = 0.003.............................................. 220 10.18 Contour plots of steady base flow streamfunctions and instan- taneous perturbation streamfunctions at t = 9T for small-scale roughness case I. Flow direction is from left to right....221 10.19 Comparison of the normal direction distribution of the streamwise velocity component of the base flow with small-scale roughness case 11 and the Blasius profiles at (a) x = 20, (b) x = 30, (c) x = 41, (d) x = 65. Roughness element located at x = 35.41, k = 0.27, Re; = 500....................................................... 222 10.20 Streamwise distribution of 6*, 9, and H of the base flow for small- scale roughness case II..........................................223 10.21 Comparison of the amplification factors of disturbance for small- scale roughness case II and smooth flat plate. Re^ = 500, w = 0.07, e = 0.003....................................................... 223 10.22 Maximum streamwise and normal direction amplitudes for mean flow distortion Uq, Vq, fundamental wave U\,vi, and first harmonic wave v>2, V2, with nondimensional downstream distance for small- scale roughness case II and smooth flat plate. Re^ = 500, w = 0.07, e = 0.003....................................................... 224 xxi 10.23 Contour plots of steady base flow streamfunctions and instan- taneous perturbation streamfunctions at t = 9T for small-scale roughness case II. Flow direction is from left to right.......225 10.24 Streamfunction contours of the base flow for medium-scale rough- ness case I. Flow direction is from left to right.................228 10.25 Comparison of normal direction distribution of the streamwise velocity component of the base flow with medium scale roughness case I and Blasius profiles at (a) x = 20, (b) x = 40, (c) x = 52, (d) x = 85. Roughness element location at x = 44.0, k = 0.5, Re; = 663.12.................................................... 229 10.26 Streamwise distribution of 6*, 0, and H for medium-scale rough- ness case 1..........................................................230 10.27 Comparison of the amplification factor of disturbance for medium- scale roughness case I. Re; = 663.12, w = 0.0968, e = 0.004. . . 230 10.28 Maximum streamwise and normal amplitudes of mean flow distor- tion u0, Vq, fundamental wave Hi, Vi, and first harmonic wave u2, v2, with nondimensional downstream distance for medium-scale roughness case I and smooth flat plate. Re; = 663.12, u> = 0.0968, e = 0.004...................................................... 231 10.29 Contour plots of instantaneous perturbation streamfunctions at t = 9T for medium-scale roughness case I. Flow direction is from left to right..................................................232 10.30 Contour plots of steady base flow streamfunction for medium-scale roughness case II. Flow direction is from left to right........232 10.31 Comparison of the normal direction distribution of the streamwise velocity component of the base flow with medium-scale roughness case I and the Blasius profiles at (a) x = 20, (6) x = 40, (c) x = 52, (d) x = 85. Roughness element located at x = 44.0, k = 0.50, Re; = 663.12................................... 233 xxn 10.32 Streamwise distribution of Â£*,0, and H for medium-scale rough- ness case II...........................................................234 10.33 Comparison of the amplification factor of disturbance for medium- scale roughness case II and smooth flat plate. ReJ = 663.12, = 0.0968, e = 0.004................................................ 234 10.34 Maximum streamwise and normal amplitudes of mean flow dis- tortion uq, Vo, fundamental wave Vi, and first harmonic wave u2,v2, with nondimensional downstream distance for medium-scale roughness case II and smooth flat plate. Re = 663.12, ui = 0.0968, e = 0.004...................................................... 235 10.35 Contour plots of instantaneous perturbation streamfunctions at t = 9T for medium-scale roughness case II. Flow direction is from left to right.....................................................236 10.36 Streamfunction contours of the base flow for large-scale roughness case I. Flow direction is from left to right....................240 10.37 Comparison of the numerical streamwise velocity profile for large- scale roughness case I and experimental results obtained by Dovgal and Kozlov (1990) (over a square hump)..........................240 10.38 Comparison of the numerical streamwise disturbance of shape factor H and the experimental results obtained by Dovgal and Kozlov (1990) (over a square hump)................................240 10.39 Streamwise distribution of the amplification factors for large-scale roughness element case 1..........................................241 10.40 Normalized amplitude of streamwise perturbation velocity for large-scale roughness case 1......................................241 10.41 Contour plots of the base flow for large-scale roughness case II. Flow direction is from left to right............................242 XXUl 10.42 Comparison of numerical shape factor distribution for large-scale roughness case II and experimental results obtained by Dovgal and Kozlov (1990) for a square hump..............................242 10.43 Comparison of streamwise velocity profile at selected locations with experimental results. Re% = 663.12, k = 1.12................242 10.44 Streamwise distribution of amplification factors for large-scale roughness case II.................................................243 10.45 Streamwise distribution of the normalized streamwise disturbance velocity profiles for large-scale roughness case II..............243 10.46 Contours plots of instantaneous perturbation streamfunctions at t 9T, Re*Q = 663.12, to = 0.0968, e = 0.003. Nonuniform contour intervals are used, and flow direction is from left to right.243 10.47 Streamfunction contours of the base flow for multiple distributed roughness elements case (part of the domain is shown)............245 10.48 Nondimensional shape parameter H, displacement thickness 6* and momentum thickness 6 of the base flow with the presence of multiple roughness elements.......................................245 10.49 Streamwise distribution of amplification factors; (a) amplification factors for u, (b) amplification factors for v...................246 10.50 Contour plots of the instantaneous perturbation streamfunctions at t = 10T for the distributed multiple roughness elements (part of the domain is shown, and nonuniform contour intervals are used). Flow direction is from left to right..............................246 10.51 Vector plots of the base flow velocity on (a) j 2 and (b) j 3 grid surface for the small-scale 3-D roughness case. Re$ = 663.12,
ReK = 31.83, k = 0.3. Flow direction is from left to right.......253
xxiv

10.52 Vector plots of the base flow velocity on j = 2 and j = 3 grid
surfaces (from bottom to top) for the medium-scale 3-D roughness
case. Re J = 663.12, .Re* = 123.3, #e = 0.6. Flow direction is from
left to right.......................................................254
10.53 Contour plots of the streamwise vorticity on i = 50,55,60, 65,70,75
(y, z)-planes for the large-scale 3-D roughness case. TZeJ = 663.12,
ReK = 418.14, k = 1.12. Flow direction is from left to right. . 255
10.54 Vector plots of the base flow velocity on (a) j = 3 grid surface,
(b) on the y = 0.245 plane, for the large-scale 3-D roughness case.
Re*Q = 663.12, ReK = 418.14, k = 1.12...............................256
10.55 Contour plots of the base flow vorticity magnitude on the y =
0.245 (, z)-plane and j = 3 grid surface (from top to bottom)
for the large-scale 3-D roughness case. Re% = 663.12, ReK =
418.14, k = 1.12. Flow direction is from left to right............257
10.56 Contour plots of the disturbance vorticity magnitude on j = 2
and j = 3 grid surface (from bottom to top) for the small-scale
3-D roughness case at t = 5T. ReJ = 663.12, k = 0.3, e = 0.003.
Flow direction is from left to right...............................258
10.57(a) Comparison of amplification factors of u between the 2-D smooth
flat plate and 3-D small-scale isolated roughness case, (b) Com-
parison of the disturbance velocity magnitude amplification factors
between 2-D smooth flat plate and 3-D small-scale isolated rough-
ness case, (c) Amplification factors of spanwise velocity magnitude
w for the 3-D small roughness case. iZeJ = 663.12, k = 0.3, e =
0.003............................................................. 259
10.58 Contour plots of the perturbation vorticity on j = 2 and j =
3 (from bottom to top) grid surfaces, for the medium-scale 3-D
roughness case. Re5 = 663.12, k = 0.6, e = 0.003. Flow direction
is from left to right..............................................260
10.59 Distribution of the 3-D distributed roughness surface...........265
XXV

10.60(a) Vector plots of velocity fields on j = 2 and j = 3 grid surfaces
for the small-scale distributed roughness case....................265
10.60(b) Contour plots of streamwise vorticity on i = 30,42, 54,66,78,90
(y, z)-planes for the small-scale distributed roughness case. Flow
direction is from bottom-left to top-right..............................266
10.61(a) Vector plots of velocity fields on j = 2 and j = 3 grid surfaces
for the large-scale distributed roughness case....................267
10.61(b) Contour plots of streamwise vorticity oni = 30,42, 54,66,78,90
(y, z)-planes for the large-scale distributed roughness case. Flow
direction is from bottom-left to top-right..............................268
10.62 Distribution of uvelocity profiles of the base flow for the large-
scale 3-D distributed roughness case..............................269
10.63 Amplification factor of the maximum perturbation velocity for the
small-scale distributed roughness...................................269
10.64 Amplification factor of the perturbation u-velocity for the small-
scale distributed roughness.........................................270
10.65 Contours of spanwise perturbation vorticity on j = 3, 5 and 7
grid surfaces (from bottom to top) for the small-scale distributed
roughness. e2d = 0.005, k = 0.35, ReK = 43.5. Flow direction is
from left to right..................................................270
10.66 Contour plots of relative helicity of the total flow at t = 5T for
the small-scale distributed roughness. e2j = 0.015, k = 0.35, ReK =
43.5. Flow direction is from left to right........................271
10.67 Contour plots of relative helicity of the total flow at t = 5T for
the small-scale distributed roughness. e2d = 0.05, k = 0.35, ReK =
43.5. Flow direction is from left to right..........................271
xxvi

LIST OF TABLES
7.1 Relative Â£2 error norm for u and v after 13 T-S periods (calculated
in physical domain)...........................................143
8.1 Relative error Â£2-norms for the z = 0.654 plane for u and v, and
for the z = 0 plane for w.....................................160
8.2 Relative Â£2 norms of the numerical solutions at t = 7T {Re = 7500,
and a 30 x 66 X 18 grid is used for each wavelength)..........191
10.1 Flow parameters for 2-D roughness flow.......................214
10.2 Computational parameters for 2-D roughness flow..............214
10.3 Flow parameters for 3-D roughness flow.......................248
10.4 Computational parameters for 2-D roughness flow..............248

ACKNOWLEDGEMENTS
I would like to take this opportunity to thank my advisor, Professor Steve
McCormick, for the guidance and encouragement that he provided during all
phases of the preparation of this thesis. Throughout my Master and Doctoral
studies, he provided exceptional insight and guidance.
My sincerest thanks also go to Professor Chaoqun Liu, who provided me with
close supervision and who guided me step by step through algorithm development.
His outstanding knowledge of computational fluid mechanics and his continuous
support, encouragement and interest greatly enhanced the quality of this work. I
significantly benefited from the countless and fruitful discussions that I had with
him.
Thanks also go to the other committee members: Dr. John Ruge, Dr. Tom
Russell, and Dr. Randall Tagg, for their contribution to my development.
Special thanks and appreciation go to Dr. Tom Zang and Dr. Ron Joslin,
at NASA Langley Research Center, for their many contributions, technical guid-
ance, support and encouragement. Special thanks also go to Dr. Helen Reed, at
Arizona State University, and Pat Roache, at Ecodynamics Research Associates,
I would like to thank Mr. Phil Smith, the computer system administrator,
and Mr. Zhixing Liu and Mr. Wenqian Shi, for their help in the final stage of

finishing this thesis. I would also like to thank Debbie Wangerin and Stephanie
Sanchez for the help they offered.
I am also indebted to NASA for their support of this work under Research
Contract NAS1-19312.
XXIX

CHAPTER 1
INTRODUCTION
General Overview
The problem of explaining the origin of turbulent flow has been recognized for
more than a century. To. date, because of the synergistic combination of novel the-
ory with large-scale numerical computation, remarkable developments has been
achieved. But we still cannot fully understand the real origin of transition, even
the most simplest cases: planar channel and flat plate flow.
Historically, evidence has been gathered in cumulative fashion, with the bulk
of opinion about transition swaying dramatically from one misconception to an-
other. The complexities of the phenomena continued to impede a definitive pic-
ture of the transition process for quite a long period (cf. White 1974). Fortu-
nately, after a century of study, some aspects of transition are becoming well
understood. Furthermore, control of the instabilities leading to transition and fi-
nal turbulence, as well as prediction of the transition location, are now of practical
interest.
Numerical methods have been applied to theoretical studies of instability and
transition from laminar to turbulent flow since the modern digital computers
appeared. Although many approximate theoretical methods have been devel-

oped for the linear (primary) instability (cf. Orr 1907a,b; Sommerfeld 1908) and
nonlinear interaction (cf. Smith & Bodonyi 1982; Herbert 1983a,b,c,d; Pugh &
Saifman 1988) stages, as well as secondary instabilities arising from the inter-
action of three-dimensional disturbances with primary instability wave (Herbert
1988), these theories are all highly dependent on the quality and restrictions of
the numerical methods being used.
Early stages of transition starts with the appearance of a 2-D, infinitesimile
disturbance. In this stage, the Orr-Sommerfeld equation, which was obtained
from a parallel base flow assumption on the linearized perturbation form of the
Navier-Stokes equations, gives quite good prediction, and the solution of this
equation depicts the behavior of the travelling wave the Tollmien-Schlichting
(T-S) wave. Unfortunately, when the amplitude of this T-S wave reaches a fi-
nite value of about 1% of the base flow maximum velocity values, nonlinearity
becomes nonnegligible, so, the behavior of the disturbance no longer obeys the
linear theory. After the occurrence of such a nonlinearity, the three-dimensional
disturbances soon appear as the T-S waves rather quickly begin to show spanwise
variation. This tendency toward rapid achievement of three-dimensionality was
first investigated by Klebanoff & Tidstrom (1959), and KlebanofF, Tidstrom &
Sargent (1962), who show that the shear layer in the unstable region has a strong
ability to amplify any slight three-dimensional disturbance that surely must be
present in any natural-disturbance spectrum. Also, aligned three-dimensional
structures that are associated with the peak-valley splitting were investigated.
2

This type of secondary instability is named after Klebanoff and called a K-type
pattern (see Figure 1.1 (a)). The three-dimensional amplification was verified by
Benney & Lin (1960), Benjamin (1962), and Criminale & Kovasznay (1962). In
the analysis of Benney & Lin, they considered the interaction between a two-
dimensional T-S wave and a pair of oblique three-dimensional waves, and inter-
esting results are obtained by applying this idea to the inflow boundary condition
for numerical simulation. Later experiments on boundary layers performed by
Kachanov, Kozlov h Levchenko (1978) formed another type of secondary insta-
bility, which leads to the subharmonic breakdown of the laminar flow; this type
of instability was later named after Herbert and called an H-type pattern (see
Figure 1.1 (b)).
The first detailed experimental investigation of planar channel transition was
performed by Nishioka, Iida & Ichikawa (1975) and Nishioka, Asai h Iida (1980).
Only fundamental breakdown was investigated because of the special inflow con-
ditions, though the existence of subharmonic instability was also investigated
later by Ramazanov (1985).
Instability with the appearance of roughness is quite a new topic in the field
of transition research. This difficult largely neglected research area is primarily
experimental, although numerical simulations can help to clarify some features of
the mechanisms (cf. Morkovin 1990). To date, no specific instability mechanisms
have been identified experimentally or theoretically in this field. The irregularity
of the surface was believed to change the phase of perturbation waves, and, on the
3

other hand, the inevitable local separation zones are thought to lead to higher
amplification rate three-dimensional disturbances, which may develop directly
into turbulent flow by the vorticity stretching instability (Corke, Bar-Sever, h
Morkovin 1986).
Transition from laminar to turbulent flow in flat plate boundary layers and
in planar channel flows has been considered as the fundamental subject in the
area of numerical simulation because of its simple geometry and exact base flow.
Due to the lack of computer facilities, early research works were restricted to
temporal approaches. In this kind of approach, streamwise flow periodicity is
assumed. Temporal simulations follow the time evolution of a single wavelength
of the disturbance, so this kind of simplified simulation can obtain only upstream
effects only within the limitation of the periodicity. Furthermore, nonparallel
base flow is not allowed in this approach, although something like changing the
base flow in time is possible. Although this approach has obvious limitations, it
did achieve some notable successes, advancing recently to the accurate simulation
of complete transition to turbulence.
Spatial simulation, which consider the flow transition in a more realistic do-
main, became possible after the appearance of the modern supercomputers. In
this approach, a finite streamwise length is assumed, so that it is comparable to
the experiments in a wind tunnel. One of the most productive groups in spatial
simulation is Fasel et al. (cf. Fasel 1976; Fasel & Bestek 1980; Fasel, Rist &
Konzelmann 1990; Fasel & Konzelmann 1990). Later works performed by Dan-
4

abasoglu, Biringen & Streett (1989,1991), Liu, Liu & McCormick (1991), Liu
& Liu (1992a), Liu & Liu (1992b), and Joslin, Streett Sz Chang (1992) further
improved this approach to more practical applications.
Spatial simulation of 2-D boundary layer transition was first accomplished by
Fasel (1976) using a second-order finite difference method in both stream and
wall-normal directions. The governing equations he used are in the so-called
vorticity-velocity form, and the scheme is fully-implicit. In later work, Fasel
improved his method to fourth-order finite differences in space (Fasel & Bestek
1980) to allow for investigation of 2-D channel transition, and, later to 3-D flat
plate boundary layer transition ( Fasel, Rist &z Konzelmann 1987, 1990; Fasel &
Konzelmann 1990).
It is well-known that how to treat the outflow boundary conditions down-
stream of the computational domain is one of the main difficult sources of spatial
simulation. In CFD, these outflow boundary conditions for incompressible flows
are treated using extrapolation conditions from upstream, but these extrapola-
tion conditions are not adequate for transition simulation because of the high
sensitivity of these problems. In Fasels work, the streamwise wavenumber was
incorporated to specify the outflow boundary conditions. This approach worked
well at small disturbance stage, but caused reflection for large amplitudes. An-
other important work is that accomplished by Patera (1984). In his work, a
2-D planar channel is considered, and a modified Fourier spectral method cou-
pled with explicit time marching was developed. However, outflow boundary
5

conditions based on vanishing gradients along the streamwise direction caused
considerable errors. More recently, Streett & Macaraeg (1989) developed and
implemented a very efficient buffer domain method in their fully spectral simula-
tions of 2-D planar channel flow. Also, another method applied by Spalart (1989,
1990) uses periodic boundary conditions in the x-direction, but force functions
are used to suppress waves in the fringes of the streamwise domain (cf. Kleiser $z Zang 1991). A new method that used the buffer domain method with an increase in downstream viscosity was developed by Liu & Liu (1992a) and Liu &; Liu (1992b). Because of the specialty of the staggered grid they used, this method is extremely efficient for both small and large disturbances. Only two computational studies exist concerning the spatial stability of a sep- arated flow: Fasel, Bestek h Schefenacker (1977), and Danabasoglu, Biringen & Streett (1993). Only the 2-D case was considered in these studies. Qualitatively, these works are comparable to what is observed in the physics, but there is still a long way to go for satisfactory quantitative comparisons. Advances Represented by This Study The current status of transition simulation shows that, even though a lot of methods have been developed, including the highly accurate spectral methods, the efficiency of the codes developed are not as good as could be expected. Also, because of various limitations of the spectral methods, geometries are restricted to rectangular cases, like planar channel, flat plate, or steps, which are far from 6 the geometries of practical interest. The motivation of this work is triggered by demanding a very accurate and very efficient multigrid code that can simulate transition on a 3-D rough surface. Since current numerical studies are still rather limited, this motivation demands new numerical approaches. The rather new approach used in this work includes: General coordinates that are able to work for more general geometry Multilevel stretched and staggered grids that are generally accepted to have higher accuracy than that of non-staggered grids Fully implicit time-marching (implicit for both convection and diffusion) that has much better stability than that of explicit methods Fourth order discretization in all three directions that has high accuracy and can reduce the artificial viscosity and phase error Line distributive relaxation that is very efficient smoother for Navier-Stokes equations Semi-coarsening multigrid scheme for solving the large scale nonlinear al- gebraic systems arising in the implicit approach at each time step, which combines with line distributive relaxation to make fully implicit schemes comparable to explicit schemes at each time step in CPU cost New out-flow boundary treatment that efficiently eliminates all wave reflec- tions from the outflow boundary, and reduces the size of buffer domain to 7 as short as at most one T-S wavelength for both small and large magnitude of inflow disturbance, which greatly enhances the efficiency of spatial DNS Since the new approach developed in this work has high accuracy, stability and efficiency, not only some physics in transition such as spatially linear growth and secondary instability in 3-D channels and flat plates have been successfully simulated in an efficient way, but also some new numerical simulation regimes have been reached: we have successfully simulated transition over multiple 2-D and 3-D roughness elements, with some general shape and general distribution in 2-D and 3-D channels and flat plates. 8 (b) Figure 1.1. Secondary instability patterns: (a) K-type, (b) H-type. 9 CHAPTER 2 GOVERNING EQUATIONS IN PHYSICAL SPACE In this chapter, the governing equations in Cartesian coordinates are pre- sented in their conventional form in physical space. The details of the non- dimensionalization are given, and the proper boundary conditions are specified. Since both the Blasius similarity solution and the linear stability theory are fre- quently used through out this report, for completeness, they are also included in this chapter. Original Form Low speed flows can be governed by the incompressible Navier-Stokes equa- tions, in which momentum and mass conservation are represented by three mo- mentum equations and a continuity equation, respectively. Considering only Cartesian coordinates, the governing equations can be expressed as follows: x-momentum equation Du dP , . PDi = + (:U) y-momentum equation z-momentum equation Dw dP P~F>7 = K~ + fiAw + fzt Dt dz (2.3) continuity equation Dp .du dv dw. Di + p{di + di+ a7) = 0 (2.4) where, u, v, w are dimensional streamwise, normal, and spanwise velocity com- ponents, respectively, P is the dimensional pressure, p represents the density of the fluid, /*, fy, and fz are body forces in the three directions, p is the dynamic viscosity coefficient, A is the three-dimensional Laplacian, and ^ denotes the material derivative. In Cartesian coordinates, we have A = d2 d2 d2 dx2 dy2 ^ dz2 (2.5) D d d d d (2.6) Dt dt^ dx dy dz It is assumed that no density variation exists in the investigation domain, thus, Â£ = 0, Dt (2.7) and the continuity equation becomes du dv dw dx dy dz (2.8) Furthermore, assuming the effect of the body force / = (/,/y,/z) is negligi- ble, the momentum equations can be written as du du du du 1 dP .d2u d2u d2u . dt+Ud^ + Vd^ + Wd^ + ^~^d^+dy^ + d^^t * 11 (2.10) dv dv dv dv 1 dP d2v d2v d2v. m+uai + vai + wrz + -pfy-'iw + v + a?)= dw dw dw dw 1 dP . d2w d2w d2w. , . l)i + Ulte+Vlfy+Wlh+plh~V(lW+lhj*+lh*^ = 0 ^2'U^ Here, v = ^ is the kinematic viscosity coefficient. Nondimensionalized Forms for Planar-Channel-Type Flow Equations (2.8)-(2.11) can be nondimensionalized by dividing them by some reference constants appropriate to the flow, for instance, reference velocity Vo and reference length L. For planar-channel-type flow, we use the channel half height h as the reference length and the central line velocity Uo of the steady flow as the reference velocity, i.e., L = h, F0 = U0. (2.12) Using asterisks to denote the resulting dimensionless variables, we have U * w u = Uo' fe II ?> w = Uo X * y . z tj B if H h' p pUS y =h tu0 h 2 = By defining the Reynolds number v 12 (2.13) > (2.14) we finally obtain the following nondimensionalized governing equations (without asterisk here and henceforth): du du du du dP 1 ,d2u d2u d2u. dt ^ 5 dy dz dx Re'dx2^ dy2 dz2 dv dv dv dv dP 1 ,d2v d2v d2v. ----(- u--v-----(- w---1--------(-----1-----1----1 = 0, dt dx dy dz dy Re dx2 dy2 dz2 dw dw dw dw dP 1 ,d2w d2w d2w 4- U-----h V------h W------1------------( dt ^ dx^ dy dz^ dz Re'' dx2 + dy2 + dz2 ) = 0, du dv dw dx ^ dy ^ dz (2.15) (2.16) (2.17) (2.18) It is recommended to use the conservative form of the momentum equation in the numerical simulation, which can be derived from (2.15)-(2.18): du duu duv duw dP 1 ,d2u d2u d2u. ~dt + Ih + +1)7 + + dy2^ = ^ dv duv dw dvw dP 1 ,d2v d2v d2v. . ... ~dt +lh+~fy+lh~ + lfy~te('dx2 + d^ + Ih2^ ^ ^ dw duw dvw dww dP 1 ,d2w d2w d2w. , +-97 + 3r+a7&(^+aF+8?) = 0- (2'21) The major purpose of this work is to investigate the behavior of disturbances, so the perturbation form of the equations plays a more important role. The perturbation equations are obtained by decomposing the total flow into steady base flow and a fluctuating unsteady perturbation. Letting the subscript 0 denote the base flow variables and prime the fluctuating variables, we have u(x,y,z,t) = u0(x,y,z) + u'(x,y,z,t), v(x,y,z,t) = v0(x,y,z) + v\x,y,z,t), w(x,y,z,t) = w0(x,y,z) + w'(x,y,z,t), P(x,y,z,t) = P0{x,y,z) + P'{x,y,z,t). 13 (2.22) Substituting (2.22) into (2.18)-(2.21), and noting that the base flow itself also satisfies the Navier-Stokes equations, we have x-momentum equation du' d(2u0u + u'u') d(uQv' + u'vQ + u'v') dt dx dy d(uow' + u'w0 + u'w') dz + dP' 1 d2u' d2u' d2u' 2 a.,7. Qz2 / dx Re' dx2 dy2 (2.23) y-momentum equation dv' d{u'v o + Uqv1 + u'v1) d(2v0v' + v'v') dt ^ dx ^ dy ^ d(v ow' + v'w0 + v'w') dz + dp' 13V d2v' d2v'_ dy Re' dx2 + dy2 + dz2^ ~ (2.24) z-momentum equation dw' ^ d(u0w' + u'w0 + u'w1) ^ d(v'wo + v0w' + v'w') dt dx d(2w'wo + w'w') dy dz + dP' 1 .d2w' d2w' d2w' dz Re dx2 ^ dy2 dz2 (2.25) continuity equation: du' dv' dw' dx ^ dy ^ dz (2.26) For simplicity, we will omit all primes henceforth, so that all variables without any sub/superscript will denote the fluctuating part of the flow in the perturbation 14 equations. In planar channel flow, since both Vq and Wq vanish and Uq = Uo(y), the above momentum equations can be simplified to du d(2uoU + uu) duv duw dP dt ^ dx dy dz dx 1 d2u d2u d2u _ Re dx2 dy2 dz2 dv d(u0v + uv) dvv dvw dP dt dx dy dz dy 1 'd2v d2v d2v _ ^d^2 + d^2 + dl2^ dw d(u ow + uw) dvw dww dP dt ^ dx ^ dy ^ dz ^ dz 1 .d2w d2w d2w ^'dx2+ ~dy2+~d^2^ ' (2.27) (2.28) (2.29) Nondimensionalized Forms for Flat-Plate-Type Flow There are several methods to nondimensionahze the Navier-Stokes equations for flat-plate-type flow. Since the choice of reference length for this kind of flow seems to be somewhat arbitrary, we can for instance choose the distance from the leading edge to the inflow boundary, Lq, or the thickness of the boundary layer at inflow, Â£0, or the thickness of the displacement at that location, or something else, as the reference length. In this paper, we choose 6Â£ as the reference length, and the free stream velocity Ugo as the reference velocity. In this setting, the Reynolds number is defined as Re*0 = (2.30) The general form of the governing equations is the same as (2.15)-(2.18) if we 15 substitute Re by Re5. For the perturbation flow, if the flat plate is used, the base flow is two-dimensional, and w0 vanished from equations (2.23)-(2.25). Thus, the resulting governing equations are du d(uu + 2'Uo'u) d(uv + uqv + uv0) d(uw + u0w) ~dt + Tx + Fy ~~ + dz dP 1 .d2u d2u d2u _ ~di ~ lie{W + dtf + dz*> ~ dv d(uv + u0v + uv0) d(vv -f 2v0v) dt dx dy dP 1 . d2v d2v d2v ~dy ~ R^lh? + dtf + ~d? = dw d(uw + uow) d(vw + vow) dt ^ dx dy dP 1 .d2w d2w d2w !h ~ R^~d^ + Hy2+ = du dv dw + d[yw + v0w) + ~dz + dww dz + (2.31) (2,32) (2.33) (2.34) Linear Stability Theory We first briefly derive the well-known Orr-Sommerfeld equation. Unlike the equations derived in the above section, the higher order powers and products of u, v and w are neglected. Using the non-conservative form of the governing equations (2.15)-(2.18), we can write the linearized perturbation equations as du du du 0 du du 0 du du 0 dt dx dx dy dy dz dz dP 1 . + Re** dv dv dv0 dv dv0 dv dv 0 al + u a~ + u~a~ + vT" + V~E~ + + W~R~ dt dx dx dy dy dz dz dP 1 A x-^+tsAv' (2.35) (2.36) 16 dw dw dwo aT + & +u~fc dw + Wo-fl- + V ay dwg dy dw + W0 + w dz dw0 dz dP 1 A x-aI + luAw' du dv dw dx + dy+ dz (2.37) (2.38) Since uq, vq and Wq are known, they just act as coefficients, and equations (2.33)-(2.38) are linear partial differential equations for u, v, w and P. The above system can be systematically reduced to a single ordinary differential equation by assuming local parallel base flow, i.e., Vo = 0. (2.39) Further, assume that u0 u0(y), w0 u>0(y), (2.40) and that the three-dimensional disturbance (u,v,w) is a traveling wave with amplitude varying with y and moving along the wall, at an angle 6 respect to the x axis. We can express the traveling waves as = [),*.(). (2.41) where ao is the wave number of the disturbance, c the propagation speed, and w = ac the frequency. Substituting (2.39)-(2.41) into (2.35)-(2.38), we obtain the following linear ordinary system: iaou(uocos8 + wgsinQ c) + u'Q(f>v 17 (2.42) = -iaQpcos6 + ( OL1 2Q)u, iaov(vocos6 + WosinB c) - (2-) iaow(u0cos6 + w0sin9 c) + w'o^v 1 d2 = -ia0psin9 + ( al)4>vn (2.44) ia.Qucos9 + 'v + ia.Q(j>wsin9 = 0. (2.45) Here, prime denotes the derivative with respect to y. Introducing the compact notation Uo = uqcos9 + w0sin97 (2.46) u0 4>ucos9 + (f)wsin9, (2.47) and multiplying (2.41) by cos9 and (2.43) by sin9 and adding them together, we obtain a simpler system with only Uo, v and + tfv = (2.48) 1 d2 iu*(Uo c) + Ufa = -ia0p + (^ a20)0uo (2.49) 1 d2 ia0v{U0 -c) = - If Wq is further supposed to be identically zero, we find that the shape of Uo(y) is independent of 9, thus, the stability computation is also independent of 9 except for the scale factor cos9. Thus, for the oblique disturbance, the base flow Uqcos9 is slower than the disturbance parallel to uq. As Squires theorem 18 asserts, the minimum critical unstable Reynolds number occurs for the case of 2-D disturbance propagating along the same direction as base flow iio(y). If we eliminate ^ and p, the so-called Orr-Sommerfeld equation can be obtained: {( ~ o)2 ~ iRe[{a0u0 w)(^- ao) oloUq]}4>v = 0. (2.51) This equation was first derived independently by Orr (1907) and Sommerfeld (1908). The boundary conditions are specified that the disturbances u and v must vanish at infinity and any walls. From (2.48) we get the proper conditions for (2.51): Channel flows: Boundary layer: 4>v(h) = 4>'v(h) 0. *.(0) = &(0) = o, v{oo) = <#,( oo) = 0. (2.52) (2.53) Both the Orr-Sommerfeld equation and the boundary conditions are homoge- neous, so this problem is really an eigenvalue problem. There are two functional forms for the Orr-Sommerfeld equation, one corresponds to temporal growth and the other to spatial growth: Temporal: f(Re}a0,u>R,uI) = 0, 19 (2.54) Spatial: 9{Re,oc 0R,a0l,w) = 0, (2.55) where subscript R and I denote real and imaginary part of the complex variables, respectively. In temporal growth, Re and a0 (real) are given, and u> = wr + iu>i is the eigenvalue of the problem. Perturbation grows if u>j < 0, damps if u>[ > 0, and is neutrally stable if u>j = 0. Similarly, for the spatial growth, we give Re and u) (real) and solve for o = &oR + ioto,. Perturbation grows with aoj > 0, damps with aoj < 0, and is neutral with c*0j = 0. A solution v(y) which generates one of the eigenvalues in (2.54) and (2.55) is called an eigenfunction. Figure 2.1 gives the neutral curves of the Orr-Sommerfeld equation for channel and flat plate flow. The streamwise and normal disturbance velocity are then obtained by tao(*co40+z4mtfct) i OiQ V i iai)(xcoa8+zain$ct)
(Pv 6
(2.56)
Blasius Similarity Solution
In this section, we discuss briefly the Blasius similarity solution. Zero pressure
gradient is assumed. After ignored the higher order small quantity terms, the 2-D
boundary layer equations then can be obtained as
du dv
dx dy
du du
u------1- v
dx ^ dy
0,
TTdU di 2u
V + v
(2.57)
(2.58)
20

where U = U(x) is the freestream speed outside the boundary layer. Introducing
u(x,y) = U(x)f(v),
V = V\
\JL_
2vx
(2.59)
(2.60)
and considering ^ = 0 for a flat plate, we obtain
ro?)+fMnv) = o.
(2.61)
This equation was first derived by Blasius. The boundary condition is given by
No slip: /(0) = /'(0) = 0,
Far field: /'(oo) = 1.
(2.62)
System (2.61)-(2.62) can be easily solved by conventional numerical methods.
Letting U = 1, the base flow of a flat plate now can be expressed as follows:
u(z,y) =
(2.63)
(2.64)
where
Rex =
U x
v
(2.65)
is the local Reynolds number. Some other useful formulas for the Blasius flow
are:
boundary layer thickness 6:
6 5.0
x y/Rex
21
(2.66)

displacement thickness 6*:
Â£* 1.7208
x y/Rex
(2.67)
momentum thickness 6:
0 0.664
x \JRex
(2.68)
Re
()
Re*
(b)
Figure 2.1. Neutral curve for (a) planar channel, (b) boundary layer flow.
22

Initial and Boundary Conditions for Planar-Channel-Type Flow
For a planar channel, Poiseuille flow is considered as the base flow:
U(x,y,z) = (1 (y 1)2,0,0). (2.69)
The case of roughness elements will be discussed in later chapters. The initial
conditions for the perturbation flow for temporal simulation is
(u,v,u>) = e2d-fteaZ{((Â£u,<^,0)2det(a*)}
+Â£3rf_i2eo/{((^u, P = 0, (2.70)
where u, v and w are eigenfunctions for u, v, and w obtained from Orr-
Sommerfeld equation, e2d, e3t/+ and e3(j_ are amplitudes of the 2-D and 3-D
perturbation waves, aci and a2 are streamwise wave numbers for 2-D and 3-D
perturbation, and ft is the spanwise wavenumber. For spatial simulation, the
initial conditions are
(u,v,w,p) = 0. (2.71)
On solid walls, no-slip conditions are assumed:
{u,v,w)\waU = 0. (2.72)
No boundary condition for pressure is needed since we use a staggered grid for
the numerical simulation.
23

For temporal simulation, periodic boundary conditions are imposed in both
the streamwise and spanwise directions:
V(x + Lx,y,z) = V(x,y,z),
P(x + Lx,y,z) = P(x,y,z),
V(x,y,z + Lz) = V(x,y,z),
P(x,y,z + Lz) = P(x,y,z). (2.73)
For spatial simulations, the inflow boundary conditions are obtained from the
2-D T-S wave and a 3-D oblique wave pair. Because of the specialty of staggered
grid, the condition can be written as:
e2dReal{(f>u2de tt,i}
^3d- Real {(j>U2d- e*( _/3z },
e2dReal{(f>v2dei^a^-^}
e3d+Real{v3d+e^-a^z-ut)}
e3d-Real{4>v3d-^-a^2-ut)}
e3d+Real{w3d+e-a^+l3z-ut)}
e3d-Real{(j>w3d-^-a^-^-^}. (2.74)
The outflow boundary conditions are much more complicated, and will be dis-
cussed in detail in later chapters.
24
u(0,y,z,t) =
+
+
H2">y>z>t) =
+
+
M^-,y, *,t) =
+

Periodicity conditions are also used in the spanwise direction for 3-D temporal
and spatial simulations because the flow in this direction is not confined, and the
laboratory experiments also support the spanwise periodic structures.
Initial and Boundary Conditions for Flat-Plate-Type Flow
For flows over a flat plate, the Blasius similarity solution is applied as the
base flow. For rough plates, the Blasius solution will be used as the initial guess
of the base flow. Details will be given in later chapters.
Only spatial simulation is conducted for smooth or rough plates in this study.
Thus, similar to the last section, Benny-Lin type disturbances are imposed at the
inflow boundary:
e2dReal{u2de *wt}
e2dReal{v2de^-a^-ut)}
e2d-Real{^d-^-a^-^-ut)}
e3(f+ Real{W3d+ e*( +^Z-L,t)}
e3d- Real{(f>w2d- } (2.75)
u(0,y,z,t) =
+
+
+
+
=
+
where u> is the real frequency of the disturbance, /3 is a real number that rep-
25

resents the spanwise wavenumber, and a = a# + ia/ is the streamwise complex
wavenumber obtained from linear stability theory.
A no slip boundary condition is applied at the solid wall. According to the lin-
ear stability theory, the disturbances vanish at infinity, so, we obtain the bound-
ary conditions at far field given by
u(x,y > oo,z,t) = 0,
v(x,y > oo,z, t) = 0,
w(x,y > oo,z,t) = 0. (2.76)
Also, no pressure condition is needed at the inflow, solid wall, and far field since
a staggered grid is used.
The outflow boundary condition will be discussed later.
26

CHAPTER 3
TRANSFORMATION OF GOVERNING EQUATIONS
In this chapter, a conservative form of the incompressible Navier-Stokes equa-
tions is presented.
Vector Form of the Navier-Stokes Equations
The dimensionless 3-D incompressible Navier-Stokes equations in Cartesian
coordinates may be written as:
du duu duv duw 1 d2u d2u d2u dP (3.1)
~dt dx + dy 1 dz Re ^dx2 + dy2 + dz2' + dx = 0,
dv dvu dvv dvw 1 ,d2v d2v d2v dP (3.2)
dt dx + dy 1 dz Ye dx2 + dy2 + Ih2' + dy = 0,
dw dwu dwv dww 1 , d2w d2w d2w dP = 0, (3.3)
~di + dx dy + dz Re^ dx2 f dy2 + dz2' + dz
du dv dw (3.4)
dx + dy + dz = 0.
It is convenient to rewrite them to a general vector form. First, we introduce the
following vectors:

V w
vu , 0 = wu
v2 + P wv
wv w2 + P
Then equations (3.1)-(3.4) can be rewritten in a simple form:
dU_ (PE (W dG _ -
dt dx ^ dy ^ dz Re
where
A d2 d2 d2
dx2 dy2 dz2
is the Laplacian operator in the physical (x,y,z) space.
(3.5)
(3.6)
General Transformation
In order to deal with arbitrarily shaped boundaries, it is often convenient to
employ general curvilinear coordinates Â£, 77, which are related to Cartesian
coordinates by a transformation of the general form (see Figure 3.1)
* = (Â£,77,0,
V =2/(Â£.?7>0>
2 = z(Z,V, C), (3-7)
or the corresponding inverse transformation
Â£ = Â£(^,y,z),
28

v = iK*y*)i
C = C(*.y.*)-
(3.8)
y v
Figure 3.1. Transformation from the physical (x,y,z)
space to computational (Â£,*7,0 space.
Using the chain rule, the partial derivatives become
d . d a + Cx d
dx I My H UJ* II +,a? dC
d d a + Cy d
dy ~*ydt dC
d d d +c d
dz *zd( + Vz^~ 011 ''of
(3.9)
If the inverse analytical transformation (3.8) is easy to obtain, the metric Â£x,
Tjx, (x, Â£y, jjy, Â£y, , y z, can be obtained directly. Otherwise, the following
procedure is recommended. Rewriting the differential expressions in matrix form,
29

we have
1 t. tv tz dx
drj = Vx Vy Vz dy
1 V/* 1 Cx tv Cz dz
and, similarly,
dx 1 B H a H /\ 1 di
dy = vt Vv vc drj
dz zt zv ZC dt _
Defining the Jacobian of the transformation J as
d(x,y,z)
tx ty tz
Vx Vy Vz
Cx Cy Cz
from (3.10) and (3.11) we obtain
tx ty tz xt Xr, X(
Vx Vy Vz vt Vri Vt
Cx Cy Cz H ZV Z(
= J
yr,Z( y^Zr, -(air,Z( ~ X^) Xvy( Xflr,
(vtzc ~ vczt) xt.H xczt ~(xm xtVi)
ytzn VetaZt -(:C(Z ~ Xnz{] X^ X^
(3.10)
(3.11)
(3.12)
tx = J(yvz<-y
30
and, thus,

x(zrj)>
Â£z J{x r)2/( XCVt])j
*im= -J(y(zc~y
Vy ~ Jfaz< Â£Â£Â£())
Vz = -J{xtyc xcyt),
Cx = J{ytzTi ~ yr,z(),
Cy = ~J{xtzri ~ xt]z()i
Cz = J(xtyr,-xT,y()- (3.13)
After the proper transformation, a domain with general geometry can be trans-
ferred to a rectangular domain in the computational (Â£,77, Â£) space.
General Curvilinear Coordinates System
Based on the above transformation, we now define the new vectors in the
computational space:
Ui = U
J
Ei = j{E(. + F(y + Gl.,),
A = j(Et!x + FVy + &Vz)
Gi = ~j(E(x + F(y + G(z).
(3.6) then becomes
QUX dEj. 8F1 dGx 1 A r7
7T+aT+aT + ir-ife 1 ''
31
(3.14)
(3.15)

where Ai is the physical Laplacian operator transferred to the computational
(Â£>^0 space:
Ai - +1* drj + + Vx Qji + Cx g(>)
d d d... d d fl.
+{ty QÂ£ +7ldj]+ Cy d(. ){Cy d(f +Vydv+ Cy ^ )
+(6^ + ^ + +Vq^ + C.q^)- (3.16)
Though (3.15) looks very simple, it is not convenient to discretize. To obtain
equations that are more easily discretized for numerical simulation, we define the
new variables 17, V, and W in the computational (Â£, 77, Â£) space:
+ vÂ£y + wÂ£z), (3.17)
j(u7]x + VT]y + W7}z), (3.18)
j(uC* + vCy +wCz), (3.19)
and consider both variable triples u,v,w and U, V, W as unknowns in the compu-
tational space. Our objective is to develop relatively simple expressions for the
governing equations in general coordinates that are easy to discretize. Using the
above transformation, we rewrite the vectors Ui, E\, F\ and G\ in the following
form:

1 u
u uU+ ^
, Â£i =
V vU+^
w wU+Zj*
32

V W
uV + ^ uW+^
,C?i
vW + ^
wV + ^ wW + ^
The governing equations then become
du r,duU duV duW.
Â¥+J(-9r+-ar+-9r)+
.. d d 9 1.
(^9Â£ + Vm&i+ CxdC)P ~ ReAlU ~
(3.20)
(3.21)
dv r,dvU dvV dvW.
^+J(^-+^r+^r)+
{ty-Qi + Vyg^ + Cyg^)P ~ ^ = 0,
(3.22)
dw r/dwU dwV dwW.
3<-+J(^-+^r+^r)+
^*9Â£ + Vzg^ + Czg^)P ~ feAlW = 0
(3.23)
dU dV dW n
----1-----1----- 0,
dÂ£ dV d(
(3.24)
The system (3.17)-(3.19) and (3.21)-(3.24) has 7 unknowns, u,v,w, U, V, W and
P. Similarly, we can obtain the perturbation form of the governing equation as
follows:
du Tf d[u(U + U0) + u0U] d[u{V + ko) + uQV]
ai+ ( ai + +
33

8Kw + Wo) + oW]. ,/f a a a
--------g-(------) + ((-gj + V.fij + (gf)P
~TeAlU =
(3.25)
Bv a\v(u + va) + vu\ g[v(y + Vo) + voV} ,
8t+'y( 8f + Aj +
8KW+W0) + 011 ,, 9 a a
------F(------) + ^
~tAlV =
(3.26)
dw d[w(U + U0) + w0U] d[w(V + V0) + w0V]
m+J{--------Â¥(------+-------Â¥,------+
-------F<------) + (6^ + V.frj + (,F)P
-i-A.w = 0,
(3.27)
dU dV dW
di + dr} + d(
= 0,
(3.28)
where I7o, Vo and Wq are base flow solutions defined in the same way as (3.17)-
(3.19). Combined with (3.17)-(3.19), this system also has 7 equations and 7
unknowns for the perturbations. From (3.16), we define the operator Ai in the
computational space as
Al (Â£x + il + (\)-qÂ£2 + + Vy + ^0^2 + + Cy + 2
q2 q2
4'2(^x7/x -)- Â£y1]y + tzTJz) "V ^(Â£*Cs "I" Â£yCy H" CzCz)
d2 d
+2{Vx(x "f VyCy *b VzCz) (Â£x Â£yy
34

+(7?x* + Vyy d" Vzz) + (C** d" Cyy + Czz)
(3.29)
Because of the high accuracy requirement in flow transition simulation, an ana-
lytical transformation is recommended, so there will be no extra artificial errors in
the system (3.17)-(3.19) and (3.21)-(3.24) for the base flow, or (3.17)-(3.19) and
(3.25)-(3.28) for the perturbation flow introduced by the transformation. Note
that no orthogonal assumption for the curvilinear coordinates is made, so the
above governing equations can in principle handle any kind of curvilinear coordi-
nates, provided the transformation exists and is one to one. The new emphasis
here is that the system has 7 unknowns (u,vtw,U,V,W and P) instead of the
usual 4 unknowns (t/, V, W and P) suggested by others.
35

CHAPTER 4
NUMERICAL METHODS
FOR TWO-DIMENSIONAL PROBLEMS
The Navier-Stokes equations and their transformations given in Chapter 2 and
Chapter 3 are nonlinear, time-dependent and three-dimensional. Except for a few
cases, analytical solutions are impossible to obtain, and numerical methods must
be used. For numerical simulation of transition, there are several requirements
we must meet:
phase-accurate discretization of the convective terms
high resolution for regions close to the solid wall
stability of the discrete system
proper outflow boundary conditions that remain non-reflective even in the
presence of nonlinear wave interaction
In addition to the above basic requirements, efficiency, memory storage, and
vectorization of the calculations must also be considered.
In this chapter, we first present an accurate method that meets these require-
ment for two-dimensional problems. The procedure uses fourth-order central

finite differences in space and second-order backward Euler in time. The numeri-
cal procedure is conducted on a staggered grid, so the resulting schemes are quite
complicated to describe. At the outflow boundary, the buffer domain method is
applied.
u | v O P
Figure 4.1. Staggered grid structure for 2-D channel flow. For
the primary pressure point inside the computational domain,
the index is marked from i = 2, .ni 1, j = 2, , nj 1.
Discretization on a Uniform Staggered Grid
A uniform staggered grid is first employed to discretize the computational
domain in the case of flow past a smooth plate channel. Figure 4.1 depicts the
2-D uniform staggered grid structure, in which the discrete pressure P is defined
at cell centers, the discrete streamwise velocity u at centers of vertical cell walls,
and the discrete normal velocity v at centers of horizontal cell walls.
For planar channel flow, since the base flow is known (the Poisieulle flow), only
37

the perturbation equations are needed to be solved. The governing equations
described in Chapter 2 can be simplified to two-dimensional equations. The
computational domain is
x E [0 ,mLx\,
G [0,2], (4.1)
where Lx is the nondimensional wavelength of a 2-D T-S wave and m is an integer.
The base flow in this domain is the Poisieulle flow
(uo,o) = (1 ( l)j,0),
(4.2)
and the resulting simplified governing equations can be written as follows:
du
dt
+
d(uo + u)u duv
dx
dy
1 d2u d2u du0 dP _
Re dx2 ^ dy2 ^ dy ^ dx
dv d(u o + u)v dw
dt dx dy
_ J_(^1 -n
Re dx2 dy2 ^ dy
du dv
dx dy
(4.3)
(4.4)
(4.5)
Fourth-order central differences is applied here. Because of the use of staggered
grids (see Figure 4.2, which illustrates the indexing scheme we use), we need the
following two forms for the first derivatives of a generic function (x,y):
at X: i:
* i
A
^dx^-'
12A*
+ 0( Ax4),
(4.6)
38

di+lj + 274>i+u 27^_t,. + 4>i_u
(to)u -----------------------24Al--------------------+ 0(Al >'
(4.7)
The second derivative can be written as
A
Kdx2)i-^j
-ti+ij + 16i+itj 30&_i.,- + j,j-
12Ax2
+ 0(Ax4). (4.8)
The formulas for the derivatives in the normal direction are similar.
The order of the finite difference approximation is reduced to second-order at
the locations adjacent to the inflow/outflow boundaries to avoid very complicated
one-sided high order finite difference approximations. Thus, the formulas are
expressed as
h+kJ ~ h-u
() i = + 0(Ax2),
oa_ ~ v )i
2 Ax
1dxH~
5V k+U-Wi-U + h-U t
9x2 a,2'
Ax2
+ 0(Ax2).
(4.9)
(4.10)
The normal derivatives were also changed to second-order for the points adja-
cent to the solid walls. The derivatives for the y-momentum equation (4.4) are
approximated by
(jph-i = ~ 2t.? + ^ + o(V).
Aj/
(4.11)
(4.12)
Since the no-slip boundary condition must be satisfied, we do not use the image
points (located at y = ^ and y = 2+ ^ ), but derive the required second-order
derivatives directly:
39

y = % (i = 2):
,d2s 4 , . ,
\ ~Qy2 )*<2 ~~ 3^*-3 4 (4.13)
y = 2- ^ (j = nj 1):
A
\ Qy2
r0t,nj-2 4,- nj_i,
(4.14)

ui-U PiJ
*

Figure 4.2. Indexing scheme for the staggered grid, illustrating
primary variable locations for the ijth cell.
Before we express the resulting discrete system, it is necessary to clarify the
special indices used for the staggered grid. The spacings are determined by
Ax =
mLx
(ni 2)
Ay =
2
V
(4.15)
where ni and nj are numbers of grid points in the streamwise and normal di-
rections, respectively. Using (i,j) as the coordinates at cell center (see Figure
4.2), then, the cells in the computational domain is marked by i = 2, ,m 1,
j = 2, ,nj 1. Based on the above setting, we obtain the following difference
scheme written in general form:
40

x-momentum equation
Ai+liui+bi + Ai+biui+kd + Ai-biui-hi + ,i
+Ai-b3+2Ui-h,3+2 + ^i-5,J+lUi-j,J+l +
PAi-\,j-2Ui-\,j-2 Ai-\,jUi~bj Di+l,jPi+l,j
+D;-i,j-Pi-i,j + Di-2,jPi-2,j DijPij Qu .
1 Jw
(4.16)
The corresponding point for the above equation is (t |, j), which will be indi-
cated superscript of the coefficients when there is danger of confusion (Note, for
example, that A*. I3. and A!+1'3. are generally very different).
1 JJ 1 JiJ
y-momentum equation
^+2,j-pi+2,j-i + Bi+l,j-%v
+Bi-2,i_t^_2ij_i + Bu+iViJ+l + Bitj+x_viJ+ i
1 Bi,i-kvi,i-\ + Eu+
f^'i.J-lPij-1 d" Pi,j-2Pi,j-2 ~ EijPiJ = Qy. -_jl
j
j+i
with central point (z,y |).
continuity equation
(4.17)
Ri+\,iui+y + Ri+biui+ii + Ri-yui-\,i
-Ri-biui-bi + 5i.i+i^.i+i +
(4.18)
with central point (i,j).
41

Figures 4.3-4.5 depict the relevant variables for the i-momentum,
^/-momentum, and continuity equations, respectively. In those figures, variables
at other than primary locations are obtained by fourth-order interpolation. For
example, we have (see Figure 4.6 (a))
vi-u = [9Km +
+ Wi+l,j+|)]/32> (4-19)
with truncation error 0(Ax4 + Ay4).
For points adjacent to solid walls, the accuracy of the interpolation in the
y-direction is reduced to second order to avoid the need for ghost points which
may reduce accuracy. The interpolation formula in this situation is (see Figure
4.6 (b))
vi-iJ = + viJ+k + + Vi-U+^)
~{vi-2,j-% + vi-2tj+% + + ^t+ij+j)]/32* (4.20)
with truncation error 0(Ax4 + Ay2).
42

Vi~1,J+2 3"a> J+2
V-a ii+l_
o U 3 Vi-i. i _^*-a-J o 2 - 3+U 0 -
Pi-2,j Pi-1, j Pi,3 Pi+1,3
Vi-J-1 i-J.J-1
2_
Figure 4.3. Neighbor points for the discrete xmomentum equation.
C V- 3 Pi,3+1 + k
Vi-*,i-\ _L Vi-l,j-\ c V- 5 Pi,3 Vi+1,3-', Vi+2,j-\
Ui~l,3~k Vfy-i Ui+1,3~\
c t/. 1 Pi, 3-2 -1
Figure 4.4. Neighbor points for the discrete ymomentum equation.

Figure 4.5. Neighbor points for the discrete continuity equation.
vi-2 ii+| Vi-1 2 Vi J+2 Vi+1 j+!

1 Vi-1 3~\ 1 1
Vi- !,i+j Vi- (a) i.j+* v >j+2 Vi+ 'li+j
Vi-U
Vi- Vi~ U-k v U-k V<+ i.i-5
(b)
Figure 4.6. Interpolation for v at a u point (a) interior point,
(b) point adjacent to the solid walls.
44

The formulas for the coefficients for the interior grid points are listed below:
x-momentum, equation
Ai+bi
Ai+li
Ai-y
Ai-\,J =
Ai- j.J+2
Ai-bi+1
Ai-%,j-l
Ai-h,j-2

A*+i,j
Di-U
Di-2,j
Di,i
Qu. ,.

12iÂ£eAx2 12Ax
4 2
3fleAx2 3Aa^U<+*J' +
ZReAx2 + + M-P
1 1 .
12iZeAx2 12Aaj^i"*',' + ,t0i-^^
1 1
+ rrr-:V;_i
12ReAy2 12Ay f_2lj+2
4 2
3ileAy2 3Ay =
4 2
ZReAy2 + 3AyVi~^'5~1'
1 1

12J?eAy2 12Ay **
3 5 1 1 X
n D v A _2 A -.2 /
2 At 2Rey Ax2 Ay2
1
24 As
27
24Ax
1
~ 24 Ax
27
24 Ax
-4u? + iin" .
*~2'J t-2'J
2At
2yiv- i
-jij
y-momentum equation
_ 11, .
~ uRcAx2 + 12Ax{Ui+2j-* +Uow-h>*
4 2
5i+1-M = 3iZeAx2 _ +o<+1J_j),
45
(4.21)

Bi-i,M
Bi,3+1
Bi,j+k
Bi,i-I
Bi, J-J
Ei,j+i
Bi,j- 1
Ei,j-2
Eij
En-i
w 2
+ nK-ij-i + o_l _*)>
ZReAx2 ZAx
12-ReA2 12A^i-2,i_* +
1 1
12J?eAy2 + 12AyV,i+^
4 2
3ReAy2 3AyVi,j+*
4 2
3ReAy2 + 3Ay1'i*j~*
1 1
lZReAy2 ~ 12AyVij~^
3 5.1 1 .
+
2A t
1
2i?evA2 Ay2'
24 Ay
27
24Ay
1
_ 24Ay
27
24Ay
-4u?. i +r?711
.J3
2Ai
(4.22)
continuity equation
Bi+bi = -Ri-\,i =
Xi-U - Ri+5,3 2
si.i+1 = Si, 3 = hJ 2
sj-i ~ Si,3+k ~ 2
Qrrnj = 0.
24A
27
27
(4.23)
Though Qmij is zero at the finest grid level, it will be nonzero at coarser grid
46

levels, so we still keep it in our formulas.
It is easy to describe the coefficients for the second-order accuracy scheme,
but for those points adjacent to the solid wall, the coefficients are special. Only
those coefficients of u in the normal direction for the xmomentum equation are
irregular, so they are listed below:
lower wall:
Ai-,j+2 ~ Ai-$,j-2 ~ Ai-$,j-1 -
, 4 1
ZReAy2 2Ay^'i+1
3 1 / 5 4 1
2At + Re 2A2 + Ay2^ + 2Ay^li+1
upper wall:
(4.24)
Ai-\,j+2 ~ 2 ~ - 0>

Ai~\,i =
322eAy2 + 2Ay^,i*-,'-1,
4
2At + fle^A2 + Ay2^ 2Ayt,i"^'"1-
(4.25)
Discretization on a Stretched Staggered Grid
For flow over a flat plate, the computational domain is a y-direction truncated
semi-infinite domain. Since the subdomain that we are most interested in is very
close to a flat plate (scaled by boundary layer thickness), local high resolution
near the flat plate is required, and thus a y-direction stretched grid that becomes
increasingly dense near the solid wall is employed (see Figure 4.7).
47

Figure 4.7. 2/-direction stretched grid.
To meet the above requirement, an analytical mapping in the wall-normal
direction is employed, i.e.,
y(v) =
Umax^y
ymax & Vmaxijimax y)
or its inverse mapping
ymaxy{& I" Umax
ymax( where ymax is the height of the computational domain in the physical coordinate y,
ymax is the height of the computational domain in the computational coordinate
77, and o' is a constant used to adjust concentration of the grid points. Thus, we
can use the mapping
)
x = (,
y =
ymax&y
(4.26)
max& "t" ymax {ymax yi)
to transfer the stretched grid in the physical {x,y) plane to the uniform grids in
the computational (Â£,77) plane. With this transformation, the variables in the
48

computational plane are now defined as
u =
V = v,
(4.27)
(4.28)
and the resulting governing equations, if we represent them by using only U, V
and P in the computational (Â£,77) plane, are
U-momentum equation
dU 1 d(UU + 2U0U) d UpV + UVq + UV.
9t + y di + drj yv
1 ,d2U 1 d2 ,U. d ,U^ dP
Re^d? + yr,dT,>(yJ + yrVvvd7,{yJ) + yvdt ~
V-momentum equation
dV d(U0V + UV0 + UV) d{ 2V0V + VV)
Vr) dt+ di + dr)
1 d2V 1 d2V dV, dP
Re*0 V" d? + y dr,2 + Vt,Vvv dr,} + dr, ~
(4.29)
(4.30)
Continuity equation
dU dV
d(+ di ~ '
(4.31)
The required derivatives in the governing equations can be obtained analytically:
Vmaxymax&(,& "h 2/mox)
Vn =
Vw =
[l/mor1 "h ymax(j)max ^/)]^
2r,
mate "t" ymax)
(4.32)
(4.33)
ymax (& + y)3
Based on the above derivation, all terms are considered in the uniform grids of the
computational (Â£,77) plane only, and, thus, a procedure similar to that described
in the last section can be used to derive the finite difference formulas.
49

Figure 4.8 depicts the staggered grid structure used for the spatial discretiza-
tion. Note that the discrete pressure Ph is defined at cell centers, the discrete
velocity Uh is defined at centers of vertical cell walls, and the discrete velocity
Vh is defined at centers of horizontal cell walls. Moreover, the discrete continuity
equation is associated with cell centers, the discrete Umomentum equation is
associated with Uh, and the discrete Vmomentum equation is associated with
Vh.
Figure 4.8. Staggered grid for flat plate flow in computational (Â£,77) plane.
The resulting discretized equations, if denoted symbolically, can be written

50
as

+Ai_iij_2Ui_itj_2 - + Di+uPi+u
+Di-i,jPi-i,j + Di_2 Bi+2,j-Vi+2,j- \ + + Bi-l,j- j
+Bi_2ij_i_Vi_2tj_i_ + Bij+ iViij+i + BiJ+iVi)j+i
+Bi,j-- Bi,i-\Vi.i-\ + Ei,j+Ipi,j+1
+Eiij-1Pitj-i + Eitj_2Pi,j-2 EijPij = Qv. ._^, (4.35)
Ri+yUi+ij + Ri+LtjUi+i_tj +
+ Ri,i+\Vi,i+\ + Ri,mvi,j+h
+Ri,j-\Vi,i-\ ~ Ri,i-kVi.i-k = Mi>r
(4.36)
V^J+*
' Ui~j,3+1
Ui_. -* a,J O U{_ 3, Vi-\j -Ui~l3 0 Si+1>3 O
Pi-2,j Pi-1,3 Pi,3 Pi+1,3
Vi-l3~ 1 *Ui-\,3-1
V^J-2 *Ui-\,3-2
Figure 4.9. Neighbor points for U momentum equation.
Since the process is similar to that in the last section, we only give the coef-
ficients for the U-momentum equation (Figure 4.9). For the interior grid points,
51

k[M MjM U[H
the formulas can be obtained as

^*+ k ,j
A{_3
* a iJ
Ai_i
* a *j
\2Re*0At2 + 12A(yVj +
4 2
3Re*AÂ£2 ~ 3A\$yTtj + 2Â£/o.+ ^.P
+ + 2%-l>
12J2eJAfa 12A^{Ui-^ + 2Uai-\J'
i+2
12i2eoA7/2yfJj+ayTJi
1 /y , y \ VviVvyj
12A _ U-^-+a+ 1Me*A .
,j+i
.j-i
ZRelAtfy^y^
2 fy , y \ ^yviVyyj
3Atj2/. / 2-J+1 ZRelA-qyj,.
^j+i
4
jj+i
ZRelA^y^y^
2 + y.2s,',*!"
' 1 2 A fi1'
3A*nhij-i
ZRelA-qy^
,i-2

-D+l,3
A,i

-2'J
12iZeJ AT)2yVj
12A7/i/r;._3
12iEe5A7/yJJi_3
3 5 1 1
2 At + 2Re*^Ai2 + At/2^/
n. ^
Di-u =
i~2,i 24AÂ£
27^-
24AÂ£
l^-!1. 4Z7 i . 1 ,+a%. . ,+1ET0.
-j.J -?.J , 1 r ~2J+Z V-fr.j+a g 2|3+1
2 At 12 At/ Vrij+i Vvj+i
V^L^Uo. Vi_L i_2Uo. .
g 1 3'J 1 t-f-.J-l ^ * 3'J ^
Vrii-i Vnj-3

(4.37)
Here, superscripts n and n 1 are used to indicate previous time values and
52

the superscript n + 1 is understood for the current time, which is dropped for
convenience. Similar to last section, the variables at other than primary locations
are obtained by fourth-order interpolation in the computational (Â£,7/)-plane.
At the adjacent of upper and lower boundaries (say j = 3 and nj 2) except
for the inflow and outflow boundary points, the finite difference scheme in the 77-
direction is reduced to second order, while fourth order accuracy is maintained in
the ^-direction to preserve the phase accuracy. The resulting modified coefficients
for {/-momentum equation then are:
Ai-bi+2
- 2 0>
Ai~k,}+1
Re*0Av2ynjyr,i+1 2AyyT).^ ^ +
VnsVyys
A:
i-bj-1
2Re^Ai]yJ,j+1
1 +^
Ai-y
RelA^y^y^ 2A-qylij_^
_____VrijVyyj
2Re^ArjyT)j_1
3 1 SV<-bi+\ + _ W-y-k +
2At 2A77 yv.+ k yv._h
1 r 5 2 %
* O A Â£2 A_2..2 /
Qu. 1.
'~2<1
Re*Q 2A(2 A-rj2y2.
U?~i 4 U? .
2At
1
AtI[ 2 y,
(4.38)
M
with and obtained by fourth-order interpolation (see Figure
53

4.10):
= m.j+^+w*) (Wi+Wi))/^
Vi-L.L
1 jw 3
(Wi + Wi) (Wi + Wi))A. <4-39)
and similarly for Vo , , and Vo , , and other coefficients are kept un-
j-ifj
changed.
Figure 4.10. Interpolation for V{_ij+i and Vi_kij_k
at the adjacent of upper and lower boundaries.
At the lower boundary (j=2), since a no-slip boundary condition is used, the
associated derivatives need to be modified for (4.38):
A- i 3=0,
* 2 2 9
A- i =
4 aiJ+2 y^jVyyj
2.J+1 ZRel^y^y^ 2AW^+i ' 2ReoAWw+i'
, 3 1 ^-y+} + Voi~^
'-H 2 At + 2A77 y^
1 / ^ ^ \ Vm
,*UaC2 ' A2.,2 '
(4.40)
Re*Q 2AÂ£2 A-q2y*. 2Re^Ar]
At the upper boundary (j nj1), a far field condition is assumed, which implies
that all the disturbance should have disappeared, so we can also let U = V = 0
54

there, which yields
Ai

0,
4 | +
2AyyVj_^
J______1 (Vi-bH+Voi-l;-0
2A t 2Arj yn._^
i 1 ^ 5 ^ 4 . ^7w/j
+ Ei*(2A^ + + 2RelAv'
yrgVyyi
2Re^AriyVj_1
(4.41)
Â£ 0 Â£ Loriginal Â£ Ltotal
Figure 4.11. Extended computational domain.
Outflow Boundary Treatment
Outflow boundary conditions have been the focus of study for spatial simu-
lation of flow transition by many researchers. Taking advantage of the staggered
grid leads to a fairly effective approach, which we now describe. First, a so-called
55

buffer domain technique developed by Streett Sc Macaraeg (1989) is applied to
our problem. Thus, a buffer domain is appended to the end of the original outflow
boundary to smear all possible reflections from the buffered outflow boundary (see
Figure 4.11). The problem is generally that the conventional buffer domain is too
long (usually four to eight T-S wavelengths), which greatly increases computation
costs. Our goal is to maintain accuracy in the original computational domain,
but eliminate all the possible reflection waves in a very short buffer domain. To
realize the above goal, the governing equations in the buffer domain should be
parabolized to have strictly outgoing waves. Thus, an initial buffer function 6(Â£)
is introduced here and applied to the streamwise viscous terms:
d2U
dÂ£2
d2V
d?
KO
KO
d2U
d(2
d2V
di2'
(4.42)
6(Â£) is monotonously decreasing function from 1 to 0 so that the upstream effects
of the streamwise viscous terms will gradually disappear in the buffer domain.
The essential feature here is that all damping mode disturbances at the buffered
outflow boundary become zero. To understand this, we rewrite (4.42) as
Ki)
d2U
d2U

(4.43)
Since 6(Â£) > 0 at the buffered outflow boundary, and accordingly > oo,
then we can also consider that the buffered outflow boundary is compressed
from Â£ = oo by the function ^/f>(Â£). Now, clearly, if the disturbances are stable
56

(damping modes), they will vanish at Â£ = oo. To treat growing or unstable modes,
we need a second buffer function bRe(t) that reduces the Reynolds number in the
buffer domain gradually to less than the critical (or subcritical) Reynolds number
and makes all the perturbation modes become damping
J;____^ &jfe(Â£)
Re Re
Thus, the new modified governing equations in the computational (Â£,f/) plane
become:
dU_ 1 d(UU + 2U0U) d UqV + UVq + UV
a + Vr, + dy yv
dt
bRe_^ 1* Â£ 0,Â£vv 9P
Re 3( d(2 yr,dr,^yri) + yr,Vyvdr,{yJ) + yndC
dV d(U0V + UV0 + UV) d{2V0V + VV)
Vr) dt + d( + dy
bRe d2v id2v t dV. dP
Re% ^ *Vn dP + yv dp + yr,Vyv dy + dy
dU dV
dÂ£ + dy
0, (4.45)
0, (4.46)
0. (4.47)
The buffer functions are chosen as follows:
KO
tanh(Ltotai-()
tanh(Liujfer)
1
I*original ^ t ^total >
0 ^ t ^original >
WO
/
<
L?.. (Â£ ^original)2 + 1
buffer
1
Loriginal t Rtotali
0 ^ t ^original
(4.48)
It is clear that the first function decreases from one to zero very rapidly as one
moves from the original outflow boundary to the buffered outflow boundary.
The second function increasing from 1 to c + 1 is a quadratic function that is
57

continuously differentiable at the original outflow boundary. Note that the total
effect of these buffer functions is that, toward the buffered outflow boundary:
the momentum equations become increasingly convection dominated in the
Â£direction, while the equations generally become parabolic; and
the momentum equations generally become highly diffusion dominated in
the tjdirection.
This treatment makes the outgoing waves propagate forward without reflection in
the Â£direction, and any oscillation in the tjdirection will be effectively smeared
in the buffer domain.
Finally, we need to specify the buffered outflow boundary conditions under
these modified governing equations. The parabolic character of the above equa-
tions requires only two boundary conditions. As mentioned before, we have the
disturbances tending to zero at the buffered outflow boundary (which is actually
located at Â£ = oo), so one condition is
P = 0. (4.49)
This is a very important condition since the elliptical character of pressure has
not been modified in our new governing equation. Any improper condition for P
may cause trouble. For the second condition, we use the traditional extrapolation
method for V, i.e.,
d2V
58

Though this condition may not be so accurate, accuracy of solutions in the buffer
domain is not so important, and the main concern is that there must be no
reflection wave traveling back to the original computational domain.
Referring to Figure 4.12, condition (4.49) is imposed directly on P by defining
Pm- 1J = 0. (4.51)
This is an implicit condition. With it, the discrete continuity equation associated
with Pni-i,j is then used to define U at the buffered boundary:
TT TT ^ni l,j+ 2 l,i+ 2 A t
------------------------A<-
(4.52)
Condition (4.50) is imposed by determining V at ghost points just outside the
boundary:
V
mj-
= 2V
" m-
i ,j-k
- K;
-2
(4.53)
yni -l|j+2
c Pni-l,j_

Vni-2,j-\ 1 rH i

l.
2
Figure 4.12. Buffered outflow boundary points.
Note that the above treatment is only suitable for the perturbation equations.
For the original equations, although the idea is similar, the resulting boundary
condition are treated differently as we will discuss for 3-D problems.
59

Consistency and Stability
There are a number of finite difference methods for time-dependent flow prob-
lems, like FTCS (forward time, central space), BTCS (backward time, central
space), FTFS, Lax, Lax-Wendroff, Leapfrog, etc. Each scheme has its own mo-
tives and corresponding relationships with the original governing equations, and
therefore produces different numerical effects. Thus, it is necessary to analyze
the consistency and stability of our scheme.
Consistency
We take the U momentum equation as an example:
dU_ 1 d(UU + 2U0U) d UqV + UVq + UV
+ y di + drj yr,
1 ,d2U 1 d2 ,U, d ,U dP
Ret( d? + yr, dif(yr,) + ^ dr,{ y^ + ^
= 0.
(4.54)
Now suppose all variables at the required discrete grid points are obtained exactly.
Using (i,j) as the index of U (for simplicity, here we will use (i,j) instead of
(i |,j) in the above subsections for the central primary U points, and the
indices for other variables are also modified, see Figure 4.13), our scheme (similar
to the so-called BTCS scheme) can be expressed as
3t^-4 VZ + VTJ' 1
2At---------+ + 2U>+*jDW
+8 (U + 2C/0)i+i1jt/i+iij 8(17 + 2Uo)i-ijUi-itj + (U + 2Uo)i-2,jUi-2,j]
, 1 r (v + ^>)t,j+2 TT
[-------------Ui,j+ 2
12 At;
2/tj,+s
60

..(V + HOy+ip n(V + F.ki-i (V + V'h-,
Vm+i Vrij-i y^j-i
! 1 r U0i,j+2Vi,j+2 + gPojj.n^.J+1 g^.y-l^.j-1 + t/0.,,-a^rt,j-2j
12A77 2/ijj+j y^i+i y*iji vvj-2
+ Ee* [l2A^(_C/i+2,i + 16t/i+1-J 30^*.j + 16^1I,j -
+------+ 16^1 30^ +
Vr,i 12A772 yi+2 2/tjj+1 Ihv Jfc-i w-j
12A77 3/rjj+a 2/rjJ+i 2/iy-i Vui-i
+^(~F^ + 27P;HJ 27F-~U + = * (4-S5)
Vi,j+2 ~*Ui,j+2
Vi,j+1 ~Uij+1
2 Pi-ij 2 0 Pi-kj Khj -Â£/. 0 P. 1 rt+ I.j Pi+l,i
Vi,i-1 -Uij.,
Vij-2 cs 1 Â£
Figure 4.13. Neighbor points for BTCS scheme at Uij point.
Based on (Xi,yj,tn+i), using Taylor series, we can rewrite the above discrete
Unite difference equation term by term in a modified form. Since
Un = Un+1-UtAt+\uttAt2-lutuAt3 + ---,
z 0
Utl-i = t/n+1 2UtAt + 2UuAt2 \utttAt3 -\-, (4.56)
6
61

we get
3E7?, + 1 40?, + E/?,"1 ,dU.n+1 , 1__ A 2
----2Af---- = <
+
\n+l
hi
(4.57)
For the spatial terms, we use
Uii,j Uij (dx )uAx + 2 (dx2 )i,jA*2
, i .5*0, A 3 i .a40. A 4
6(5x3)iljAa: +24(5*4)i,iAa;
, 1 ,d5U, A b 1 ,360, A e
12n(/3^)^Aa: + 72(\( dx*^iAx +"
Ui2,j Uitj 2( ^ ).jAa: + 2( 5a,2 )*.iA2
, 8 33U A 3 16,040. A 4
Vdz3^* + 24 ( dx^ijAx
, 32 ,0B0X 5 64 ,060, A 6
120 ^ 9s5 ^ X + 720 ^ dx6 ^ +
Note that
0+i,j 0i-i,j -
0+2,j Ui-2,j -
n/dU, A l,d30, A 3
+60(a?')yAl + '"
t,dU, A 8 d2U A 3
4(dz )ljA*+3(9a;3 )i,jAa:
32,350, A B
+ 60(5^)iiAa5 +"'
(4.58)
(4.59)
Now Richardson extrapolation yields
1
12A2
(-0i+2ii + 80i+1|j 8Ui-i'j + 8Ui-2,j)
,8U. r 1 d50 A 4
+ [ 30 dx5 As
Thus, for the streamwise convection term, we obtain
1
12AÂ£y,
-[(0 + 20o)t+2,j0+2,j + 8(0 + 20o)i+i,j0i+i1j
(4.60)
62

8(U + 2Uo)i-itjUi-itj + (U + 2Uo)i-2,jUi-2,j
ld(UU + 2U0U) 11 d\UU + 2U0U) i
Vr, dÂ£ kj + [ 30^ dp f ^ 1 }
Similarly, in the normal direction,
1 r (V + V0)i,j+2 .AV + V0)i>j+1
------r--------Ui^2 + 8:r
12A, yw ' yw ~"i+1
{V + Vo)i,i-l JJ {V + Vo)ij-2 TJ U0iij+1 Vitj+2
Vrij-i Vvi-1 2/tjj+j
U0i,j+1 Vi,3+1 ^Uoij-iVij-l +
2/^j+i Vvi-i Vvi-i
d ,uv + U0V + uv0^
a_v .. )i,j
-8-
i
+8-
dv Vr,
,r 1 d5 ,UV + UoV + UVo,A 4 , ,
+[-3W(----------S------)A +-k-
For the viscous terms, with the above Taylor series, we obtain
Uj+ij ~ 2Uj,j + Ui-i,i
A(2
,d2Us 1 ,d4U, aa2 1 ,d6Us aa4
(^2ki+12(^4)''.iA^ ^ 360 dp '3 +
(4.62)
Uj+2,i 2Uj,j + Ui-2,j
AAP
,d2U, 4 ,8*11 A,2 16 ,d6U. A ,4 ,
( 5Â£2) + 12( dp ^i,j ^ 360 3Â£6 i,J ^ +
(4.63)
Therefore, we have
-Ui+2tj + 16Ui+lij 30Utj + 10Uj -i ,j ffi-z.i
12AÂ£2
,d2Â£/\ ,r 156J7aMi ,
^ dÂ£2 ^ + * 90 dp + i,J'
(4.64)
----(_^1 + i6^l 30^ + 16^ ^t)
Vm 12Arl2 Vvi+2 VfU+l Vvi 3&W-1 Vrtj-2
63
Similarly,

ynjVyyj
1 ( ^.J + 2 | Q^j'+l
ft.*-:
!2At/ ^.+a 3/rjJ+1 2/,_
_ ^IhizL +
Ifru-a
= ^4(^)L+[-^S(^)V+---k"
For the pressure term,
_ l.W, A> 11,02P, a,2 11.53P, ...
piy pi d( )iJ + 2 22 ^ di2 k,i 6 23 dÂ£3 ^
1 1 d4P aa4 11 ,dBP. A,B
+ 04 o4 ( )*.i 1 on OB ( )*J AÂ£ ^ >
24 24Va^4
120 2B v 0Â£B
(4.65)
_ 3,dP, AA 1.3.,, d2P. A.2 l,3,3,5aP. A
P*\,j ~ Pi ^ 2( di ^2 2 3^2 ^62 di3
-l_-L(_wd ~P). .aÂ£4 -AÂ£b +
+ 24l2M 8? h'3 120 2 dÂ£B Jlj +
(4.66)
and, therefore,
-P,
V*i-
+U + 27P<+w 27pi-y + pi-1,
5P
- ^ )i.i + ( 64QTJ
24AÂ£
3 dBP
Vv
A Â£4+ )*.;
(4.67)
The modified partial differential equation corresponding to the scheme now can
be written as
au i djuu + 2Z70Z7) a cry + Er0y, + ffvr0
+ yv di + dri y
1 ,d2U 1 92 ,17, d U dP
ReS( <9Â£2 + y, cfy2(y) + Vr,Vvy dV{yj^ + Vv di
1 ^ At2 1 1 g5(^ + 21f0lf) 4
3 df3 30 y dÂ£B
1 95 + U0V + UV0
30 dif

2/tj
4 1 1 ,d*U. ,4 1 d6 4,
)A?1 + + -^(-)A71 )
PeS 90 di*
Vrjdjj6 yn
1 yrfyy 95 (_))An4 y A^4 + = 0
Re5 30 drf^yj* V y"m d? +
64
(4.68)

The terms in the last three lines of (4.68), involving the A quantities, represent
truncation error, denoted by R(AÂ£, A77, At). Evidently, if all unknown functions
and their combinations have at least six-order bounded derivatives, then
lim R(Ai,AV,At) = 0. (4.69)
AÂ£,At;,At>0
Thus, we find our scheme is consistent to the original governing equation.
We can also prove that the discrete Vmomentum and continuity equation
are also consistent to their original governing equations. The process is similar,
so not included here.
Stability
For the unsteady problem, numerical stability of the finite difference scheme is
the most important issue. For the problem of transition, since what we are going
to study is the physical instability, this issue becomes even more critical. Since
numerical instability and physical instability are usually confused in a numerical
calculation process, it is very difficult to identify them individually, and thus we
require the difference scheme be very stable over a quite wide range of regimes.
The governing equations are definitely nonlinear and strongly coupled, and so
exact stability analysis is almost impossible. Noting that the solver we used for
solving the discrete system is the linear iteration method, like Jacobi relaxation,
Gauss-Seidel relaxation, or SOR method, then linear analysis is still helpful,
although the original governing equations are still too complicated. Thus, in
65

this subsection, we will only give the analysis for a one-dimensional convection-
diffusion model problem.
The model equation can be expressed as
du du d2u
----h a-----b---= 0.
dt dx dx2
(4.70)
where a > 0, b > 0. The Von Neumann-Richtmeyer method is used for the
analysis. First, we study the stability property of a second-order backward Euler
and second-order central difference. The scheme can be written as
Ui+l ~U
3un+1 4un + un_1 -n+1
----------------1- a
2 At
.a?~ 2ri+ai
A2
n+l
t-1
2A
= 0.
(4.71)
Define
ai
a2
aAt
~Ax
bAt
Ax2
Then (4.71) becomes
3u"+1 4u" + u-1 + ~ "-i1)
-2a2+11 2<+1 + lift1) = 0. (4.72)
Now suppose there is some numerical disturbance with some special wavenumber
k introduced into our scheme. We denote this disturbance at time nAt and grid
point i as e" = Anetkxi, and the amplification factor by
Gi = c^V*.
66

After eliminating the original term, and ignoring the subscript i (to avoid the
confusion between the index and the imaginary part of the complex variables
used later), we obtain
3G 4 + G-1 + o1G(ee e~ie) 2 a2G(ei9 -2 + e~iB) = 0,
je -is
iS
-iS\
where, for simplicity, we omit the subscript i, and where 9 = kAx.
equation can finally be simplified to
Q
(3 + 8a2sin2 + i2a,\ sin6)G2 4 G +1 = 0,
a
where
Letting
G =
2 yjl 8a2sin2| i2a\ain9
3 + 8a2sin2| + i2a\sin6
E = 2 y 1 8a2sin2- i2aisin6,
F = 3 + 8a2ain + i2aiA..,
A
what we need to determine is when |G| < 1. Defining
p = ^(1 8a2ain2^)2 + 4 a2ain29,
. 1 8a2sin2%
cosp =--------------,
p
2a,isin9
we have
E = 2 (^)i,
67
(4.73)
The above
(4.74)

or
E = 2 (pe^+V).
We will show for the first case that the scheme is unconditionally stable. For the
second case, the proof is the same and is therefore omitted .
For the first case,
\E\2 = (2 pho3~)2 + psin2^
i /3
= 4 4p*cos~ + p.
z
On the other hand,
l*T
q
= (3 + 8 a2sin2-)2 + Aa\sin26
z
, ,e
= p2 + 8 + 64a2.si7i2-.
In order to satisfy
we need
\E\2 < |F|2,
i0 0
4 Ap*cos^- + p < p2 + 8 + 64a2ai7i2-,
z z
or equivalently
id 6
Ap*cos^- < p2 p + 4 + 64a2si7i2
Noting that the right-hand side can be rewritten as
(p ~ \f + Y + +64a2^ra,2^ >
68

we now only need to show
4p*cos < p2 p 4 + 64 Recalling that
~(3 1 + cos/3
cos =------------,
2 2
we have
4pacos^ = 4/jj .
1 1 8a2sz7i2Â£
^2+ 2-p
and, thus, we only need to show
1 1 8a2sm2Â§ , . . 20n2
) < (p p + 4 + 64a2szn2-)2,
16^2+ 2p
or
8/o + 8 64 z
+2(p2 p + 4)64a2sz7i2- + 642a\sinA-,
Z Z
since
(/>2 p + 4)2 8/j 8
= 2/>3 + 9/?2 16^> + 8
= (,-l)V + 8)
> o,
8p + 8 <(p2 p + 4)2,
69
we have

so we only need to show that
-64a2sin2^ < 2(p2 p + 4)64a2sira2^ + 642a2sin4^,
2
But this is true since

so RHS > 0 and LHS < 0. Therefore, this scheme is unconditionally stable.
For fourth-order central differences, the scheme is still unconditionally stable.
The scheme can be written as
3un+i 4un + un~x + Su^1 8ujix + ui+2
2At +0 12Ax
, -&? + Hugj1 30<+1 + 16<: <Â£
12Ax3
(4.75)
and the associated amplification equation becomes
d G a G 6
3G 4 + G-1 -|]r-i(8sin9 sin29) H(32sin2- 2sin29) = 0.
o u 2
We thus have
a2
.ai
[3 + -^(32sm2- 2sin29) + i-^(8sin9 ain29)]G2 4G + 1 = 0, (4.76)
and
G =
2 3 + ^-(32sin2| 2sin29) -f i^-(8ain9 sin29)
(4.77)
Note that
u2
,6
(32sin2- 2sin29)
n A n9 1 .
8a2{-sin sm 9)
o 2 ^ / 1 . A
= 8a29in -(- + sm -)
70