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

All rights reserved.

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,

Inc., for their helpful discussions.

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:

iao

u(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*)}

+e3d+Real{{cl>u, v, w)ad+eax+M}

+Â£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}

eu+Realifad+eW*-^}

^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}

eu+Realifad+eW*-^}

eu-Realifad-e*-132-^},

e2dReal{v2de^-a^-ut)}

etd+Realifad+e*-^13*-^}

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