7“ nuwycn

rr.y

5‘

I

v.

 

o. 043. .‘hl'a 7 41;"

.-' ‘ ﬁ-v'u

#53:; - 399135552”ij . ‘-
(2‘15": .51.?

'b'lA4
)j':

I"...

151;} \‘Irl‘lm
’33:; m. $29“
1":

I“: ‘<
‘3 A t,;2_:
..

"um

”i

1‘ r523... 31%; ’50}
“3’;”r...c 92w ;.5.

'25: 5
.. ' r‘ﬂj
.5335

52"“355'52
"(5' '55‘3’ '

«'55. .

'3; #52135?!“ -—n~

“$2.5,ou y“ ’ "
t' ‘ :II

'- l

~2».

ﬁrm

OLVI
we??? ‘
l
.

I ‘3'“;

‘
'12.".‘u
. .
8}" H"

51"?

”PM" ' .
H'-

. W?
J 2 WHY! ‘JMX

. 22' . h ’35. . 52-5:

'._v.;,."§m{h “7;” n ﬁr' 3’3“, ”i; ‘ﬁv‘ "12"", "W

”I“ 5;“. “. $741;
“"2

.55.
.1 $4 2.12.»...

. W...
L rid", ..

 

'...
I“$5.”;

PM“: '51-: 55- ‘ ﬁrm“

.. .. mu

4.
y’a-
WK :4.‘ ,t‘

W" 0‘" t“:
Luggagl‘lk» :3:
~91 We?“ ‘2

$.13, 35651.. 5%
‘7

vnvlul

"'",“:\' ¥

7.557

 

J Méé‘fl ,{L
W .3555; {5552.555

¢
,
~35”!

74325 ,_{:r&§,._

1“
3333143.“,«9 gifts-
'1;'

,. "r.
"'J i

5595;... ,
2U" w‘f‘i‘
owl”!

vESlURASRTYBRIE

1111 11111121111 111

 

 

 

 

 

111111111111111111111111111111

31293 0091

This is to certify that the
dissertation entitled

THE NUMERICAL SOLUTION OF CONFINED LAMINAR
FLOW PAST A MOVING BOUNDARY

presented by

Bashar Sulaiman AbdulNour

has been accepted towards fulﬁllment
of the requirements for

Doctor of Philosophy degree in Mechanical Ergneering

Maggie

Major professor

Merle C. Potter

MS U is an Afﬁrmative Action/Equal Opportunity Institution 0-12771

 

r LIBRARY
Michigan State
University

\_ I

T‘

 

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1F
11 11

 

 

 

 

 

 

 

 

 

 

 

 

ﬁlT—T

MSU I: An Affirmative Action/Equal Opportunity Institution
emulation”:

 

 

 

 

THE NUMERICAL SOLUTION OF CONFINED LAMINAR
FLOW PAST A MOVING BOUNDARY

By

Bashar Suleiman AbdulNour

A DISSERTATION
Submitted to
Michigan State University
in partial fulﬁllment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Mechanical Engineering

1990

ABSTRACT

THE NUMERICAL SOLUTION OF CONFINED LAMINAR
FLOW PAST A MOVING BOUNDARY

By

Bashar Suleiman AbdulNour

The solution of the time—dependent, viscous, incompressible, two-dimensional,
conﬁned flow past a moving boundary is considered. The problem is modeled as a
laminar ﬂow in the entry region of a rectangular channel with uniform inlet velocity

upstream and a fully-developed parabolic profile downstream.

A numerical solution is obtained from the continuity and Navier-Stokes
equations and the appropriate boundary conditions. It provides an exact solution over
the entire domain, in the sense that it utilizes all terms in the describing equations.
The vorticity-stream function formulation is selected for the treatment of the problem.
After discretizing the non-dimensional governing equations in finite-difference form,
the Altemating-Direction Implicit method is used to solve the vorticity equation while
the iterative Successive Over-Relaxation method is used to solve the stream function
equation. The vorticity wall condition is derived by expanding in a Taylor series. This
second-order scheme maintains the diagonal dominance of the system of equations at

high Reynolds number, thereby expanding the domain of convergence.

The vorticity, stream function, and velocity components are calculated at each
grid point for Reynolds numbers ranging between 10 and 5000. The plane channel
ﬂow is first investigated. This steady-ﬂow solution is obtained as the asymptotic
solution for large time of the unsteady ﬂow. The results conﬁrm two primary
conclusions: the local maxima in the streamwise velocity profile near the wall is found
to be minimal, and the inviscid-core region is approximately one-ﬁfth of the entrance
length with a proﬁle-development region accounting for the rest. The use of
accelerating parameters to control the rate of convergence has signiﬁcantly improved

numerical stability and reduced the time of computation.

When a ﬁnite step is ﬁxed to the lower boundary, however, distinct separated
regions with pronounced eddies of recirculating fluid are developed on both sides.
Motion of part of the smooth lower boundary is then included; both steady and
oscillatory motion are studied. The size of a separated region, formed when the
boundary is moving opposite to the flow, depends on the Reynolds number, Strouhal
number, and velocity of the boundary.

Cepyright by
Bashar Sulaiman AbdulNour
1990

To my father Sulaiman and my mother Lydia

for their inﬁnite love and unequivocal support.

ACKNOWLEDGMENTS

Foremost, the author wishes to express his enormous gratitude and
indebtedness to his major professor, Dr. Merle C. Potter. Mere words cannot convey
this since Merle was not only an academic advisor, but also a true friend and a
signiﬁcant inﬂuence. Most importantly, the author greatly appreciates the remarkable
interest and the invaluable educational experience gained from his association with

Professor Potter.

The author is also much obliged to his guidance committee members Professor
James V. Beck and Professor David H. Y. Yen for their time and contributions. In
addition, sincere thanks are extended to Dr. Craig W. Somerton for his constructive

comments and to Dr. Reinier B. J. Bouwmeester for his early participation.

Above all, without the strong encouragement and understanding that the author
received from his family during his college career, this effort would have not been
possible. Finally, the author’s faith in the Lord Jesus Christ was a primary source of

inspiration.

TABLE OF CONTENTS

LIST OF TABLES ............................................................................................................ ix
LIST OF FIGURES ........................................................................................................... x
NOMENCLATURE ....................................................................................................... xvi

1. INTRODUCTION .................................................................................................... l

2. BACKGROUND ...................................................................................................... 8

2.1 Entry Flow in a Pipe ........................................................................................ 9

2.2 Entry Flow in a Channel ................................................................................ 14

2.3 Internally Separated Flows ............................................................................ 20

2.4 Flow in a Driven Cavity ................................................................................. 23

2.5 Externally Separated Flows ........................................................................... 27

2.6 Time-Dependent Flows .................................................................................. 28

2.7 Numerical Treatments of Flow Problems ...................................................... 3O

3. DESCRIPTION OF THE PROBLEMS ................................................................. 32

3.1 The Flow Problems ........................................................................................ 32

3.2 Mathematical Formulation ............................................................................. 35

3.3 The Problem for Solution: Vorticity-Stream Function Formulation ............. 39

3.4 The Steady Flow Problem .............................................................................. 42

3.5 Initial and Boundary Conditions .................................................................... 43

4. THE NUMERICAL METHODS ............................................................................ 51
4.1 The Computational Domain ........................................................................... 52

4.2 The Alternating-Direction Implicit Method for the Vorticity Equation ....... 52

4.3 The Successive Over-Relaxation Method for the Stream Function

Equation ......................................................................................................... 61

4.4 Initial and Boundary Conditions .................................................................... 67

5. THE COMPUTATIONAL SCHEME .................................................................... 72
5.1 Accelerating Parameters ................................................................................ 72

5.2 Outline of the Computational Procedure ....................................................... 75

5.3 Computing Environment ................................................................................ 80

6. ANALYSIS ............................................................................................................. 81
6.1 Error Estimate for Vorticity Boundary Conditions ........................................ 81

6.2 Consistency and Errors of a Numerical Solution ........................................... 83

6.3 Numerical Stability of the Differenced Navier- Stokes Equations ................ 85

7. DISCUSSION OF RESULTS AND CONCLUSIONS .......................................... 89
7.1 Entry Flow in a Plane Channel ...................................................................... 89

7.2 Channel Flow with a Finite Step .................................................................... 95

7.3 Channel Flow with a Moving Boundary ........................................................ 98

7.4 Channel Flow with an Oscillating Boundary ............................................... 101

7.5 Accelerating Parameters .............................................................................. 104

7.6 Conclusions .................................................................................................. 108
APPENDIX A - COMPUTER PROGRAM .................................................................. 170
BIBLIOGRAPHY .......................................................................................................... 188

viii

Table 1.

Table 2.

Table 3.

Table 4.

Table 5.

LIST OF TABLES

The inviscid-core length Li, the proﬁle-development length L d,
and the entrance length L e for various Reynolds numbers. ........................ 91

Comparison of the ratio L¢/(I-I-Re) as reported in different
channel studies for Reynolds number = 200. .............................................. 92

Length of the aft separated region L2 for the channel ﬂow with a
ﬁnite step. .................................................................................................... 96

Optimum accelerating parameters, CPU time, and number of outer
iterations for the channel ﬂow with a uniform mesh of size 0.05 at
various Reynolds numbers; the solution is based on one half the
height of the channel. ................................................................................ 105

Optimum accelerating parameters, CPU time, and number of outer
iterations for the channel ﬂow with a uniform mesh of size 0.05 at
various Reynolds numbers; the solution is based on the full height
of the channel. ........................................................................................... 106

Figure 1.

Figure 2.
Figure 3.
Figure 4.
Figure 5.
Figure 6.
Figure 7.

Figure 8.

Figure 9.

Figure 10.

Figure 11.

Figure 12.

LIST OF FIGURES

A schematic representation of the flow problem including the

moving boundary. ....................................................................................... 33
Laminar entry ﬂow in a two-dimensional channel. .................................... 34
Schematic representation of the sequence of entry ﬂow problems. ............ 36
Dimensionless boundary conditions. .......................................................... 48
The ﬂow problem with the ﬁnite step placed on the lower boundary. ....... 49
Fume-difference representation of the computational domain. .................. 53
Flow chart of the computational procedure. ............................................... 79

Development of streamwise velocity proﬁles through the entry
region of a two-dimensional channel at various Reynolds numbers. ....... 111

Streamwise velocity proﬁles at locations near the inlet section for
various Reynolds numbers. ....................................................................... 114

Development of streamwise velocity proﬁles through the entry
region of the channel using the boundary-layer equations at Re =
200. ........................................................................................................... 115

Development of streamwise velocity proﬁles through the entry
region of the channel using a cascade inlet condition at Re = 200. .......... 115

Development of streamwise velocity proﬁles through the entry
region of the channel using an actual experimental inlet velocity
proﬁle at Re = 200. ................................................................................... 115

X

Figure 13.

Figure 14.

Figure 15.

Figure 16.

Figure 17.

Figure 18.

Figure 19.

Figure 20.

Figure 21.

Figure 22.

Figure 23.

Figure 24.

Figure 25.

Figure 26.

Figure 27.

Streamwise velocity development along the channel at different
transverse locations across the channel for Re = 200 and 2000. .............. 116
Centerline velocity along the channel for various Reynolds numbers. 117
Development of normal velocity proﬁles through the entry region

of a two-dimensional channel at various Reynolds numbers. ................... 118

Development of vorticity proﬁles through the any region of a
two-dimensional channel at various Reynolds numbers. .......................... 119

Vorticity at the upper boundary along the channel at various

Reynolds numbers. .................................................................................... 120
Vorticity at the upper boundary along the channel calculated from
different formulas for Re = 200. ............................................................... 120
Vorticity contours of the ﬂow in the entry region of a two-
dimensional channel at Re = 200. ............................................................. 121
Streamlines of the flow in the entry region of a two-dimensional
channel at Re = 200. ................................................................................. 121
Normalized pressure gradients at the wall along the channel for
various Reynolds numbers. ....................................................................... 122
Normalized pressure gradients at the centerline along the channel
for various Reynolds numbers. ................................................................. 122
Normalized pressure gradients at the wall and centerline along
the channel for Re = 200. .......................................................................... 123
Streamwise velocity as a function of the mesh size at selected
locations; Re = 200. .................................................................................. 124

Flow separation due to a ﬁnite step of various width; X1 = 1.0, B
= 0.25, Re = 20. ......................................................................................... 125

Flow separation due to a finite step of various height; X1 = 1.0, D
= 0.05, Re = 20. ......................................................................................... 127

Flow separation due to a ﬁnite step of various width; X1 = 2.0, B
= 0.20, Re == 200. ....................................................................................... 129

Figure 28.

Figure 29.

Figme 30.

Figure 31.

Figure 32.

Figure 33.

Figure 34.

Figure 35.

Figure 36.

Figure 37.

Figure 38.

Figure 39.

Figure 40.

Flow separation due to a ﬁnite step; Re = 200, X1 = 2.0, B = 0.25,
D = 0.05. ................................................................................................... 130

Flow separation due to a ﬁnite step of various height; X1 = 1.0, D
= 0.05, Re = 200. ....................................................................................... 130

Effect of step location on the ﬂow separation due to a finite step;
X1=8.0,B=0.20,D=0.05,Re=200. ................................................... 131

Velocity vectors of the flow in the presence of flow separation due
toaﬁnitestep;Xl=1.0,B=0.20,D=0.05,Re=20. ............................ 131

Development of streamwise velocity proﬁles along the channel in
the presence ofafinite step;Xl = 1.0,Re=20. ....................................... 132

Development of streamwise velocity proﬁles along the channel in
the presence of a ﬁnite step; X 1 = 2.0, B = 0.2, Re = 200. ....................... 133

Development of streamwise velocity proﬁles along the channel in
the presence of a ﬁnite step; X1 = 1.0, B = 0.25, D = 0.05, Re 2

200. ........................................................................................................... 134

Vorticity contours in the presence of ﬂow separation due to a
ﬁnite step; X1 = 1.0, Re = 20. ................................................................... 134

Streamlines of the ﬂow with the lower boundary section moving
opposite to the inlet velocity at UB; X1 = 1.0, X2 = 2.0, Re = 20. ........... 135

Velocity vectors of the flow with the lower boundary section
moving opposite to the inlet velocity at UB = -10.0; X1 = 1.0, X2 =

2.0, Re = 20. .............................................................................................. 136

Streamlines of the ﬂow with the lower boundary section moving in
the direction of the inlet velocity at UB; X1 = 1.0, X2 = 2.0, Re =
20. ............................................................................................................. 137

Streamlines of the ﬂow with the lower boundary section moving
opposite to the inlet velocity at UB; X1 = 1.0, X2 = 2.0, Re = 200. ......... 138

Streamlines of the ﬂow with the lower boundary section moving in
the direction of the inlet velocity at UB = 1.0; X 1 = 1.0, X2 = 2.0,

Re = 200. ................................................................................................... 139

Figure 41.

Figure 42.

Figure 43.

Figure 44.

Figure 45.

Figure 46.

Figure 47.

Figure 48.

Figure 49.

Figure 50.

Figure 51.

Figure 52.

Streamlines of the ﬂow with the lower boundary section moving
opposite to the inlet velocity at UB = -1.0; X1 = 1.0, X2 = 5.0, Re

= 200. ........................................................................................................ 139

Streamlines of the ﬂow with the lower boundary section moving
opposite to the inlet velocity at UB; X1 = 4.0, X2 = 5.0, Re = 200. ......... 140

Streamlines of the ﬂow with the lower boundary section moving
opposite to the inlet velocity at UB; X1 = 4.0, X2 = 5.0, Re = 200.

(A close up of Frg. 42.) ............................................................................. 141

Streamlines of the ﬂow with the lower boundary section moving at
U13;Xl =2.0,X2=4.0,Re=200. ........................................................... 142

Streamlines of the ﬂow with the lower boundary section moving at
UB; X1 = 2.0, X2 = 6.0, Re = 200. ........................................................... 143

Streamlines of the flow with the lower boundary section moving at
UB; X1 = 8.0, X2 = 10.0, Re = 200. ......................................................... 144

Streamlines of the ﬂow with the lower boundary section moving
opposite to the inlet velocity at UB = -0.5; X1 = 10.0, X2 = 14.0,

Re = 1000. ................................................................................................. 145

Streamlines of the flow with the entire lower boundary moving in
the direction of the inlet velocity (Couette-Poiseuille ﬂow) at UB

= 1.0; Re = 500. ........................................................................................ 145

Development of streamwise velocity proﬁles along the channel
with the lower boundary section moving at UB; X1 = 2.0, X2 = 4.0,

Re a 200. ................................................................................................... 146

Development of streamwise velocity proﬁles along the channel
with the lower boundary section moving at UB; X1 = 2.0, X2 = 6.0,

Re a 200. ................................................................................................... 147

Development of streamwise velocity proﬁles along the channel
with the lower boundary section moving opposite to the inlet
velocity at UB = -0.5; X1 = 10.0, X2 = 14.0, Re = 1000. ......................... 148

Development of streamwise velocity proﬁles along the channel
with the entire lower boundary moving in the direction of the inlet
velocity (Couette-Poiseuille ﬂow) at UB . ................................................. 148

xiii

Figure 53.

Figure 54.

Figure 55.

Figure 56.

Figure 57.

Figure 58.

Figure 59.

Figure 60.

Figure 61.

Figure 62.

Figure 63.

Figure 64.

Vorticity contours with the lower boundary section moving at UB;
Xl =1.0, X2 = 2.0, Re = 20. ..................................................................... 149

Vorticity contours with the lower boundary section moving at UB;
X1=1.0,X2=2.0,Re=200. ................................................................... 150

Vorticity contours for the entire lower boundary moving in the
direction of the inlet velocity (Couette-Poiseuille ﬂow) at UB =

1.0; Re = 500. ............................................................................................ 150

Streamlines of the ﬂow at different time with the lower boundary
section oscillating at UB = Acos(t); X1 = 1.0, X2 = 2.0, Re = 20,

St= l.0,A= l.0,n=8. ............................................................................. 151

Streamlines of the ﬂow at different time with the lower boundary
section oscillating at UB = Acos(t); X1 = 2.0, X2 = 4.0, Re = 200,
St= 1.0,A=2.0,n=25. ........................................................................... 152

Streamlines of the flow at different time with the lower boundary
section oscillating at UB = Acos(t); X1 = 2.0, X2 = 6.0, Re = 200,

St = 0.5, A = 2.0, n = 14. ........................................................................... 153

Streamlines of the flow at different time with the lower boundary
section oscillating at UB = Acos(t); X1 = 1.0, X2 = 2.0, Re = 200,

St= 1.0,A= l.0,n=21. ........................................................................... 154

Size of the separated region as a function of the Strouhal number;
X1 =1.0, X2 = 2.0, Re = 200, A = 1.0, t = (2n+1)n, UB = -1.0. ............. 155

Size of the separated region as a function of the Strouhal number;
X1 = 2.0, X2 = 4.0, Re = 200, A = 1.0, t = (2n+1)rt, UB = -1.0. ............. 156

Size of the separated region as a function of the Strouhal number;
X1: 2.0, X2 = 6.0, Re = 200, A = 0.5, t = (2n+ l)rt, UB = -0.5. ............. 157

Flow sueamlines showing the size of the separated region as a
function of the amplitude of oscillation of the boundary; X1 = 1.0,

x2 = 2.0, Re = 20, St = 2.0,t=(2n+1)1t. ................................................. 158

Size of the separated region as a function of the amplitude of
oscillation of the boundary; X1 = 1.0, X2 = 2.0, Re = 200, St = 1.0,

t= (2n+l)rt. .............................................................................................. 159

Figure 65.

Figure 66.

Figure 67.

Figure 68.

Figure 69.

Figure 70.

Figure 71.

Figure 72.

Figure 73.

Figure 74.

Figure 75.

Size of the separated region as a function of the amplitude of
oscillation of the boundary; X1 = 2.0, X2 = 4.0, Re = 200, St = 2.0,

t = (2n+ l)x. .............................................................................................. 160

Size of the separated region as a function of the amplitude of
oscillation of the boundary; X1 = 2.0, X2 = 6.0, Re = 200, St = 1.0,

t= (2n+l)1t. .............................................................................................. 161

Streamlines of the ﬂow with the lower boundary section
oscillating at UB = Acos(t); Re = 1000, A = 1.0, t = (2n+l)1t, UB
= -1.0. ........................................................................................................ 162

Effect of the location of the oscillating section of the boundary on
the size of the separated region; X1 = 5.0, X2 = 6.0, Re = 200, St

=1.0,A= l.0,t=(2n+l)1t,n=28,UB =-l.0. ........................................ 163

Effect of the length of the oscillating section of the boundary on
the size of the separated region; X1 = 1.0, X2 = 3.0, Re = 200, St

=1.0,A=1.0,t=-=(2n+1)rt,n=26,UB =-1.o. ........................................ 163

Effect of the location and length of the oscillating section of the
boundary on the size of the separated region; X1 = 8.0, X2 = 10.0,

Re=200,t=(2n+l)x. ............................................................................. 164

Development of streamwise velocity proﬁles along the channel at
the instant when the boundary is moving at maximum velocity in

the direction of the ﬂow; Re = 200, t = 2n1t. ............................................ 165

Development of streamwise velocity proﬁles along the channel at
the instant when the boundary is moving at maximum velocity in

the direction opposite to the flow; Re = 200, t = (2n+1)1t. ...................... 166

Vorticity contours for the case of the lower boundary section
oscillating at UB = Acos(t); X1 = 1.0, X2 = 2.0, Re = 20, St = 1.0,

A= l.0,n=8. ........................................................................................... 167

Vorticity contours for the case of the lower boundary section
oscillating at UB = Acos(t); X1 = 1.0, X2 = 2.0, Re = 200, St =
1.0,A=1.0,n=21. .................................................................................. 168
Vorticity contours for the case of the lower boundary section
oscillating at UB = Acos(t); X1 = 8.0, X2 = 10.0, Re = 200, St =
l.0,A=1.0,n=24. .................................................................................. 169

XV

gvazz_rargrwrmmauw>

NOMENCLATURE

amplitude of oscillation of boundary velocity

height of the ﬁnite step

length of the ﬁnite step

over-relaxation factor for the stream function equation
gravitational acceleration

channel height

mesh size (equal in the x- and y-directions)

time step

proﬁle-development length

entrance length

inviscid-core length

number of grid points in the x-direction

number of grid points in the y-direction

number of cycles of oscillation for convergence
pressure

kinetic pressure

Reynolds number

Strouhal number

time

uniform velocity at the inlet

velocity of the moving section of the boundary
streamwise velocity component (in x-direction)
normal velocity component (in y-direction)

velocity component in the z-direction

upstream location of the oscillating boundary or the step
downstream location of the oscillating boundary or the step

streamwise coordinate

xvi

y = normal coordinate

z = coordinate perpendicular to the xy-plane
Greek symbols

or = f'mite-difference operator for the x-direction
[3 - ﬁnite-difference operator for the y-direction
y weighting factor for the stream function

A incremental increase (in space or time)

8 tolerance error

i; z-component vorticity

v kinematic viscosity

p density of the ﬂuid

0' = weighting factor for the vorticity

\l’ stream function

or = frequency of oscillation of the boundary
Vector symbols

i unit vector in the x-direction

i = unit vector in the y-direction

k = unit vector in the z-direction

V = velocity vector

x position vector

Subscripts

B = nodes on the lower or upper boundaries

CL = centerline nodes

i = x-location index of an interior node

j = y—location index of an interior node

M = row of nodes on the upper boundary

N = column of nodes on the downstream boundary
Superscripts

m = index of inner iterations for the stream function
k = index of time step

(") = relaxed value

(') = weighted value

( )* = dimensional variable

xvii

INTRODUCTION

Exact analytical solutions are available only for a limited number of ﬂuid
dynamics problems. In general, simplifying assumptions are applied to the describing
equations of motion to facilitate a solution. Unfortunately, for most complicated ﬂow
problems such assumptions are not suitable. Consequently, these ﬂows can be
approached experimentally or by solving the equations of motion either numerically or

by means of approximate analytical techniques.

In complex ﬂow problems it is necessary to determine the detailed velocity and
pressure ﬁelds before an integral quantity of interest, such as an energy requirement
or a wall shear force, can be calculated. From such calculations we can then predict
how the integral quantity is inﬂuenced by the ﬂow ﬁeld. One fundamental way to

achieve this task is to obtain a numerical solution to the mathematical model of the

physical problem at hand.

Advancements in high speed digital computing have resulted in a widespread
interest in numerical solutions of the Navier-Stokes equations. Harlow and From
(1965) remarked that computer calculations can give rapid and inexpensive, as well

as accurate, solutions and that it is possible to perform sophisticated experiments on

1

2
the computer alone. Recent developments and advances in computational ﬂuid
dynamics were reviewed by Krause (1985). A comprehensive insight into the
interdisciplinary nature of computational ﬂuid dynamics and its application to carry out

numerical experimentation was also explored by Aref (1986).

The ﬁnite-difference approximation is one among a number of numerical methods
that have been developed and studied extensively. In principle, furite-difference
solutions can be obtained for ﬂows having arbitrary Reynolds numbers. At high
Reynolds numbers, however, convergence of the numerical solution depends on the
numerical method, spatial grid resolution, and size of the time step. Difﬁculties in
attaining convergence is usually associated with the loss of diagonal dominance of the

system of difference equations approximating the ﬂow ﬁeld.

The accuracy of the numerical solution is also dependent on the grid size and the
error associated with the numerical approximation and the computation. For the
majority of practical ﬂows involving relatively high Reynolds numbers, the ﬂow ﬁeld
must be represented adequately in order for a numerical solution to be sufﬁciently
accurate. Consequently, the number of nodal points may become prohibitively large
from the practical standpoint of computer time and storage. In addition, the lack of
uniqueness in the numerical solution, due to the nonlinearities of the Navier-Stokes
equations or inappropriate numerical representation of the boundary conditions, may

prevent meaningful results.

The present study is motivated by the potential importance of specifying
solutions to ﬂow problems comprising moving boundaries as well as investigating

ﬁnite-difference methods effective in analyzing hi gh-Reynolds-number laminar ﬂows.

3
Flow ﬁelds containing a moving boundary have many important applications in
unsteady ﬂuid mechanics. These applications include oscillating airfoils, helicopter
rotors, impellers, compressor blades, bafﬂe conﬁgurations, vessel geometries,
washing machine agitators, agitated tanks, and industrial mixers. The ﬂow patterns
in such flows are usually complicated and affected by the ﬂuid, dynamics of the ﬂow,
ﬂow ﬁeld geometry, and motion and shape of the moving boundary.

Most importantly, moving boundaries can produce flow separation. The
phenomenon of separation is one of the important features of ﬂuid ﬂow. It is an old
subject in the literature of ﬂuid mechanics. Using ﬂow visualization techniques,
Leonardo da Vinci observed and sketched recirculating eddies formed in the ﬂow over
objects with various conﬁgurations. Still, a thorough understanding of the subject of

ﬂow separation has not been obtained despite this long history.

Laminar ﬂow separation and reattachment of various kinds occur frequently in
complex ﬂows. One speciﬁc example is the laminar separated regions occurring in
regimes of low-Reynolds-number ﬂow around airfoils. Another example is the
separated bubble generated by the ﬂow near the leading edge of a thin airfoil.
Unpredictable rapid growth of these separated regions can cause stall. The problem
of ﬂow separation is rather general, however, but it is the separation problem due to a
moving boundary that is being addressed here. This will be illustrated by a steady
and an unsteady motion of the boundary.

The effect of boundary motion on a two-dimensional, viscous, incompressible,
conﬁned ﬂow is considered in this study. The problem is modeled as a time-

dependent, isothermal, laminar ﬂow of a Newtonian ﬂuid past a moving boundary

4
section placed within the entry region of a semi-inﬁnite, straight, rectangular channel.
The channel has a uniform velocity at the inlet and fully-developed Poiseuille ﬂow far
downstream. Both steady and oscillatory boundary motions are investigated. Our
main objective is to numerically simulate the ﬂuid ﬂow when a moving boundary
exists thereby developing a numerical technique suitable to treat problems that

contain time-dependent boundary conditions.

The solution can be obtained numerically by considering the continuity and the
Navier-Stokes equations along with the appropriate boundary conditions. The full
Navier—Stokes equations do not invoke the boundary-layer assumptions and,
therefore, represent an exact solution in the sense that no terms in the x- and y-
momentum equations are neglected. The vorticity-stream function formulation is
selected for the treatment of the problem. The ﬁnite-difference equations are
obtained from the partial differential equations of motion using second-order
discretizations. After expressing the initial and boundary conditions in terms of the
stream function and the vorticity in ﬁnite-difference form, a solution to the system of
equations can be pursued to yield the vorticity and the stream function. From this

information the velocity ﬁeld, which is of primary interest, can be calculated.

The Alternating-Direction Implicit (ADI) method is used to solve the vorticity
equation while the iterative Successive Over-Relaxation (SOR) method is used to
solve the stream function equation. Finite-difference approximations used when
modeling physical problems often lead to the solution of large-order sparse linear
systems. For large Reynolds numbers, however, it is essential to ensure the stability
of the numerical solution by maintaining the diagonal dominance of the system of

equations. An upwind differencing applied to the convective terms in the vorticity

5
equation helps avoid numerical instability. Second-order accuracy is then recovered
at convergence via a correction technique. The power of this ADI/SOR computational
scheme lies in the structure of the difference equations, which for all Reynolds
numbers, yields a diagonally dominant system of linear algebraic equations. The
present numerical solution is second-order accurate in its entirety since both the

difference equations and boundary conditions have second-order accuracy.

The use of accelerating parameters can further improve the convergence
characteristics of the numerical procedure. Although somewhat artiﬁcially inserted
into the difference equations, such parameters expand the domain of convergence and
improve the computational speed as well. The general notion of accelerating the
solution is well known. Nevertheless, a systematic determination of the optimum
values of the accelerating parameters suitable for the nonlinear set of simultaneous
equations resulting from the discretization of the Navier-Stokes equations has not
been established. Therefore, a search must be made for the optimum values of the
accelerating parameters in order to minimize the computing time and identify the
domain of convergence. The optimum values of these parameters depend on the
shape of the computational domain, mesh size, and the Reynolds number. They are

determined primarily by numerical experimentation.

The development of the steady ﬂow in the entry region of a two-dimensional
channel is ﬁrst investigated using the full Navier-Stokes equations. Conduit ﬂows
are important to many engineering and biological applications, such as piping, heat
exchanger design, short passages leading to nozzles and diffusers, and blood ﬂow in
arterial systems. Furthermore, the entry ﬂow problem provides a prelude for

developing laminar ﬂow dreary and serves as an excellent example for testing

6
numerical schemes for solving the equations of motion. The streamwise velocity
proﬁle undergoes a change from an assumed inlet proﬁle to that of a fully—developed,
parabolic proﬁle at a location downstream. The two regions making up the entrance
ﬂow are identiﬁed; the inviscid-core region, in which a viscous layer is assumed to
exist near the wall, is followed by the proﬁle-development region, in which viscous
effects inﬂuence the entire channel ﬂow. The results conﬁrm two primary conclusions:
the overshoot in the streamwise velocity proﬁle near the wall is found to be minimal,
and the inviscid-core region is approximately one-ﬁfth of the entrance length with the

proﬁle-development region accounting for the rest.

A second problem of entry ﬂow in a channel constricted with a step is also
discussed. When a ﬁnite step is ﬁxed to the lower boundary, distinct separated

regions with pronounced eddies of recirculating ﬂuid are developed on both sides.

A third class of problems involves an extension of the ﬁrst; the channel ﬂow is
investigated when part of the lower boundary is allowed to move in its own plane,
with either a constant or an oscillating velocity. Interest is focused on the problem of
boundary motion that leads to the formation of a separated region during the upstream
motion of a section of the lower boundary; the effect of the velocity, length, and
location of the moving boundary on the size of the separated region is of particular

interest.

Prediction of such separated ﬂows was previously done by patching the
numerical solution in the vicinity of the separated region to the reminder of the ﬂow
ﬁeld. Advancements in computing power make it possible to obtain a numerical

solution for the entire ﬂow ﬁeld, thus avoiding difﬁculties associated with patching the

7
local numerical solution to the reminder of the ﬂow. Moreover, a complete numerical
solution is attractive for ﬂows in which the boundary—layer theory is questionable for
the range of Reynolds numbers considered, as is the case with separated ﬂows,
where the boundary-layer equations do not properly describe the separated portion of
the ﬂow. The present scheme provides a numerical solution over the entire domain

utilizing all terms in the Navier-Stokes equations.

It is the purpose of this study to also use different approximations in expressing
the wall boundary conditions on the vorticity, examine several techniques for
evaluating the coefﬁcients of the convective terms in the vorticity equation,
investigate the inﬂuence of several terms in the describing equations of motion, and
study the effect of the mesh size and Reynolds number on the numerical solution.

Numerical solutions are obtained for Reynolds numbers ranging ﬁom 10 to 5000.

Furthermore, the over-relaxation and weighting factors for the present
computational domain are determined. Their use in controlling the rate of convergence
has signiﬁcantly improved the numerical stability. The computer time is also reduced

by up to a factor of three using the optimum values of these accelerating parameters.

Reviews of work pertinent to this study are highlighted in Chapter 2. The
present ﬂow problems and their mathematical formulation are outlined in Chapter 3.
Chapter 4 contains the details of the numerical methods used while Chapter 5
illustrates the mechanics of the computational scheme. Important considerations in
the solution of the Navier-Stokes equations are discussed in Chapter 6. Finally, the

results, along with a discussion and conclusions, are presented in Chapter 7.

BACKGROUND

In this chapter, it is our intention to highlight works relevant to the ﬂow problem
considered in this numerical study. First, a historical perspective into the laminar
entry ﬂows in pipes and channels is given. The ﬂow in the entry region of pipes and
channels is a classic problem in ﬂuid mechanics and one that has important
applications as well. In general, entry ﬂows are studied by analytical methods, which
include integral methods, series expansions, and linearization techniques; numerical
methods; and experimental methods. Most of the analytical and numerical analyses
involve some form of boundary-layer approximation in which a boundary layer exists
near the walls together with an inviscid flow at the core near the inlet section.
Bibliographies of laminar entry ﬂows in pipes and channels, classiﬁed according to the
method of solution, are given by Schmidt and Zeldin (1969), Lew and Fung (1970),

Van Dyke (1970), and Fargie and Martin (1971).

Secondly, a variety of numerical studies, in ﬁnite-difference form, pertinent to
the ﬂow problems of this study are presented. In particular, numerical considerations
of internal and external separated ﬂows, in addition to time-dependent ﬂows of

similar nature to the present investigation and some special numerical treatments to

9
ﬂow problems, are discussed. Although numerical solutions to ﬂow problems
somewhat similar to parts of the present study, with some using correlated numerical
techniques, are presented in the literature, no work could be identiﬁed with direct

application to ﬂows with moving boundaries or time-dependent boundary conditions.

A brief survey of the speciﬁc numerical methods used in this study will be

presented along with the numerical methods in Chapter 4.

2.1 Entry Flow in a Pipe

Due to its resemblance in the physical setting and mathematical approach to the
entry ﬂow problem in a channel, the analogous problem in a circular pipe is included in
this review. A variety of approximate analytical and numerical methods as well as
experimental studies have been employed for the determination of the ﬂow
characteristics in the entry region of a pipe, as reported in the literature. Although not
strictly considered a boundary-layer problem (Schliehting, 1979), this axisymmetric
ﬂow has been studied with the aid of methods similar to boundary-layer solutions.

An assumed inﬁnitesimally thin boundary layer on the wall of a pipe at the inlet
section, grows in thickness as the ﬂow passes downstream until it becomes equal to
the radius of the pipe. The ﬂuid at the core is practically uninﬂuenced by the
viscosity. The fully-developed parabolic velocity distribution is approached
asymptotically at a distance well beyond when the boundary layer thickness becomes
equal to the pipe radius. At a large distance downstream from the inlet, the velocity

distribution approaches the theoretical parabolic Poiseuille ﬂow over the cross-

section of the pipe.

10

Prandtl and Tietjens (1934) presented a summary of early analyses of the ﬂow
in pipes and channels. Independently, Hagen and Poiseuille (Prandtl and Tietjens,
1934) were the ﬁrst to experimentally investigate the ﬂow in straight tubes and
include the inﬂuence of viscosity. Primarily, they were concerned with the application
of the equations of hydrodynamics to estimate the pressure drop. Their ﬁndings,
however, apply only to ﬂow at a sufﬁcient distance from the inlet section. On the
other hand, using a perturbation of the parabolic distribution, Boussinesq (1891) was
the ﬁrst to make theoretical studies of the phenomenon, but his calculations for the

velocity distribution did not agree with the experimental results.

In the entry region, it is clear that no parabolic distribution exists. Prandtl
(Prandtl and Tietjens, 1934) suggested a velocity proﬁle approximated by two
parabolic curves near the walls, with zero velocities exactly at the walls, joined by a
ﬂat middle part near the pipe axis. At the inlet section, the width of the parabolic
curves is zero. It then increased as the ﬂow advances downstream until the two
curves are united into a single parabola. The ﬂow in the uniform velocity core was

assumed to satisfy Bemoulli’s equation.

Schiller (1922) was the ﬁrst to use an approximate method based on an
equilibrium condition between momentum, pressure drop, and viscous drag to form an
equation similar to the momentum integral equation. Schiller dropped the viscous
terms and employed Bernoulli’s equation in his analysis to substitute for the pressure
term in his integration of the Navier-Stokes equations. He designated the entrance
length to be the streamwise location where the two parabolic curves are combined.
Carrrpbell and Slattery (1963) corrected the original derivation of Schiller (1922) to

account for the loss of energy due to viscous dissipation within the boundary layer.

ll
Employing both integral and differential momentum equations, Fargie and Martin
(1971) simpliﬁed the treatment of the previous authors. They presented a thorough
review of previous studies and compared their results as well.

Lately, Mohanty and Astlrana (197 8) investigated the laminar ﬂow in the entry
region of a pipe using an integral method. They integrated the boundary-layer
equations in the region closer to the inlet section according to the von KArmén-
Pohlhausen scheme. Using fourth-degree polynomial approximations to the velocity
proﬁles, they solved the Navier-Stokes equations with an "order-of-magnitude"
analysis in the region further away from the inlet. They also showed that the laminar
viscous layers meet at the pipe axis, when the inviscid core terminates, at a distance
of approximately one-quarter of the entrance length. The previous literature did not
indicate the location of such a meeting of the viscous layers in an entry ﬂow, although

the proﬁles of Schlichting (1979) appear to give such an impression.

The unpublished work of Atkinson and Goldstein for the ﬂow near the inlet of a
pipe, presented by Goldstein (1965), used a variation of the method due to Schlichting
(1934) for the corresponding problem in a two-dimensional channel. From the
equations of motion in cylindrical coordinates, a solution is obtained for large
Reynolds number by generalizing Blasius’ solution of Prandtl’s two-dimensional
boundary-layer equation. A numerical integration of the equation resulting from
applying a power series expansion yields the velocity distribution. Results agreed
with the calculations of Schiller (1922) except near the wall. Away from the inlet,
Atkinson and Goldstein continued the solution by calculations similar to those of
Boussinesq (1891). The combination of these two solutions produced more accurate

results than those of Boussinesq (1891) and Schiller (1922).

12
Langhaar (1942) analyzed the ﬂow in the why region of a straight tube by
means of a linearizing approximation to the Navier-Stokes equations. He deﬁned the
velocity proﬁles by Bessel functions and obtained the entrance length. He also
determined the pressure drop in the entry region by applying the energy equation and
the computed velocity ﬁeld. Throughout his work, the inlet velocity was constant, the
pressure was assumed to be a function of the streamwise coordinate only, and the

second derivative with respect to the streamwise coordinate was neglected.

The laminar ﬂow development in the entry region of ducts of any cross section,
including circular tubes and rectangular ducts, were also analyzed by Sparrow et al.
(1964). Their analytical solution, which is modiﬁed from the method due to Langhaar
(1942), involved a stretched axial coordinate and a linearization of the inertia terms
of the equation of motion; hence, the boundary layer model does not need to be
invoked. A correction term was used to compensate for the accumulated wall shear

stresses exerted by the developing ﬂow relative to the fully-developed ﬂow.

Experimental velocity measurements in the entry region of tubes with well-
rounded inlets were ﬁrst made by J. Nikuradse and published by Prandtl and Tietjens
(1934). Results showed that the velocity distribution around the core close to the
inlet section is independent of friction and that the parabolic reduction in velocity
toward the walls is justiﬁed. After the boundary layers meet at the axis, the ﬂat
velocity at the core ceases to exist. Nikuradse’s tests checked well with the
theoretical values obtained by Boussinesq (1891) in the region away from the inlet
section and were in good agreement with the calculations of Schiller (1922) only in
the ﬁrst third of the entrance length. The complete transition to fully-developed ﬂow

according to Schiller (1922) seemed to be rather premature. This discrepancy is

13
apparently due to the acceleration of the ﬂuid near the axis which causes the pressure
drop in the inviscid core to increase compared to that of a fully-developed ﬂow. At
larger distances from the inlet section, Schiller’s calculations predicted a high rate of

increase of centerline velocity. As a result, his entrance length was rather short.

The numerical solution for laminar entry ﬂow in a pipe by Hornbeck (1964) used
essentially the same procedure as Bodoia and Osterle (1961), which is discussed in
Section 2.2. Christiansen and Lemmon (1965) numerically analyzed the uniform entry
ﬂow in a pipe on the basis of linearized inertia terms without postulating a boundary-
layer model. They also compared their results with previous theoretical and
experimental results. The work of Vrentas et al. (1966) used the boundary-layer
equations and the method of Wang and Longwell (1964), as outlined in Section 2.2.

Laminar ﬂow in the entry region of a pipe at low and moderate Reynolds
numbers (in the range 0 to 500) was also studied numerically by Friedmann et al.
(1968) without the boundary-layer approximations. They observed the overshoot
phenomenon. To determine whether the bulges in the velocity proﬁles might not
merely be a consequence of singularity at the leading edge, they treated the limiting
case of vanishingly small Reynolds number analytically and modiﬁed the axial velocity
proﬁle near the wall to eliminate the singularity. Their results showed that the bulges
persisted in the region near the wall. The entry ﬂow into blood vessels was modeled
by Lew and Fung (1970) as steady axisymmetric ﬂow in circular tubes. A numerical
solution, without the boundary layer approximation, was provided using a series of

two sets of eigenfunctions for Reynolds numbers between 0 and 100.

14

2.2 Entry Flow in a Channel

The laminar, incompressible ﬂow in the entry region of a straight channel is an
example of two-dimensional boundary-layer ﬂow. Owing to viscous friction,
boundary layers are formed on both walls, and their thickness increases as the ﬂow
propagates downstream. The development of the ﬂow is primarily concerned with
changes in the streamwise velocity component over the entrance length. Due to the
inﬂuence of viscosity, the velocity proﬁle gradually changes to assume the Poiseuille

parabolic distribution at a sufﬁciently long downstream distance.

Since the pioneering work of Prandtl, boundary-layer theory has been the
centerpiece to all analytical solutions of the entry ﬂow. Although considered as a
powerful tool, it is well known that the boundary-layer assumptions are invalid in the
vicinity of the leading edge of a ﬂat plate as is the case near the inlet section to a
channel. In this region the derivative azu/ax2 is not negligible relative to aZu/ayZ.
In addition, the y-component velocity is not negligible and the pressure gradient in
the y-direction is not necessarily small so that the y-component momentum equation

does not vanish.

In his early work, Schlichting (1934, 1979) analyzed the ﬂow in the entry region
of a two-dimensional channel using the boundary-layer model. The entry region was
considered to be composed of two zones: ﬁrst, a section of two boundary layers
developing on the walls with an inviscid core near the centerline followed by a section
where the velocity proﬁle gradually develops into the parabolic distribution of
Poiseuille ﬂow. At small distances from the inlet section, he assumed that the two

boundary layers on the lower and upper walls will grow in the same way as they

15
would along a ﬂat plate at zero incidence. The resulting velocity consists of two
boundary-layer proﬁles near the walls joined by a ﬂat proﬁle at the core. Since the
ﬂow rate at each and every cross section must be the same, the reduction in the ﬂow
rate near the walls due to friction must be compensated for by a corresponding
increase near the centerline. Mathematically, Schlichting used the continuity equation
and the deﬁnition of the displacement thickness to deﬁne the velocity outside the
boundary layer as a function of the streamwise coordinate. Following a procedure
similar to Blasius’ solution for laminar boundary-layer flow over a ﬂat plate, he
expanded the stream function in a series and integrated in the downstream direction
to ﬁnd the boundary-layer growth and the velocity distribution. The boundary-layer
growth in this case is distinguished from Blasius’ solution by the effect of an
accelerated external ﬂow. The boundary-layer solution is then matched with the

uniform ﬂow at the core.

In the region after the two boundary layers meet, where the ﬂow is nearly fully-
developed, Schlichting perturbed the Poiseuille velocity proﬁle by performing an
asymptotic power series expansion leading to progressive deviation from the
parabolic distribution as the integration proceeds in the upstream direction. The two
solutions, expressed in series form with a sufﬁcient number of terms, are patched at
some appropriate streamwise location where they are simultaneously valid. Thus,
the ﬂow in the entire entry region is determined. The entrance length was calculated

by Schlichting as L,/H = 0.04 Re, where H is the height of the channel and Re is the

Reynolds number.

Roidt and Cess (1962) improved Schlichting’s (1934) solution by the inclusion
of additional terms in the perturbation series. Van Dyke (1970) introduced an

16

asymptotic solution that is uniformly valid for large Reynolds numbers. He corrected
Schlichting’s technique by adding second-order terms for the stream function to
account for the displacement thickness in the boundary-layer region near the inlet;
due to the displacement effect of the boundary layer, the inviscid core near the inlet
does not resemble a uniform ﬂow. The ﬁrst approximation of Van Dyke’s expansion
was the leading edge solution for a semi-inﬁnite plate, which was calculated earlier
by Davis (1967). Van Dyke also pointed out that the vorticity at the inlet is not zero
because of its upstream diffusion. A similar procedure was used by Wilson (1971);
Wilson (1969) considered the transition in the proﬁle-development region to a fully-
developed ﬂow in more detail.

In their works, Collins and Schowalter (1962, 1963) showed that the boundary-
layer' solution approaches the ﬁnite-difference solution (Bodoia and Osterle, 1961) by
considering more terms than had previously been retained by Schlichting (1934).
Abarbanel et al. (1970) used a similarity solution to check whether the overshoot in
the velocity proﬁle of channel ﬂow is part of the exact solution and not a numerical
consequence of the truncation error. They investigated the Stokes’ ﬂow over a half-
plane with ﬂat inlet velocity proﬁle. Their velocity proﬁles exhibited the same sort of
bulges as had been found in the numerical solutions. Moreover, the ﬂow pattern in

the vicinity of the leading edge agreed with that of channel ﬂow.

Recently, Li and Ludford (1980) examined the entry ﬂow analytically to
determine whether the overshoot in the velocity proﬁle near the channel walls, found
numerically, is real; both uniform and special irrotational inlet conditions were
considered. They used a second-order boundary-layer analysis, up to a Reynolds
number of 1000, and concluded that the axial velocity must have an overshoot at all

17
Reynolds numbers. This overshoot was found to lie inside the boundary layer for the

uniform inlet condition and at the edge for the special irrotational inlet condition.

Wiginton and Wendt (1969) simpliﬁed the transformation used in the
linearization solution of Sparrow et al. (1964) and extended their techniques to solve
the eigenvalue problem for two-dimensional ﬂows. Also, Wiginton and Dalton
(1970) modiﬁed the technique used by Sparrow et al. (1964) to analyze the laminar
ﬂow in the entry region of a rectangular duct with different aspect ratios; they noted
that the pressure gradient reached a fully-developed state in about half the distance
of the velocity proﬁle. An analytical method for determining the pressure drop due to
ﬂow development in the entry region of ducts of arbitrary cross sections, without
solving for the velocity ﬁeld, were formulated by Lundgren et al. (1964). They used
the momentum and energy equations with correction factors to account for the change

in momentum between the inlet and downstream sections.

An experimental investigation of the development of the velocity and pressure
ﬁelds for a laminar ﬂow in the entry region of a rectangular duct was carried out by
Sparrow et al. (1967). Reynolds numbers were in the range 1000 to 5000. An
interesting phenomenon was the higher velocity seen in the lower half of the duct as
compared to the upper half. Comparison of their data with the results of other authors
yielded a close resemblance. Feliss et al. (1977) also examined the incompressible
channel flow experimentally and found that the entrance length, as a linear function of
the Reynolds number, correlated well with the analytical predictions of Schlichting
(1934).

l8

Numerically, Bodoia (1959) and Bodoia and Osterle (1961) were the ﬁrst to
obtain a solution for the development of plane Poiseuille and Couette ﬂows of an
incompressible ﬂuid in a straight channel. They used a ﬁrst-order ﬁnite-difference
approximation to Prandtl’s boundary-layer equation. Their solution was for Reynolds
number of 75. Wang and Longwell (1964) studied the steady-state laminar ﬂow in
the entry region of a channel at 11 Reynolds number of 150 with both uniform and
cascade inlet conditions. They employed Southwell’s relaxation method (Allen and
Southwell, 1955) to solve the complete Navier-Stokes equations expressed in terms
of the vorticity and the stream function, and used a ﬁrst-order backward-difference
approximation to stabilize the numerical solution. Their vorticity boundary condition,
borrowed from Thom (1933), was ﬁrst-order as well. The results were in better
agreement with Schlichting (1934) in the region downstream of the inlet. The
overshoots in the velocity proﬁle near the inlet section were nonexistent for the

cascade inlet condition.

A solution for Reynolds number of 250 was obtained by Gillis and Brandt (1964)
and echoed in their subsequent work (Brandt and Gillis, 1966) in which they extended
the method to the equations of magrretohydrodynamics, speciﬁcally to the problem of
the development of Hartrnann ﬂow. A fourth-order equation in the stream function
was used for the linearized Navier—Stokes equation. As in the case of Wang and
Longwell (1964), the uniform inlet ﬂow develops proﬁles with local minima on the
channel axis and a pair of symmetrically placed maxima on either side of it. These
kinked proﬁles persist for a distance downstream of the order of magnitude of the

channel height.

l9

McDonald et al. (1972) solved the Navier-Stokes equations in the entry region
of both a tube and a channel for Reynolds numbers between 0.75 and 37 500. They
developed a stable numerical scheme based on the Alternating-Direction Implicit
method m solve the vorticity and stream function equations; their representation of
the vorticity boundary condition was second-order accurate. An overshoot in the
velocity proﬁle existed for the uniform inlet case but was less pronounced for a zero-
vorticity inlet condition. Using quasi-linearized Navier-Stokes equations in terms of
the primitive variables, Morihara and Cheng (1973) solved the entry ﬂow between
parallel plates; their quasi-linearization involved a truncated Taylor series expansion
of the ﬁrst-order. A numerical solution for Reynolds numbers in the range 0 to 1000
was provided. The effect of the upstream condition on the laminar ﬂow development
in a channel for Reynolds numbers between 1 and 1000 was studied numerically by
Sparrow and Anderson (1977). Their results showed that, due to the upstream
interaction at low and intermediate Reynolds numbers, the sharply turning ﬂuid near
the walls tends to overshoot at ﬁrst then readjust as it progresses downstream.

In all the numerical solutions mentioned above, either some simplifying
assumptions were made, such as the assumptions inherent in boundary-layer theory
(that is, both the streamwise derivative azu/ax2 and the pressure gradient normal to
the walls, ap/By, were neglected) and in the linearization of the Navier-Stokes
equations, or a first-order vorticity boundary condition was used, or the numerical
solution of the full Navier-Stokes equations was presented for relatively low

Reynolds numbers.

The recent numerical work by Abdulla (1987) used a successive-over relaxation

procedure including accelerating parameters to solve the N avier—Stokes equations for

20
the may ﬂow for Reynolds numbers up to 2000. Although he used a second-order
vorticity boundary condition, his difference equations were only ﬁrst-order accurate;
he used forward] backward diﬂ‘erencing to improve the convergence of the Navier-

Stokes equations expressed in terms of the vorticity and the stream function.

A study of transient magnetohydrodynamic duct ﬂow was conducted by
Lundgren et al. (1961). They considered the parallel ﬂow of an electrically conducting
viscous, incompressible ﬂuid with transverse magnetic ﬁeld and determined the exact
solutions for the velocity and magnetic ﬁelds in the form of convolution integrals.
Using a ﬁnite-difference solution, Hwang and Fan (1963) also investigated the
development of velocity proﬁles in the entry region of a ﬂat rectangular duct for an

electrically conducting ﬂuid and transversely applied magnetic ﬁeld

The subject matter of the next few sections is concerned with further literature
review of work relevant to this study. A wealth of analytical, numerical, and
experimental studies, for steady and time-dependent ﬂows, are cited in the
literature. A comprehensive survey is bound to be cumbersome and, therefore, will
not be attempted. However, a number of investigations (predominantly numerical),
with pertinence to the ﬂow problem and numerical procedure of this study, are

reviewed. The studies discussed hereafter are conveniently classiﬁed by ﬂow type.

2.3 Internally Separated Flows

Internal ﬂows mainly containing separation and wakes are of interest to this

study. These include ﬂows with sudden contractions or expansions, ﬂows in

21
constricted channels and pipes, ﬂows in conduits with obstructions, ﬂows in driven

cavities, and ﬂows past blunt or streamlined bodies placed in conﬁned ﬂow ﬁelds.

Greenspan (1969a) studied the steady ﬂow in a channel with both a forward-
facing and a ﬁnite step for Reynolds numbers up to 10 000; the approaching ﬂow was
considered to be parabolic. He used the Successive Over-Relaxation method with
accelerating parameters and ﬁrst-order boundary vorticity to solve the vorticity and
stream function equations. An approach consisting of a numerical procedure similar to
that of Greenspan (1969a) for low Reynolds numbers and modiﬁed with upwind
differencing for moderate and high Reynolds numbers, was used by Friedman (1972)
to study the laminar ﬂow in a channel with a ﬁnite step. Optimum accelerating

parameters were also used to reduce the computer time needed for convergence.

Taylor and Ndefo (1970) used the method of splitting, which was originally
developed by Yanenko (1971), to compute the ﬂow in a two—dimensional channel
with straight walls and a backward-facing step for Reynolds numbers up to 100.
Abdulla (1987) applied the same over-relaxation procedure that he used for entry
flow (which is also similar to that of Greenspan (1969a)) to predict the laminar ﬂow
separation in a channel constricted with a forward-facing step, a backward-facing

step, and a ﬁnite step.

An experimental and numerical investigation of a two-dimensional channel ﬂow
over a single backward-facing step was performed by Armaly et al. (1983). They
reported the downstream velocity distribution and reattachment length for Reynolds
numbers between 70 and 8000. Laser-Doppler measurements revealed the presence

of separation regions downstream of the step and on both sides of the channel, in

22
addition to the primary recirculating zone. The numerical results obtained showed
good agreement with the experiment. Results for the time-dependent impulsively
starting backward-facing step mounted in a two-dimensional duct were presented by
Durst and Pereira (1988). They used a ﬁnite volume numerical approach for Reynolds
numbers 10, 389, and 648. Experimentally, they showed that no recirculating ﬂow

region existed on the wall opposite to the step for Reynolds number S 420.

A numerical solution of the N avier—Stokes equations for steady, laminar ﬂow in
a channel constricted with a semi-inﬁnite forward-facing symmetric step on both
walls, was provided by Dennis and Smith (1980). They studied separation for
Reynolds number up to 2000 using different size mesh. Comparisons of the upstream
separation position and the wall vorticity showed good agreement with the
asymptotic solution of Smith (1979). The effect of an asymmetric downstream
disttn'bance, in the form of a constriction, branching, bend, or ﬂow injection, in a
channel was also studied by Smith (1977). The boundary-layer equations were
solved numerically. The ﬂow becomes separated on the wall ahead of the disturbance
and the ﬁnite disturbance prohibited a smooth joining with the original Poiseuille ﬂow.

Two different types of internally separated ﬂows were studied by Kumar and
Yajnik (1980). The ﬂows in a channel with a symmetric sudden expansion and in a
channel "with a base" at large Reynolds numbers were solved by means of a small
perturbation in the eigenfunctions of the Poiseuille ﬂow; this reduces the problem to
solving nonlinear ﬁrst-order ordinary differential equations. Results for the problem
with the sudden expansion compared well with the ﬁnite-difference solution obtained
from an iterative procedure. The separated ﬂow due to an obstruction placed in a two-

dimensional channel was explored by Nallasamy (1986). The characteristics of the

23
separated ﬂow behind the obstruction were obtained from the ﬁnite-difference

solution of the N avier-Stokes equations for Reynolds numbers up to 1500.

A solution of the N avier-Stokes equations of a conﬁned wake problem was
obtained numerically by Paris and Whitaker (1965). The two—dimensional problem
was formed by the merging of two Poiseuille ﬂow streams. The Altemating-Direction
Implicit method was used in the solution of the vorticity and stream function
equations for Reynolds numbers as high as 647. At high Reynolds numbers, the
boundary layers equations, without an order of magnitude analysis (i.e., the y-
momentum equation was not neglected), were used. The results were in agreement

with the solution of the complete equations.

The ﬂow in a circular pipe with recessed walls was numerically investigated by
Friedman (1970). Separation and vortices were studied for Reynolds numbers up to
500. Also, the steady, laminar, axisymmetric ﬂow through vascular stenoses was
modeled in the study of Deshpande et al. (1976) as ﬂow in a rigid tube with a
contoured constriction. They used the full N avier-Stokes equations in cylindrical
coordinates. For the Reynolds numbers at which the ﬂow is known to be laminar,
their numerical results converged without difﬁculty and compared well with

experimental results.

2.4 Flow in a Driven Cavity

Because of the nonlinearity of the Navier-Stokes equations governing the ﬂow
in a driven cavity, and because of the presence of regions in which the solution

changes rapidly, a large number of authors used the driven cavity problem to analyze

24
their numerical techniques. The importance of the problem to ﬂuid mechanists,

however, stems from its relevance to separated ﬂows.

One of the earliest attempts at solving this problem was made by Kawaguti
(1961) who studied the ﬂow in a rectangular cavity with different length/depth ratios.
He integrated the Navier-Stokes equations numerically and found out that the
relaxation procedure failed to converge at high Reynolds numbers. The Navier-
Stokes equations were solved iteratively by Mills (1965) for the ﬂow generated in a
driven rectangular cavity at a Reynolds number of 100. Mills used the numerical
procedure advocated by Thom (Thom and Apelt, 1961).

Burggraf (1966) assessed the effects of nonlinearities and of noncircular
boundaries by considering numerical solutions of the full N avier-Stokes equations for
a square cavity. First, his analysis, based on a linearized theory restricted to an eddy
having a circular boundary, revealed the development of a recirculating eddy. He then
conﬁrmed the validity of the linear analysis as a description of separated eddies by
the results of the numerical solution. These solutions were carried out by means of a
relaxation procedure for Reynolds numbers in the range 0 to 400; the limiting cases,
as the Reynolds number goes to zero or infinity, were solved analytically. At higher
Reynolds numbers an inviscid core develops, but secondary eddies are present in the
bottom corners of the square at all Reynolds numbers. Results obtained numerically
as well as experimentally by Pan and Acrivos (1967) for a rectangular cavity,
confumed the presence of a single inviscid core of uniform vorticity with viscous
effects being conﬁned to thin shear layers near the boundaries. Using a successive

over-relaxation method with accelerating parameters, Greenspan (1967, 1968,

25
1969b) also studied the ﬂow in a driven cavity. His numerical solutions were for

Reynolds numbers between 50 and 100 000.

Numerical solutions of the elliptic transport equations of a laminar ﬂow with
heat transfer for a two—dimensional driven cavity problem were presented by Runchal
et al. (1969). Their solution featured unsymmetrical formulation of convection terms,
Gauss-Seidel iteration procedure, and a grid with a nonuniform mesh. The inﬂuence
of mesh size and the vorticity boundary condition on the accuracy and convergence of
the solution were considered for a wide range of Reynolds numbers. Their
computation results agreed with previous knowledge and showed that secondary
eddies appear in the lower corners of the cavity, even for creeping ﬂows. Bozeman
and Dalton (1973) obtained solutions for the ﬂow in rectangular cavities by solving
various implicit approximations of the Navier-Stokes equations. They determined the
effect of the Reynolds number and mesh size on the different schemes and also
examined the effect of vorticity boundary conditions; the behavior of a vortex center

and secondary vorticity were exhibited as well.

Several numerical techniques for the treatment of the viscous, incompressible
ﬂow in a driven cavity were reported by NASA (Rubin and Harris, 1975). Their
results showed that convergence is more rapid and accurate with the vorticity-stream
function formulation compared to the solution of the equations of motion in terms of
the primitive variables. The authors also attributed the higher-Reynolds-number
solutions to the stability of the Alternating-Direction Implicit method. In particular,
Morris (1975) used the Alternating-Direction Implicit method with a ﬁrst—order

vorticity boundary condition to solve the vorticity equation.

26

The two—dimensional driven cavity ﬂow remains as a model problem in the
study of different numerical schemes. Roache (1975) used the driven cavity problem
to compare three different numerical methods and also to examine the effect of
boundary conditions on his solution. This was also the case with Rubin and Khosla
(1977) in their development of a higher-order approximation for viscous ﬂow
calculations. A review of many numerical methods, complete with error and
convergence comparisons, for computing recirculating ﬂows in the cavity problem, was
presented by Tuann and Olson (1978). Benjamin and Denny (1979) discussed
convergence properties of several ﬁnite-difference schemes for solving the two-
dimensional ﬂows in a cavity at Reynolds numbers as high as 10000. They used a
relaxed Alternating-Direction Implicit procedure in their second-order accurate
solution of the vorticity and stream function equations. Ghia et al. (1982) used the
cavity ﬂow problem to test the effectiveness of their numerical method; they used a
ﬁne mesh in hi gh-Reynolds-number ﬂows. In addition, the method of false transient
was used by Mei and Plotkin (1986) in a ﬁnite-difference scheme applied to the

driven ﬂow in a square cavity for Reynolds numbers up to 5000.

Lastly, Duck (1982) studied the ﬂow inside a square cavity generated by an
oscillating motion of the upper boundary. He solved the unsteady Navier—Stokes
equations for Reynolds numbers as large as 600. His results showed that the steady
component of the ﬂow comprises two pairs of equal counter-rotating eddies situated
at the upper and lower halves of the cavity.

27

2.5 Externally Separated Flows

External ﬂows such as the ﬂow over a ﬂat plate and past cylindrical, spherical,
and other different shape objects are associated with separation. Two of these

studies for ﬂows with boundary-layer separation are of particular interest.

heal and Acrivos (1969) obtained a numerical solution for the Navier-Stokes
equations for the steady ﬂow over a ﬂat plate in which a constant ﬁnite slip velocity in
the direction opposite to the stream is imposed. Their study was concerned with the
structure of the closed streamline ﬂows within a boundary layer. These streamlines
were found to be embedded in the viscous layer adjacent to the plate at all Reynolds
numbers; Reynolds numbers were up to 80. The parabolic boundary-layer equations
were found suitable for the solution of the ﬂow including the recirculating zone near

the plate.

Briley (1970, 1971) analyzed a two—dimensional, viscous, incompressible ﬂow
containing separation bubbles by patching the numerical solution in the immediate
vicinity of the bubble to the boundary-layer solution and the inviscid solution of the
outer ﬂow. He presented numerical solutions for ﬂows with Reynolds numbers of the
order of 50 (based on momentum thickness). Briley used a fourth-order Lagrangian
polynomial for interpolating the stream function to derive a formula for the vorticity at
the boundary based on a special relationship for the x-component velocity at the
points adjacent to the boundary, thus providing higher order coupling between the
vorticity and the stream function at the wall. The Altemating—Direction Implicit

method, with an iterative procedure and variable time steps, was used to solve both

28
the vorticity and stream function equations; he preferred the marching approach in
solving his steady-ﬂow problem.

2.6 Time-Dependent Flows

A general treatment of the time-dependent, two-dimensional ﬂow of viscous,
incompressible ﬂuids can be found in many publications. Several time-dependent flow

solutions are outlined below.

A numerical method for time-dependent problems concerning the ﬂow of
viscous, incompressible ﬂuids was applied by From and Harlow (1963) to study the
development of a vortex street behind an impulsively accelerated plate placed in a
channel. The Reynolds number range was between 15 and 6000. The numerical
scheme employed was discussed in more detail by From (1963, 1964). It is an
explicit scheme with an iterative procedure based on Liebmann’s method, as outlined
by Thom and Apelt (1961); a special fust-order vorticity boundary condition was
derived. Fromm also analyzed the stability of the vorticity equation. In his later
work, From (1972) described his method for the numerical computation of the
Navier-Stokes equations at high Reynolds numbers and its application to several
time-dependent, viscous incompressible ﬂow problems. He also presented

alternative discretization schemes for the convective derivatives.

A general ﬁnite—difference computational procedure for time-dependent viscous,
incompressible ﬂow problems was described by Pearson (1964, 1965a, 1965b). It
involves solving a nonlinear, fourth-order equation in three independent variables and
the use of a second-order approximation for the boundary vorticity. The method was

29
applied to examples such as the injection of ﬂuid into a rectangular region, the time-
dependent rotating disk problem, and vortex formation and break up in a two-
dimensional region. He also discussed the choice of difference equations, relaxation
procedures, approximation to the boundary conditions, and the resulting

computational stability, speed, and accuracy.

Pao and Daugherty (1969) studied the time-dependent, two-dimensional,
viscous, incompressible ﬂow past a ﬁnite ﬂat plate. For the vorticity equation, they
used both an explicit method and the Alternating-Direction Implicit method with
forward] backward differencing, while for the stream function equation, they used the
iterative Successive Over-Relaxation method. The vorticity boundary condition was
approximated by the second-order expression due to Jenson (1959). Also, three
types of coordinates were used to construct the computational grid. The details of the
ﬂow ﬁeld around the plate were computed for Reynolds numbers equal to 4 and 993.
It was found that the velocity at the outer edge of the boundary layer exceeds that of
the free stream and that the Alternating-Direction Implicit method was two to three
times faster than the explicit method. The accuracy of the different numerical methods

and the effects of the far ﬁeld boundary conditions and mesh size were also studied.

A numerical procedure for a two-dimensional, compressible, unsteady boundary-
layer ﬂow was presented by Kwon et al. (1988). They studied the sudden
acceleration of an incompressible ﬂuid over a semi-inﬁnite ﬂat plate and the steady
sinusoidal oscillation of a plate in a semi-inﬁnite ﬂuid, better known as Stokes’
second problem. Their results were in good agreement with the available theoretical

solutions.

30

2.7 Numerical Treatments of Flow Problems

A numerous number of numerical techniques have been applied to ﬂow problems
of various sorts. As a closure to this review, a number of these numerical treatments
will be brieﬂy highlighted; they are mostly studies related to the development of the

present numerical scheme.

In their integration of the Navier-Stokes equations for steady, two-dimensional
ﬂow, Dennis and Chang (1969) studied the convergence of the iterative procedure
used to solve the steady ﬂow past a ﬁnite plate and elliptic cylinders. The vorticity
and stream function were used as dependent variables. They discussed the ﬁnite-
difference approximation of the vorticity equation necessary to secure the diagonal
dominance of the system while still retaining the higher-order accuracy of the central
differences. Gosman et al. (1969) also analyzed the use of upwind differences to

provide convergence to numerical solutions of ﬂow problems.

Several numerical methods for inviscid and viscous ﬂows were compared by
Taylor et al. ( 1972). For viscous ﬂows, they concluded that implicit schemes are
preferable to explicit or Dufort-Frankel differencing and that for the inertial terms, any
scheme which is stable for inviscid ﬂow and has controllable artiﬁcial viscosity at low
Reynolds numbers, may be used. A more thorough review of methods used in the
numerical simulation of viscous incompressible ﬂows was given by Orszag and Israeli
(1974). They discussed the choice and accuracy of available numerical
approximations and cited an enormous number of references. They also studied the
error associated by several formulas (i.e., those of Briley (1971), From (1963),
Jenson (1959), Woods (1954), and Thom (1933)), used in evaluating the vorticity on

31
the no—slip boundaries and showed that results given by the different formulas can
differ.

Another extensive comparison between results given by the above various
formulas in the case of a driven cavity ﬂow has been made available by Gupta and
Manohar (1979). To identify possible causes of instability, the effect of the ﬁnite-
difference approximations to the no-slip conditions on the linear stability of the
numerical solution of unsteady ﬂows was analyzed by Taylor (1970). All conditions
used were shown to be stable without the presence of suction and conditionally
stable if there is suction; in this case, the choice of mesh size is restricted. The

selection of the proper integration method was shown to be critical.

Khosla and Rubin (1974) suggested a diagonally dominant second-order
accurate implicit scheme to replace the ﬁrst-order upwind differencing usually used to

stabilize the solution of the nonlinear Navier-Stokes equations.

Finally, the use of optimized over-relaxation factors was discussed by
Thompson (1969). The number of iterations required for convergence for the time-
dependent Navier-Stokes equations, with iterations performed at each time step,
was reduced by more than one order of magnitude over Gauss-Seidel iteration. The

restrictions placed on the maximum size of the time step were also removed.

DESCRIPTION OF THE PROBLEMS

The sequence of the basic ﬂow problems and their mathematical modeling are
discussed in this chapter. The describing equations of motion, the Navier-Stokes and
the continuity equations, can be expressed in different forms. Here, the vorticity-
stream function formulation is adapted for the solution. The initial and boundary

conditions required in the solution are also presented.

3.1 The Flow Problems

A number of ﬂow problems are investigated in this study. The ﬂow past a
moving boundary section placed within the entry region of a rectangular conduit is of
primary interest. The physical situation is illustrated by the two-dimensional,
isothermal, incompressible, laminar ﬂow of a Newtonian, isotropic, homogeneous ﬂuid
with constant density and kinematic viscosity in the entry region of a straight
channel. The channel is made up of two semi-inﬁnite, ﬂat plates lying parallel to one
another in the xz-plane; the x-axis of the rectangular Cartesian coordinate system
coincides with the lower wall and the y-axis is normal to the ﬂow. At the inlet where

x = 0, a uniform ﬂow with velocity U is assumed across the channel of height H. A

32

33

y t Uniform Profile Parabolic Proﬁle

 

 

 

Figure l. A schematic representation of the ﬂow problem including the moving
boundary.

parallel Poiseuille ﬂow is considered to exist far downstream. The deﬁning geometry
is shown in Fig. 1. The moving section of the boundary, which is assumed to be
moving either steadily or oscillating periodically with a constant frequency and
speciﬁed amplitude, has a velocity UB and is deﬁned to exist between the streamwise

locations x = XI and x = X2.

The steady ﬂow in the entry region of a plane two-dimensional channel is
considered ﬁrst. This entry ﬂow gives rise to a growth of boundary layers on the two
walls. The most detailed description to this problem has been that of Schlichting
(1934, 1979). The ﬂuid enters the channel with a uniform velocity across the inlet
section and its velocity is assumed to be zero at the wall for all positive x. Due to the
inﬂuence of friction, the ﬂuid lying near the walls is retarded; the viscosity makes
itself felt in a layer near the wall. The boundary layer grows in thickness as the ﬂow

34
progresses downstream. The inner part of the ﬂow is inviscid and its velocity is
essentially uniform. Since the continuity of mass requires the ﬂow rate at all cross-

sections to be constant, the ﬂuid near the centerline of the conduit is accelerated.

At some streamwise distance from the inlet section, the two boundary layers
merge together at the centerline (when the boundary-layer ﬁlls half the channel’s
height). At this stage, the pressure drop per unit length is higher than that for the
fully-developed ﬂow. The ﬂow will continue to accelerate until equilibrium between
viscous forces and pressure forces is ﬁnally reached at a location downstream. As a
result, the initially uniform velocity proﬁle is transformed asymptotically into the
parabolic velocity distribution of a Poiseuille ﬂow. Since Bulax at 0 in the entry
region, the ﬂow is not one-dimensional, and the velocity components depend on the

streamwise coordinate x as well as the normal coordinate y. The entrance length,

, t Invlscid Core Parabolic Profile

    

 

 

4 L. :L: r.d __.

l

Inv'ncid-Core Length Profile-Development Length
L _

Entrance Length

 

 

 

 

Figure 2. Laminar entry ﬂow in a two-dimensional channel.

35
which includes the inviscid-core length and the proﬁle-development length, is
considered to exist whenever the velocity at the centerline is less than 99% of the
parabolic centerline velocity (Prandtl and Tietjens, 1934). The ﬂow and its distinct

regions are shown schematically in Fig. 2.

The channel is then constricted by placing a ﬁnite step on its lower boundary.
As a result, the ﬂow becomes separated before and aft of the step. A moving section
introduced on the lower boundary is also studied. When the boundary section moves
uniformly in its own plane with a constant velocity in the direction opposite to the
ﬂow, a separated region exists in the near vicinity of the moving section. Motion of
the boundary section in the direction of the ﬂow, however, causes the neighboring
ﬂuid to accelerate without the presence of any separation; the problem of Couette-
Poiseuille ﬂow results when the entire lower boundary moves in the direction of the
flow. The two cases of boundary motion with and against the direction of the ﬂow
also occur in the case of an oscillating boundary section where the ﬂow becomes time-
dependent. A schematic representation of the sequence of ﬂow problems is shown in

Fig. 3.

3.2 Mathematical Formulation

Viscous, incompressible ﬂows are described by the Navier-Stokes and the
conservation of mass equations. To facilitate the solution, the differential equations
of motion are expressed in non—dimensional form. The following characteristic

quantities are selected for normalizing the various variables of this study:

length = H, velocity = U, time = %. pressure = pUz (1)

rig

 

a) Entry ﬂow in a two-dimensional channel.

 

b) Entry ﬂow with a step.

 

 

c) Entry ﬂow with a moving lower boundary section.

 

 

<—>

d) Entry ﬂow with an oscillating lower boundary section.

Figure 3. Schematic representation of the sequence of entry ﬂow problems.

37
where H is the distance separating the lower and upper boundaries, U is the uniform

inlet velocity, u) is the frequency of oscillation of the moving boundary section, and p

is the density of the ﬂuid. With the above characteristic quantities, the following

 

 

dimensionless variables result:
* it: t
x = x , t = 0) t* ’ v = V , p = —L2— (2)
H U pU

where the boldface refer to vector quantities and the asterisk is used to denote a
dimensional variable. In the above deﬁnitions, x is the position vector (xi + yj +
zk ), t is the time, V is the velocity vector (ui + V} + wk), and p is the pressure; x, y,
and z are the three spatial coordinates; u(x,y,t), v(x,y,t), and w(x,y,t) are the
velocity components in the x-, y—, and z-directions respectively; and i, j, and k are
unit vectors in the x—, y-, and z-directions, respectively. Hereafter, all variables and

equations are expressed in dimensionless form unless otherwise stated.

Applying the transformation of variables indicated above, the dimensionless

Navier-Stokes equation in vector form is written as

3V V V V V 1 v2‘
at ( ) pk Re ( )

where V is the gradient operator and V2 is the Laplacian operator, both in non-
dimensional form. The Strouhal number, St, is deﬁned as

St = —- <4)

and the Reynolds number, Re, based on the channel’s height is given by

38

Re == —— (5)

where v is the kinematic viscosity. In Eq. (3), since the body forces do not inﬂuence
the velocity ﬁeld of conﬁned ﬂows, the effect of body forces is ignored by including the

body-force term in the pressure term to form the kinetic pressure, i.e.,
_ 2* (6)
where p is the pressure, g is the gravitational acceleration, and z“ is the vertical
coordinate; gz‘lU 2 is the dimensionless Froude number.

Furthermore, the conservation of mass equation of an incompressible ﬂuid is

v.v = o . (7)

It is evident that the ﬂow is inﬂuenced by two dimensionless parameters,
namely, the Reynolds number and the Strouhal number. The Reynolds number
represents the ratio of inertial to viscous forces, while the Strouhal number reﬂects

the effect of oscillatory boundary conditions.

For a plane ﬂow independent of the z-coordinate, we have w = 0 and em = 0.
Then the three dimensionless momentum equations and the continuity equation, as

given by Eqs. (3) and (7), respectively, in component form are reduced to

.22. _a_u_ super. _1 __a2u _82u
St at+u3x+vdy 3x+Re (8x2+3y2) (8)

 

.31. ﬂ. ELLE}; 1 32v _32v
St dt+udx+vay 3y+Re (3x2+ay2) (9)

39

an 3v
= 1
3x + ay 0 ( 0)

 

where Eqs. (8) and (9) are the dimensionless x- and y-component Navier-Stokes
equations, respectively, and Eq. (10) is the two-dimensional continuity equation.
The above three equations are the describing equations of motion in terms of the

three primitive variables u(x,y,t), v(x,y,t), and p(x,y,t).

3.3 The Problem for Solution: Vorticity-Stream Function Formulation

In general, solving Eqs. (8)- (10) with the appropriate boundary conditions
yields the velocity and pressure ﬁelds directly. The continuity and the Navier—Stokes

equations, however, can be presented in several alternative forms for the solution.

For two-dimensional ﬂows, it is often desirable to introduce the vorticity and
the stream function as dependent variables. This approach is attractive since it
reduces the number of equations to be solved from three to two; and, second it
ensures that the continuity equation is automatically satisﬁed. Failure to satisfy
continuity using the primitive variables can lead to numerical instability (Mallinson
and de Vahl Davis, 1973). On the other hand, the advantage gained by eliminating
the pressure and its boundary conditions from the governing equations shows up as a

difﬁculty in expressing the boundary conditions on the vorticity.

By differentiating Eq. (8) with respect to y and Eq. (9) with respect to x and
subtracting and invoking Eq. (10), we thus eliminate the pressure terms from the
describing equations and reduce the number of momentum equations from two to one.

There results

4o

3 an an 3 ii. ﬂ. .1. .32.- it.
8' dt(ay dx)+u3x(3y 8x)+v8y(8y ax)
= 1 82 Bu __§v_ + 32 On av (11)
Re 8x2 3y 8x 3y2 3y 8x °

In plane ﬂows, only the z-component C(x,y,t) of the vorticity vector exists; it is

 

 

deﬁned as

3v Bu
= — - —. 12
C ax 3y ( )

Substituting the above deﬁnition into Eq. (11) yield the z-component vorticity

 

equation
a; a; at. 1 32: __§_32 ,
8‘ at H ax H ay Re 8x2 + ay2 (13)

The above equation is referred to as the vorticity transport equation. It states that
the local and convective rates of change of vorticity are balanced by the rate of

dissipation of vorticity.

To fm'ther simplify the problem statement, the x- and y-components of velocity

are expressed in terms of the stream function w(x,y,t) through

u=—aN’—. v=-—a-‘JL. (14)
By 8x

The above deﬁnition satisﬁes the continuity equation identically.

41
It follows that the stream function is related to the vorticity by a Poisson-type
equation derived from substituting the deﬁnition of Eqs. (14) into Eq. (12). The result

is the stream function equation

32 82

Upon substituting for u and v in terms of \y from Eq. (14) into Eq. (13), there

 

 

results
a; 81 BQ _ 31$: 1 821; 32g
St at + By 8x 8x 8y Re (8x2 + By2 . (16)

The above equation is the z-component vorticity transport equation in terms of the
vorticity and the stream function. Equations (15) and (16) are the describing

equations that represent this two-dimensional ﬂow. The result is a nonlinear system

of two partial differential equations which involve the two unknowns C and w.

Because of the nonlinearity of the problem, the solution is by no means
straightforward. Hence, a simultaneous solution of Eqs. (15) and (16) with a
prescribed set of initial and boundary conditions in the region bounded by the
upstream condition, the far ﬁeld downstream condition, and the lower and upper
boundaries is required to give the solution in terms of C and v. Having found v, the
velocities in the x- and y-direction can then be determined using the deﬁnition of the
stream function given by Eq. (14). The above system of equations describing the
motion does not lend itself to an analytical solution, hence, a numerical solution is

attempted to yield the desired information.

42
For the sake of completeness, it is essential to mention that another way of
expressing the describing equations is through a single fourth-order equation in terms
of the stream function alone; this equation can be derived by further substitution of
Eq. (15) into Eq. (16). The advantage of solving two second-order simultaneous
differential equations in the vorticity and stream function instead of one fourth-order
equation lies in the much greater rate of convergence of Poisson’s equation compared

to the fourth-order equation (Woods, 1954).

3.4 The Steady Flow Problem

The steady ﬂow is, of course, a special case of the time-dependent problem.
When the ﬂow is steady, the problem has no independent characteristic time. In this
case, the characteristic time is different ﬁom that given by Eq. (1) and is selected as
H/U. Then, the dimensionless time becomes

 

*
t= if? . (17)

The above deﬁnition allows the time-dependent equations to be reduced to
steady-state equations. The two momentum equations given by Eqs. (8) and (9) are

then rewritten as

 

 

.22. an. _3“_=-_§L 1 32.. a2.
at +0 8x +v 3y 3x + Re (3x2 + Byz) (18)

 

 

_a_v_ .22. ALLA 1 82v 82v
at +u 3x +v 3y 3y + Re (3x2 + Byz) (19)

43

where the Strouhal number has been eliminated. The continuity relation of Eq. (10)
remains unaltered. It is noteworthy that even for the steady ﬂow, the time
derivatives are retained in the equations since the marching in time numerical
methods require the presence of the time derivatives. In such a case, the time is
ﬁctitious. The asymptotic solution of the unsteady equations for large time is
considered to be the solution to the steady-ﬂow problem. This is referred to as a
pseudo-unsteady method (Peyret and Taylor, 1983) or a false transient method
(Mallinson and de Vahl Davis, 197 3).

In terms of the vorticity and the stream function, the describing equation of
motion for the steady ﬂow are represented by Eq. (15) and the following modiﬁcation1
of Eq. (16):

_3£.+§1_§£_-§1_i§_=_1_ ﬁre—2!... (20)
at 3y 3x 8x 8y Re 3x2 8y2

3.5 Initial and Boundary Conditions

To obtain a solution to the above describing equations, we must prescribe the
appropriate initial and boundary conditions. In non-dimensional terms, the inlet is
identiﬁed to exist at x = 0, the lower boundary at y = O, and the upper boundary at y =
1. The initial condition is assumed to be the quiescent state since the ﬂuid is at rest

fort S 0. In terms of velocities,

 

1 When the term 32C/3x2 is not included, Eq. (16) becomes equivalent to Prandtl’s boundary-
layer equation.

44

u(x,y,O) = 0, v(x,y,O) = 0. (21)

The velocity boundary conditions are the uniform velocity upstream, parabolic
proﬁle far downstream, and the no-slip boundary conditions on both the lower and

upper boundaries. For t > 0, these are written as

u(0.y.t) = 1. V(0.y.t) = 0 (22)
u(xZLe,y,t) = 6y - 6y2, v(xZLc,y,t) = o (23)
{ 0 (ﬁxed)
u(x,0,t) = v(x,0,t) = 0 (24)
U3 (moving, XleSXZ)

u(x,1,t) 0, v(x,l,t) = 0 (25)

where U]! is the velocity of the moving section of the lower boundary. The entrance
length L° is deﬁned as the distance from the inlet to the point where the centerline

velocity reaches 99% of the fully-developed parabolic centerline value. At

streamwise locations beyond Le the ﬂow is independent of the x-coordinate. To

ensure that the solution is not inﬂuenced by the downstream boundary condition, we

select the location of the downstream boundary sufﬁciently larger than L.B such that
the solution at the downstream end does not depend on x; the numerical solution is
limited to the region inﬂuenced by x, i.e., x 5 Le. The choice of the x-location of the

downstream boundary condition, given by Eq. (23), is validated by subsequent

examination of the computed solutions.

45

The assumption of a constant inlet velocity was investigated in the experiments
of J. Nikuradse (Prandtl and Tietjens, 1934). Nikuradse’s tests demonstrated that a
uniform velocity proﬁle is justiﬁed for tubes with well-rounded inlets. Furthermore,
the inﬂuence of an actual experimental proﬁle for the inlet velocity, obtained by ﬁtting
a polynomial near the wall to a uniform central section using measurements from a
hot-wire anemometer, on the numerical solution was examined by Abdulla (1987).
The difference in the streamwise velocity distributions along the channel due to the

experimental proﬁles was found to be insigniﬁcant.

The vorticity-stream function system of equations needs two initial conditions
and eight boundary conditions. In terms of the stream function and the vorticity, the

initial conditions are given by
v(X.y.0) = y. C(X.y.0) = 0. (26)

The boundary conditions on C and \y are derived from the known velocity

boundary conditions. From Eqs. (15) and (16), it is apparent that the stream function
equation requires four boundary conditions on \y while the vorticity equation demands
four boundary conditions on C. On the upstream and downstream boundaries, the

stream function for all y-locations is derived from the known velocity by integrating

Eq. (14), i.e.,
w(x,y,t) = ] udy. (27)
Upstream, where the inlet velocity is uniform, we obtain

WOJJ) = Y- (28)

46
On the far ﬁeld boundary downstream, where the streamlines are parallel to the
walls, Eq. (23) and Eq. (27) lead to

V(XZLe,y,t) = 3y2 - 2y3. (29)

Obviously, even with the presence of a moving section, no ﬂuid can ﬂow across
the lower boundary; hence, it forms a streamline which can arbitrarily be assigned a

stream function with zero value, i.e.,
v(X.0.t) = 0. (30)

On the upper boundary the velocity components are zero and the stream function

v is given a value which is equal to the net ﬂow rate; hence, we require that
v(x,l,t) = 1 (31)

where the dimensionless not ﬂow rate is equal to 1.

The vorticity boundary conditions both upstream and downstream are derived
from the vorticity deﬁnition of Eq. (12). Since the y—velocity component is zero on
either of these boundaries, Eq. (12) is then reduced to

___§_u_
C- ay' (32)

It follows that the vorticity on the upstream and downstream boundaries are given by

C(0.Y.t) = 0 (33)

47

C(x2L0,y,t) = 12y - 6.

(34)

On the lower and upper boundaries, the vorticity boundary conditions are not

known explicitly. From Eq. ( 15), however, we can relate the boundary vorticity to the

stream function in order to obtain a meaningful expression of the boundary vorticity.

Since ‘1’ is constant on these boundaries, all its derivatives with respect to x,

including azwaxz, are equal to zero. Equation (15) then gives

82
C(X,0,t) = ‘ -ainl_

-231,

C(X.1.t) ay,

(35)

(36)

The above two relations for boundary vorticity are rather general and need to be

modiﬁed for subsequent use in the numerical solution. In the present form, however,

no modiﬁcation is necessary to accommodate the moving boundary; motion of the

boundary shall be imposed when Eq. (35) is modiﬁed for the numerical solution.

The boundary conditions, in terms of the stream function and the vorticity, used

to solve Eqs. (15) and (16) are summarized below and displayed in Fig. 4.
Upstream: v(0.y.t) = y. C(0.y.t) = 0

downstream (x 2 Le): w(x,y,t) = 3y2 - 2y3, C(x,y,t) = 12y - 6

lower boundary: v(x,0,t) = 0. 9090") = " $322

(37)

(38)

(39)

48

upper boundary: v(x,l,t) = 1, C(XJJ) = - %§-:g- (40)

For the steady-ﬂow problem, boundary conditions are dependent on x and y only
despite the presence of ﬁctitious time in the describing equations. Nevertheless, the
boundary conditions are still given by Eqs. (22)-(24) in terms of the velocities and
by Eqs. (37)-(40) in terms of the vorticity and the stream function; one exception is
the plane entry flow problem where the boundary motion, given in terms of UB in Eq.
(24), vanishes. The inclusion of steady motion of the entire lower boundary in the
direction of the ﬂow give rise to a Couette-Poiseuille ﬂow problem.

112mm! (221mm

I=0,V=0 n=6y~6y2,v=0
2

v=1oC=-—g-§- V=3yz-2y’.c=12!-6
y

 

 

 

i=1 v=0 0
’ I: v=0
V=Y9C= {

Figure 4. Dimensionless boundary conditions.

49
In the presence of a ﬁnite step, as shown in Fig. 5, boundary conditions on the
lower boundary must be modiﬁed. Speciﬁcally, boundary conditions in C and w on the
step sides facing upstream and downstream must be accounted for. While ‘I’ on the
entire lower boundary remains equal to zero, as given by Eq. (30) for a straight lower

boundary, the vorticity boundary conditions on the upstream and downstream faces of
the step becomes

§(an9t) = T 933%", x = X19x2;05ySB (41)
X

since Nut/3y2 is equal to zero on these faces. Vorticity boundary conditions on the

lower boundary segments both before and after the step as well as on the top surface

of the step remain to be represented by Eq. (35).

 

 

 

Figure 5. The ﬂow problem with a ﬁnite step placed on the lower boundary.

50
To further examine the effect of the inlet velocity proﬁle on the downstream
development of the steady entry ﬂow, an imaginary upstream extension to the
channel walls, in the form of streamlines may be assumed. Taking the velocity
distribution at the inlet section (x = O) to be unknown, the following boundary

conditions are imposed in the region x < 0:

u(-H,y,t) = 1, v(-H.y.t) = 0 (42)
$6.0.” = 0, v(x,0,t) = 0 (43)
3y
113,14) = 0, v(x,l,t) = O (44)
3y

where the uniform ﬂow initiates upstream at x = -H instead of directly at the inlet.
This is referred to as a cascade condition as it resembles an inﬁnite cascade of
equally-spaced plates in a uniform oncoming stream. In terms of the vorticity and the

stream function, the above boundary conditions in the assumed region upstream of the

inlet,become
W(-H.y.t) = y. C(-H.y.t) = 0 (45)
v(x.0.t) = 0. C(x.0.t) = 0 (46)
v(x,l,t) = 1, C(x,l,t) = o. (47)

Moreover, an inlet condition representing an experimental velocity distribution
obtained from hot-wire measurements is used to compare the development of the

streamwise velocity along the channel to that obtained using a uniform inlet velocity.

THE NUMERICAL METHODS

The main purpose of this investigation is to provide a detailed description of the
velocity ﬁeld by solving the nonlinear equations of motion. The ﬁnite-difference
approximation, which is commonly used in solving partial differential equations, is

selected to solve the describing equations numerically.

A rectangular computational grid is constructed with a mesh size carefully
selected since the total number of nodal points is limited by the available computer
capabilities. The describing equations, expressed in terms of the vorticity C and the
stream function \y, are ﬁrst approximated by ﬁniteedifference equations. In ﬁnite-
difference form, the vorticity equation is approximated using the Alternating-Direction
Implicit method while the stream function equation is approximated by the iterative
Successive Over-Relaxation method. The present numerical solution is second—order
accurate in its entirety since the numerical methods and boundary conditions have
second-order accuracy. The objective is to compute the vorticity and the stream
function at all the grid points for each time step. The two methods used in computing
our time-dependent, viscous, incompressible, plane, ﬂow problem are presented and

discussed in details in the following sections.

51

52

4.1 The Computational Domain

For this finite-diﬂ‘erence scheme, a rectangular grid in the xy-domain, with grid
lines parallel to the x- and y-axes, is superimposed to exactly ﬁt the geometry of the
ﬂow problem including the moving boundary. As shown in Fig. 6a, the channel is
divided into a ﬁnite number of increments such that the total number of nodes
(including boundary nodes) is M in the x-direction and N in the y-direction; this
results in (M-l) equally sized intervals in the x-direction and (N-l) equally sized
intervals in the y—direction. Each node denoted by the subscripts ”i, j " represents a
point (x,y) in the flow located at x = (i-1)Ax and y = (j-1)Ay where Ax and Ay are
the mesh size in the x- and y-directions, respectively. The inlet is taken to exist at i
= l, and the downstream condition is assumed to exist at i = M; the lower and upper

boundaries are placed at j = 1 and j = N, respectively.

Figure 6b illustrates the ﬁve nodes (x,y), (x+h,y), (x,y+h), (x-h,y), and
(x,y-h) with h being the mesh size, equal in the x- and y-directions (h = Ax = Ay is
to be used later in writing the difference equations); this representation for a typical
interior node is used throughout our work to express the derivatives of the describing

equations in discrete form.

4.2 The Alternating-Direction Implicit Method for the Vorticity Equation

After a thorough search of the literature, it was decided that the Alternating-
Direction Implicit (ADI) method would be used to solve the vorticity equation; the
reasons are outlined below. The ADI method was ﬁrst introduced simultaneously in
two fundamental papers by Peaceman and Rachford (1955) and Douglas (1955);

53

,1 _,| I, |,_ noun

 

I 2 i-l 1 1+1 M-l M x

a) ﬁnite—diﬂ'erence grid structure and ﬁeld boundaries.

 

 

(Kym) . i9j+l
(It-ha) (Xi!) (XHIJ)
. . .
i-l,’ Li “I”
(xix-h) 0 hi1

b) ﬁnite-difference cell and notation for a typical interior grid point.

Figure 6. Finite-difference representation of the computational domain.

54
further details were given by Douglas and Rachford (1956). The principles, variants,
and improvements of the method were also discussed in numerous articles such as
Birkhoff and Varga, ( 1959), Birkhoff et al. (1962), Mitchell and Fairweather (1964),
Fairweather and Mitchell (1967), and Lynch and Rice (1968); a generalized form can
be found in Douglas and Gunn (1964). A general description of the method is
presented in many texts including Amos (1977), Mitchell (1969), and Varga (1962).
Moreover, the method of Fractional Steps (or method of Splitting), which is
somewhat similar, and in some cases identical, to the ADI method, can be found in

Yanenko (1971).

In addition, early applications of the ADI method in the numerical solutions to
thermal-ﬂuid problems were given in many publications, for example, Wilkes and
Churchill (1966), Samuels and Churchill (1967), and Aziz and Hellums (1967) used
the method in their calculations of natural convection. Kuskova (1968) used what he
called the method of Variable Directions in his solution of the ﬂow in a rectangular
cavity. The speciﬁc implementation of the method to the Navier-Stokes equations
along with detailed discussions can be found in Roache (1982) and Peyret and Taylor
(1983). Indeed the ADI method continues to be used in all forms of ﬂow calculations.
The ADI method is currently the most popular approach to solve viscous ﬂow
problems (Roache, 1982). Most recently, Caughey and Iyer (1989) used the ADI
method in their computation of the supersonic ﬂow of two-dimensional inlet ﬂow

ﬁelds, using Euler’s equation.

The ADI method is a multistep, multidimensional, second-order accurate
implicit scheme used to solve partial differential equations. It is characterized by

unconditional stability when used with the proper boundary conditions; unconditional

55
stability in the case when the Dirichlet-type boundary conditions are not time-
dependent can be expected. It also provides a reduction in computing time by a factor
of two or more (Roache, 1982) since it allows larger time steps than the explicit
methods. In the ADI method, each time step is split into two halves resulting in two
ﬁnite-difference equations each of duration At/ 2 where At is the time step. The ﬁrst
equation is implicit only in the x-direction while the second equation is implicit only in
the y-direction. Equations are used in turn over successive time steps. Moreover,
the ADI method permits the construction of an efficient and implicit scheme that
requires only the inversion of a tri-diagonal matrix; full implicit schemes yield a tri-

diagonal system of equations only for one-dimensional problems.

In the x-equation, all implicit terms associated with x-derivatives are evaluated
at t = (k+l/2)At, where k denotes the number of time steps, while all terms

associated with the y-derivatives are evaluated at t = k At. In the y-equation,

however, all terms associated with the x-derivatives are evaluated at t =
(k+1/2)At, while those implicit terms associated with the y-derivatives are
evaluated at t = (k+1)At. Using the notation Cit.) to represent C(x,y,t), the two ADI

equations for an interiornode (i = 1, 2,..., M-1;j = 1, 2,..., N-l) in discrete form are

 

 

 

(k+1/2) (1!)
St ‘11 ' Cir
At/2
(k) (k) (hr/2) (k+1/2) (k) (k) (k) (k)
Vi.j+r ' V1.14 Ci+r.j ' Curd _ i+1,j ' ‘Vi-r.j Ci.j+r ' Cid-1
2Ay 2Ax 2Ax 2Ay

(43)

 

 

(RH/2) (k+l/2) (k+1/2) (k) (k) (k)
l [ Cm; ' 2911.1 + Ci—Lj + Cider ' 29m + Cigar ]
2
Re (Ax)2 (Ay)

(k+l) (k+-1’2)
C - C

 

 

 

 

 

id in‘
3‘ At/2
(k+l/2) (k+1I2) (k+1I2) (k+1/2) (k+1/2) (k+ll2) (k+l) (k+l)
+ i,j+l ' i,j-l Ci+1,j ' Ci-IJ _ Vi+l¢ ' Wing Ci,j+1 ' Cid-1
2Ay 2Ax 2Ax 2Ay
(ht/2) (k+1/2) (Int/2) (k+l) (k+1) (k+l)
= 1 Cm.) ' 29m +Ci-l,j + Ci.j+r ' 29rd + Cid-l . (49)
Re (Ax)2 (Ay)2

In Eqs. (48) and (49), the coefﬁcients of the nonlinear terms, namely the

velocity components u and v, are given in terms of the stream function using the

deﬁnition of Eq. (14). In difference form, they are given by

 

 

urn- = mi; 3"“ (50)
Vid’ = Vi'l'jzaxmﬂ'j ' (51)

All derivatives in the above equations are expressed using second-order
central—difference discretizations. In computing C i, j at a new time increment, ‘Vi.j and
its ﬁrst-order space derivatives are considered taken from the previous time step and
thus treated as known coefﬁcients. The vorticity equation then becomes a linear,
parabolic difference equation for C. This effectively decouples the vorticity and the

stream function equations.

The ADI method has a formal truncation error of o [ (1102, (Ax)2, (Ay)2]

(Anderson et al., 1984). Ideally, in order to preserve the second—order accuracy in

. . . 1 .
time, the coefﬁcrents of the nonlrnear terms must be evaluated as “if; I2) and via? 1n

57

it: 1/2) (k+1).

and vi]-

Eq. (48) andu 1n.Eq (49). Since u and v are given in terms of w,

which is calculated from Poisson’s equation, this evaluation procedure requires the
coupled solution of the vorticity and stream function equations at both time steps

(k+1/2) and (k), which is impractical. Several alternatives and their associated error

were discussed by Roache (1982). When the old values do?) and v."‘.’ are used In both

141-]

equations, the formal accuracy deteriorates to O[At, hz] (Son and Hanratty, 1969).

Briley (1971) calculated (if j)and vf" and linearly extrapolated forward tou 9‘1“”) and
1 -
via: I2) using stored values of “(1:11) and v3?) ; the procedure was stable and second-

order accurate but requires additional storage. The second—order accuracy is also
retained (Roache, 1982) when the stream function derivatives are reevaluated after

(HI/2)

the ﬁrst half of the time step by solving the stream function equation for Vi j . The

“(ht/2) (k+1/2) .
and v in

nonlinear coefﬁcients then become u:k j)and vfk j )in Eq. (48)an U M

Eq. (49).

Using a uniform mesh with spacing h = Ax = Ay, Eqs. (48) and (49) become

St— 2( C(k+-1’2) -

(k l/2)_ C(k 1/2)
+(‘l’i,j+l Vi,j 1)(ci+l,j C131 ) (vi+)1,j Vi- (k1,j) (cm i,j+l Cid-1

(k+1/2)_ C(k+l/2)+ (k+l/ 2) (k)
[ ( c1341 Ci- 1 J )+ ( C1j+l- Cid) +ci,j-l ] (52)

58

(k+1) C(kJ+l/‘2)
St— 2.,(t

(k+1l2) (RH 12) (k+1/2)_c(k+1/2) (k+ll2)_ WGHIIZ) (k+1) Cowl)
+ (VLj-t-l V111 )(Ci-l-IJ Ci- 1,j ) (Wi+l,jVi-1j )( Ci ,j+1 ci .j-l

(k+lf2)_ C(kJH/Z) C(kH/Z) (k+1) C-(km (1m)
—R4—c-[( Ci‘i'loj Ci- -,lj )+ +(Ci .-j+1 + Cid-1 ] e (53)

(HI/‘2)

The ﬁrst equation is solved for the intermediate value Cij (assuming all

variables at time step (k) are known) which is, in turn, used in the second equation;

this intermediate value can be considered as an approximation of the exact solution2

(Peyret and Taylor, 1983). The second equation then gives the solution Cg?” at the

end of the time interval. Each of Eqs. (52) and (53) represents a system of linear
algebraic equations which leads to a tri-diagonal matrix. The two sets of
(M-2)x(N-2) equations are solved simultaneously at each time step. A number of
suitable direct and iterative methods are available for solving such systems.
Tri-diagonal systems, however, can be readily solved using the Thomas algorithm,
which is a simpliﬁed adaptation of the Gaussian elimination method in the form of a

recursion formula.

 

2 lnusingguraalmuldstepmethodadnquesdonofhowwobmindrevaluesofthe
intennediatequantitiesontheboundaryarises. FortheADlmethod,inwhichtheequationusedin
eachnepiswnsiswntwimdwmighmlequadon,ahrgemcadmamrwuhmspeamdme
appears fa points adjacent to the boundary (Fairweather and Mitchell, 1967, Mitchell, 1969) if the
exacttime-dependentboundaryvalueisusedtodeterminetheprovisimalvalunggl’z)onthe
boundariesparallelmthey-axis,requiredtosolveEq.(52). Onewaytoovercornethisdifﬁcultyis
todeﬁneaboundaryvaluesuchthatthecombinationofbothtimestepsyieldssecond—orderaccuracy
atpointsadjacenttotheboundariesaswellasirrnerpoints. Thiscanbeaccomplishedifthevalue
C931”) atthcbomdariesisdeducedfromaveraging its valueattimcsteps(k)and(k+1). Other
alternatives are presented in Peyret and Taylor (1983).

59
Due to the diagonal dominance of the coefﬁcient matrix, the ADI method can
signiﬁcantly improve stability at high Re and, hence, extend the domain of
convergence. Nevertheless, time—dependent boundary conditions can destroy the
unconditional stability of the ADI method (Peyret and Taylor, 1983). In this case, a

limitation on the time step is necessary to ensure the condition of diagonal dominance

of the coefﬁcient matrix for Cid" The diagonal dominance is preserved if we have luid-l
S 2/(Ax Re) in Eq. (48) and lvid-l S 2/(Ay Re) in Eq. (49). If these conditions do not

hold, however, the constraints on the time step are given by (Peyret and Taylor, 1983)

 

2112 2112
At = . At = _
luiJle " 2/RC lvi’ley - 2/RC (54)

 

These satisfy the condition of diagonal dominance. The effect of the boundary
conditions on the stability of the method and the coupling between the vorticity and
the stream function on the boundary has been studied by Bontoux et al. (1980).

Although the limitation on the time step is not very restrictive at low Reynolds
number, for ﬂows with moderate and high Reynolds numbers the stability criteria
require such extremely small time steps and/or mesh size that the computation
becomes costly; the maximum magnitude of the eigenvalues exceeds unity so that the
von Neumann stability condition cannot be satisﬁed for any At/h. To remedy the
situation, the ﬁrst-order derivatives of C in the nonlinear convective terms in the
vorticity equation are represented with forward/backward differences instead of the
central differences used in Eqs. (48) and (49) (Forsythe and Wasow, 1960). The
resulting scheme is stable. Consequently, the second-order accuracy of the central-

difference approximation in the space variables is lost since the forward] backward

60
differences have only ﬁrst-order accuracy. By correcting the signs of the nonlinear
terms, this technique maintains the diagonal dominance of the system of equations at
high Reynolds number and, hence, helps avoid the instability of the numerical
solution. Schemes using such an approximation are found in the works of Barakat and
Clark (1966), Torrance (1968), Greenspan (1969a); also, Lilly (1965) used this

technique to stabilize his numerical solution of geophysical problems.

The technique, ﬁrst suggested by Lelevier (Richtrnyer and Morton, 1967),
requires the use of backward differences for positive velocity components It and v
(expressed in terms of derivatives of v) and the use of forward differences for
negative velocity components. In so doing, the stability criterion becomes weakly
dependent on Reynolds number (Pao and Daugherty, 1969). The stability criterion for
this forward] backward difference approximatioh is

4 lvi.'l [UL-I :1
—T + —-'— + —-l— At 5 l.
[ Reh h h (55)

The second-order accuracy is recovered at convergence (Peyret and Taylor,
1983), however, by alternating the direction of the forward] backward differences,
representing the first derivatives of C, at each time step (according to the signs of the
coefﬁcients of the nonlinear terms, i.e., the velocity components u and v). This is
achieved by considering upwind differencing for the implicit part and downwind
differencing for the explicit part. During the transient stage, the scheme is only ﬁrst-
order accurate in time and at steady-state the truncation error becomes of 0[ (At h),
hz]; this technique gives a centeral-difference approximation at convergence.

Implementation of the scheme associated with this upwind/downwind non-centered

61
approximation requires the replacement of the following differences in Eqs. (52) and
(53):

(CH-1,j ' Ci-lJ) = (1' a)(§i+1,j ' Cu) + (1 '1'“) (Cid ' CL”) (56)

(Ci.j+l ' C1,j-l ) = (1" B)(Ci,j+1 ' Chj) + (1+B)(Ci,j ' Cidq) (57)

where a = sign(ui.j) and B = - sign(vi.j) for the x-equation, and a = - sign(ui,j) and
B = sign (Vid) for the y-equation. The unconditional stability is preserved (Peyret
and Taylor, 1983) so that at each time step, the matrix is diagonally dominant. The
above technique was proposed by Peyret (1971) for the solution of the incompressible
Navier—Stokes equations in the primitive variables, and was applied by Daube and Ta

Phouc Lee (1979) and by Ta Phouc Lee and Daube (1981) for the vorticity-stream
function formulation.

Finally, the solution to the steady-ﬂow problem is obtained as the asymptotic
solution for large time of the unsteady equations (thus, the term ac/ax is retained in
the vorticity equation). Time steps are ﬁctitious and each is merely considered
equivalent to an iteration. The use of this pseudo-unsteady method has been
introduced by Burggraf (1966) in association with an explicit-implicit discretization

with respect to time.

4.3 The Successive Over-Relaxation Method for the Stream Function Equation

A commonly used technique of solving a system of simultaneous equations is

the iterative Successive Over-Relaxation method. It has been successfully utilized in

62

the solution of many important problems of mathematical physics including the
subclass of elliptic, linear, two-dimensional partial differential equations, such as the
Laplace and Poisson equations with Dirichlet or N eumann boundary conditions on a
closed contour. The method, introduced independently by Frankel (1950) and Young
(1954), is thoroughly described by Forsythe and Wasow (1960), Greenspan (1965),
Wachspress (1966), and Varga (1962). General applications of the method to ﬂuid
ﬂow problems can be found in Chow (1979) and Arpaci and Larsen (1984).

Mostly, the SOR method is desirable since, due to its simplicity, it can be
readily implemented. Furthermore, it does not require large computer storage, which
is a critical problem when a large number of nodes are involved. In the SOR method,
however, inner iterations are required at each time step. An initial guess at the
solution is improved with subsequent successive approximations until the solution
converges to within speciﬁed accuracy after a ﬁnite number of iterations. The stream
function for iteration (k+-l) is combined with the kth iteration using a relaxation
factor. The use of a relaxation factor reduces the number of iterations required to
achieve a solution, and therefore, improves the convergence behavior of the SOR

method over other iterative procedures.

The stream function equation, given by Eq. (15), in ﬁnite-difference form is

 

Yum ' 2‘Vi.i+ Vi-l,j + Vi,j+l ' 2‘1’i.j + Vij-l
(AK)2 (Ay)2

 

= - Cr,- (58)

where second-order central-difference approximations are employed to approximate

the second-order derivatives at all interior mesh points with an error of O[(Ax)2,

63

(Ay)2]. With the vorticity Cid. considered as a known constant, the above represents

a linear, elliptic partial difference equation for Win"

In the explicit point-iteration method (i.e., the Jacobi method) we march through

the grid from node to node solving Eq. (58) explicitly for the new value “(31m using

old values VIE-l) at adjacent nodes known from the previous iteration; here (m)
indicates the number of iterations. This is often referred to as Richardson’s iteration.
To speed up convergence, as in the Gauss-Seidel method, old values of adjacent
nodes are replaced by new ones as soon as they become available. Rearranging Eq.

(58), and using a uniform mesh of size h, we obtain Liebmann’s iterative formula

(n+1) = _1_ (m1) (m) (n+1) (m) (m)
U 4 1: “J + Vi+l.j + Vin-1 + Vi.j+1 + h2§i.i :1 (59)

As a reﬁnement to the above methods to further improve the speed of
convergence, the SOR method employs the correction formula in which the stream
function for iteration (m+1) is combined with its value from iteration (m) using the
over-relaxation factor FS

(if? 1’ = FS mg?” + (1 - FS) (if? (60)
where (riff 1’ is a relaxed value of the new value Vi?” obtained from Eq. (59); this

relaxed value is an adjusted (or extrapolated) value to be used in lieu of V331“)

explicitly in the next iteration; vi? is the value from the previous iteration as

modiﬁed by a previous application of this formula if over-relaxation is being applied

64
successively at each iteration. In the above equation the over-relaxation factor F S
controls the rate of convergence. It has a value 1 S FS S 2 and must be optimized for

a given computation in order to maximize the rate of convergence (i.e., reduce the

number of iterations) of the SOR method.

For Poisson’s equation, the optimum value of FS depends on the mesh size, the
shape of the domain, and the boundary conditions. For the present Dirichlet problem
in a rectangular domain with uniform grid of size (M-1)h by (N-1)h, it has been
shown (Chow, 1979) that the optimal value of the relaxation factor is calculated from

 

 

the analytical expression
_ 2 (l -/ 1- §2 )
F809, _ g2 (61)
where
g :3 COS[I/(M-l)] :- COSIE/(N‘ln (62)

With the determination of the optimum over-relaxation parameter FS, it is
possible to reduce the required number of overall iterations in the solution by more
than one order of magnitude from that required by Gauss-Seidel iteration (Thompson,
1969); in addition, we may reduce the restriction placed on the maximum size of the
spatial step imposed by the Gauss-Seidel technique. Utilizing Eq. (60), we can then
rewrite Eq. (59) as

65

.. (m+1) _ _F§_ (m+1) (m) (m+1) (m)
Vii ‘ 4 [ i-1.i + i+l.j + ‘Vioi-l + Vi.i+l + “2 Cid]

+ (1 - Fs) Vi? (63)

The above is an iterative method which improves the approximation one point at

atime.

The iterative SOR by blocks (or lines) method is used to implicitly solve Eq.
(15) for the stream function. The time of computation can be reduced by iterating on
an entire block of points; a block represents a column or a row of node. For the
present inner iteration (m+1), the stream functions for a block of nodes, considered to
be columns, are determined from the values at adjacent nodes. At each column i = 2,

3,..., M-1, Eq. (40) is rewritten as

W331) _ 4Worm) + (n+1) = - [Won-+1) + “[010 + h2 Cm." ] . (64)

i,j i,j+l i-l,j i+l.j i,j

. . l l 1
A column of nodes wrtlr stream functron values vi? ), v(im; ),. . ., Vim; ) for a

given i constitutes a system of (N - 1) simultaneous linear algebraic equations with a

tli—diagonal matrix. One iteration for each column i will implicitly solve for WES-1+1)

values at all nodesj = 2, 3,..., N-l. This procedure is repeated for all i = 2, 3,..., M-

1 , one column at a time, successively updating known Vi? values from the previous

+1)

I) Relaxed values ‘17:“;

iteration by new values Vi“? . are computed using Eq. (60)

l
as soon as tug-H) become available. These relaxed values are used in the

computation at the next adjacent column and in the subsequent iteration as well.

66
Solution of the system of equations associated with columns i = 2, 3,..., M-1, yields

the values of the stream function at all interior nodes.

This approach is expected to speed convergence since the implicitness makes
boundary conditions at j = 0 and j = N felt at all interior nodes along a given column
(Arpaci and Larsen, 1984). The old incorrect values in adjacent columns and the
added complexity of the solution of the tri-diagonal matrix, still slow the
convergence. Ames (1977) indicates that, compared to the Gauss-Seidel iterations,

the SOR by lines method requires only 1N2 as many iterations to reduce the initial

errors by the same amount.

In the SOR method, the solution is improved with successive iterations until it

converges to within a speciﬁed accuracy after a ﬁnite number of inner iterations.
1

Except for the points where Vita)“ ) = 0, the iteration process is terminated as soon as

the relative error criterion

(n+1) (m)

. . - \ll. .
1.1 1.1
(m+l) S E (65)

i,j

 

is satisﬁed at all interior nodes, where e is a preassigned tolerance error.
Conversely, in order to save computing time, if the above is not satisﬁed, the process
shall be terminated when a preselected maximum number of iterations is attained,
even if the results have not yet reached the desired degree of accm'acy. In general,
the maximum number of inner iterations required for convergence is directly

proportional to the square of the number of equations to be solved (Roache, 1982),

i.e., mm ~ [(M-2)x(N-2)]2.

67

4.4 Initial and Boundary Conditions

In order to solve the difference equations for the vorticity and the stream
function, initial and boundary conditions on C are needed for Eqs. (52) and (53) and on
\y for Eq. (64) (or Eq. (59)). The initial conditions are given in Eq. (26). The
boundary conditions are given by Eqs. (28)-(31) for \y and Eqs. (33)—(36) for C. The

ﬁnite-difference form of these initial and boundary conditions are presented herein.

First, the initial conditions (at time t = 0) on the stream function and vorticity

canbeputintlreform
‘l’i,j = (I'Dh. chj = 0° (66)

In the case of the stream function, stating the boundary conditions in ﬁnite-

difference form is straightforward. At all time levels, these boundary conditions are

asfollows:

upstream: \le = (j-1)h (67)
downstream: VMJ = “(1'1)“2 - 2[(.i-1)h]3 (68)
lower boundary: (rm = o (69)
upper boundary: ‘Vi.N = 1. (70)

Likewise, the vorticity boundary conditions on the upstream and downstream

boundaries are directly written as

upstream: CLj = 0 (71)
downstream: CM.j = 12[(j-1)h] - 6. (72)

In addition, an explicit vorticity boundary condition is necessary to replace the
second-order derivatives in Eqs. (35) and (36). Vorticity is generated if there is a
no-slip condition at the boundary; thus it must be incorporated in the numerical
solution. Expressing these conditions in ﬁnite-difference form is complicated since
the no—slip boundary conditions are difﬁcult to handle computationally (Pearson, 1964,

1965a). The motion of the lower boundary also enters these conditions.

The vorticity at each boundary point is determined by expanding the stream
function at the node adjacent to the lower boundary in a Taylor series, in terms of the

stream function at the boundary, and keeping the third-order term, i.e.,

 

..—.. h_ay_ 1‘2 32‘! i131- H...73
“'2 “’14- (3Y)i.r+ 2 (3Y2 1.1+ 6 3Y3 i.r+ OT ( )

where derivatives of w on the boundary in the x-direction are zero. Neglecting terms

of higher order, and noting that for the moving section of the boundary (av/8y)“ =

ui.l = UB and that, from Eq. (15), (32W/3y2)i.1 = - C“, we get

2 3
Via = ‘Vi,r + 11 U3 ' “1'15— “ - 316—635.)“. (74)

Upon substituting the approximation

(ﬁll 3 Ci ‘Cir (75)

By h

into Eq. (74) and rearranging, we obtain the second-order (error is of O[h2])
"Woods’ condition" (Woods, 1954)

 

(k) (k) (k) (k+”
‘9‘“) __ 3(Vi,l ‘ Wm) _ i,2 3UB (76)
"1 ‘ h2 2 + h

which gives the vorticity on the lower boundary in terms of the stream function and
vorticity obtained from the previous iteration; the superscripts (k+1) and (k) denoting
the present and previous time steps, respectively, are added to indicate the time of
the available information; since the value of boundary vorticity depends on time, it is
recalculated at every time step. It is also clear that motion of the boundary enters
the problem through the vorticity boundary conditions. The velocity of the moving
portion of the boundary can be a constant as is the case with steady motion or a
known function of time expressing the oscillation of the lower boundary (U3 = 0 on
the ﬁxed portion of the boundary). Woods’ condition has been used in particular by
Bozeman and Dalton (1973).

An equivalent formula can be derived for the upper boundary, with does not
include boundary motion; that is

 

(k) (k) (k)
C(1+1) 3 (th ' VLN-l ) C._._12 (77)
i.N = h2 - 2 .

Another frequently used second-order formula for determining the vorticity at

the boundary was introduced by Jenson (1959). It is derived in a manner analogous

70

to Eq. (76). A linear combination between the expressions of the stream function at

two points adjacent to the boundary allows us to evaluate az‘wayz. The result is

 

(k) (k) (k) (k+”
(k+!) = 7V}; - 8WB+1 + VB+2 + ﬂL. (73)
B 2112 h

Here the subscript "B” represents a boundary point and "B+1" a neighboring internal
point in a normal direction, regardless of the orientation of the boundary. The same
expression can be obtained from approximating \y near the boundary by a third-order
polynomial. The above formula was used by Pearson (1964, 1965a) and Wilkes and
Churchill (1966). The stability of Woods’ condition was considered by Taylor (1970).

Briley (1971) considered a formula of the same type which gives very little
change in CB but converges faster since it permits larger time steps. However, extra
programing effort is required since the resulting matrix is not tri-diagonal. Briley’s
rational and derivation are reproduced in Roache (1982).

An extensively used ﬁrst-order condition for the boundary vorticity was
proposed earlier by Thom (1933). After ignoring the second-order term in Eq. (73)
and following a similar analysis to that leading to Eq. (77), the following expression

 

is obtained:
(k) (k) (k+1)
(k+!) _ 2 (VB ' VB-t-l + ZUB (79)
B ’ h2 h '

An alternative ﬁrst-order formula was derived by From (1963, 1964) to give

the boundary vorticity in terms of the stream function. Of course, expressing the

71

vorticity on the boundary in terms of the stream function can also be obtained from a

first- or second-order ﬁnite-difference representation of C = - azwayz.

A comparison of the errors associated with Eqs. (76) and (79) is presented in
Section 6.1. Because C at the boundary is calculated using previous values, it can lag
behind the values in the internal region by At (Roache, 1982). For large At, the
calculation can become inaccurate and an iteration procedure is required. An
alternative technique is the use of weighting factors to slow down the solution in the

interior region; this is described in Chapter 5.

Inthepresenceofaﬁxedsteponthelowerboundary,avorticityboundary
condition equivalent to Eq. (76) is used to replace Eq. (41) on the side surfaces of the

step. Also, for the problem with a cascade inlet condition, the condition Cid = 0 is

used at all boundary points having zero vorticity as required by Eqs. (46) and (47).

A closing remark concerning the vorticity boundary condition is in order.
Methods employing creation-type conditions were discussed by Anderson (1989). In
the methods based on vorticity creation a constraint, in the form of an integral, is
imposed on all boundary points. These constraints are equal in number to the number
of unknown boundary values and, therefore, are used to close the system of
equations. The value of the vorticity on the boundary is determined when the time
derivative of the constraint vanishes as the vorticity is evolved from the solution of
the vorticity equation. Anderson demonstrated that vorticity creation boundary
conditions, when used in numerical solutions, are generally equivalent to those which

exploit the relationship between the vorticity and the stream function.

THE COMPUTATIONAL SCHEME

The mechanics of the computational scheme adapted for this study is disclosed
in the following sections. The computational scheme, used in conjunction with
accelerating parameters, is comprised of the Alternating-Direction Implicit (ADI)
method and the Successive Over-Relaxation (SOR) method. A brief description of

the implementation of the scheme and the computing environment are also given.

5.1 Accelerating Parameters

Fume-difference approximations used when modeling physical problems often
lead to the solution of large-order sparse linear matrices. It has now been
established that central-difference approximations provide convergent and accurate
numerical solutions only for low Reynolds numbers. For large Reynolds numbers
(greater than 200), however, upwind differencing is essential to ensure the stability of
the numerical solution by maintaining the diagonal dominance of the system of

equations.

As a supplement to the numerical methods, the vorticity values obtained from

solving Eqs. (52) and (53) and the stream function values obtained from solving Eq.

72

73
(64), are modiﬁed by introducing certain weighting factors into the numerical scheme.
These factors are often referred to as accelerating parameters. Another type of
accelerating parameters is the time step for the vorticity equation. The over-
relaxation factor used in conjunction with the SOR method in the solution of the
Poisson-type elliptic stream function equation is also considered as an accelerating

parameter (Thompson, 1969).

In general, accelerating parameters are used to minimize the time of
computation. The use of accelerating parameters, however, have a signiﬁcant effect
on the convergence characteristics of the difference equations; they provide numerical
stability by maintaining the diagonal dominance of the system of difference equations,

hence, expand the domain of convergence.

The idea behind the introduction of weighted averages is to avoid any possible
divergence that results from any over-correction to the value of the vorticity or stream
function; the numerical solution is slowed down whenever over-correction occurs.
This procedure stabilizes the numerical solution by compensating for the coupling
between the vorticity and stream function equations. More so in the case of the
vorticity equation, the function of the weighting factor is to slow down the solution in
the interior region in order to reduce the lag in vorticity boundary conditions behind
the solution in the interior region, since the vorticity boundary condition is calculated

from the stream function and vorticity values of the previous time step.

Weighting factors are used to modify the values of C and v in the interior
region. They allow for different fractions of the old and the new values obtained from

subsequent time steps to be used during the matrix solution. Smoodrin g formulas are

74
introduced for this purpose (Greenspan, 1969a, 1974). The weighted averages of t;

and v at the end of each time step are taken as linear interpolations of the new
values calculated at the end of the present time step (k+1) obtained, respectively,
from Eqs. (53) and (64) and the weighted values of the previous time step (k); that is,

CE?” = (1 - a) Q?" + 0 cf]? 0 $051 (80)
4M) _ (k+» (k)
i,j — (l - 'y) Vi.j + 7 Vid 0 57s 1 (81)

where o and y are the weighting factors, and CE?” and ﬁg“) are the weighted

averages of the vorticity and the stream function, respectively. The values of the

above weighting factors fall within the range of 0 to l and are not necessarily equal for

the vorticity and the stream function.

A similar formula can be used to give a weighted average for the value of the
vorticity at the boundary (Peyret and Taylor, 198 3). However, the weighting factors

corresponding to the interior points and the boundary points are not necessarily equal.

The effect of the values of the accelerating parameters on the convergence of the
solution and the computing time depend on the shape of the computational domain,
mesh size, and the Reynolds number; in essence, it depends on the number of nodal
points in the computational domain. Although the objective of using these parameters
is well understood, the determination of their optimum values is not straightforward.
Optimum values of the accelerating parameters (At, FS, 0, and y), with the
exception of the over-relaxation factor F S for the stream function equation, are

determined primarily by numerical experimentation. The optimum value of FS used in

75
the iterative SOR solution of Poisson’s equation can readily be calculated from Eq.

(61). It should be noted that the stability and rate of convergence of the numerical

solution is dependent on the size of the time step At .

5.2 Outline of the Computational Procedure

The details of the computational scheme used to implement the numerical
methods, discussed in Chapter 4, for computing our time-dependent, viscous,
incompressible, two-dimensional ﬂow problem are illustrated in this section. The
objective is to compute the stream function \y and the vorticity C at all the grid points
for each time step by solving the governing equations in their ﬁnite-difference form.
The ADI method is used to solve the vorticity equation, while the SOR method is
used to solve the stream function equation. The x- and y-component velocities, as

well as the pressure gradients, are then calculated from the converged solution.

After discretizing the equations and boundary conditions, the computational
domain is deﬁned; it depends on the Reynolds number. The mesh size and the time
step are then chosen. The computation begins at t = 0 by assuming C8.) = 0 and WE?
= (j- l) h as initial values at all nodal points (Eq. (66)). Boundary conditions given
by Eqs. (67)-(72) are then calculated. The values of the vorticity {B on the lower
and upper boundaries are also computed using Eqs. (76) and (77). When boundary
motion exists, the velocity of the boundary is determined at the end of the time step

and used in the calculation of CB.

76

The vorticity C115) at all interior nodes is computed at the end of the ﬁrst time

step by solving Eqs. (52) and (53) sequentially using the ADI method and available

2)

values for the stream function VF}; new values for V3,! ate obtained from Eq. (64)

before solving Eq. (53) for the second half of the time step. Then, utilizing the known

values of the vorticity cf]? as constants, the stream function will.) is obtained by
solving Eq. (64) using the iterative SOR by columns method and employing the
relaxation formula of Eq. (60) (for the SOR by points method Eq. (63) is used). The
inner iteration process of the SOR method is continued until the assumed error
criterion, given by Eq. (65), is satisﬁed at which point the iteration is terminated; the
tolerance error e = 10" is used. If convergence does not occur for some preselected
maximum number of iterations, the computational procedure proceeds to the next time
step to compute new values for Cid and WM.

(1)

The values for cf}; and Win am used in turn to determine the boundary

conditions for the vorticity and the coefﬁcients of the nonlinear terms in the vorticity

equation for the next time step, k = 2. Subsequently, CS.) and WE? are determined.

This above computational procedure can be generalized and repeated in the same
fashion with general steps (k) and (k+1) as the old and present time steps,

respectively, to yield a new estimate of the vorticity Cit-+1) and stream function ng’tl)

at the end of each time step.

(k+l)

i,j

(k+1)

The computed values C M are corrected at the end of each time

and \y

step using the weighting factors. Incorporating the weighting factors to modify the

77
values of the stream function and vorticity in the interior region ensures convergence

of the numerical solution.

The computation continues as the computational scheme progresses in time.

The computational procedure ends as soon as the relative error criterion, given by

(RH) _ fit)

iii hi _ -6
fowl) S e — 10 (82)

i,j

f

 

selected as a condition for convergence, is satisﬁed at all interior nodes; nodes with

Chi = 0 and Via‘ == 0 are excluded. Here a is the preassigned tolerance error and f

represents either Cid. or Via" The numerical solution converges when Eq. (82) is

satisﬁed at all nodal points for two successive time steps. For the time-dependent
flow with an oscillating boundary, the above error criterion is enforced over an entire
time period before convergence can be assumed. The computational procedure is
terminated after a preselected number of time steps even if convergence is not

attained.

In summary, the numerical procedure is executed according to the following
steps:
1. Assume initial values for Ci.j and ‘Vid' at all interior nodal points as given
by Eq. (66).
2. Calculate the boundary conditions on C and W from Eqs. (67)-(72).
3. Calculate the values of vorticity CB on the lower and upper boundaries

using Eqs. ('76) and (77).

78
4. Successively for every interior nodal point:

a. Compute the vorticity Cid“ from Eqs. (52) and (53) using values for
via. from the previous time step. For the second half of the time step,
however, compute new values for V” by solving Eq. (64) using
values of the vorticity Ciu’ obtained from the ﬁrst half of the time step.

b. Compute the values of the stream function ‘Vid’ from Eq. (64) using

the already available values of CM as constants. Carry on the
iterative process for some preselected maximum number of inner
iterations or until a certain error criterion is satisﬁed.

c. Correct the values of (i,j using Eq. (80).

d. Correct the values of Wm using Eq. (81).

5. Continue to advance the solution in time until the relative error criterion of
Eq. (82) is satisﬁed at all interior nodal points.

6. Return to step (3) and repeat the process if this error criterion is not

satisﬁed.
The above computational procedure is illustrated by the flow chart of Fig. 7.

Finally, the velocity components u and v are calculated from the values of the
stream function using a ﬁve-point, ﬁnite-difference approximation of the deﬁnition of
the stream function, as given by Eqs. (50) and (51), respectively. The streamwise
pressure gradients of interest, at the wall and centerline, can also be calculated in a
similar manner; the pressure gradient can be obtained from Eq. (18) and then

expressed in ﬁnite-difference form.

79

 

New time step (iteration)

 

 

 

Figure 7. Flow chart of the computational procedure.

5.3 Computing Environment

Computer programs were written to carry out the computational procedure as
outlined in the previous section. The FORTRAN programing language is used with
double precision for the real variables. The code for the oscillating boundary problem
is reproduced in Appendix A. Since the algorithm is written in general terms, the
other codes are similar and, hence, will not be presented.

A cluster that includes a DEC VAX/VMS 11/750 mainframe computer with 8
megabytes of memory and three DEC MicroVAX-II minicomputers, each with 9

megabytes of memory have been used to perform the computations.

When the number of nodal points becomes quite large such as in the case of high
Reynolds numbers and small mesh size, both the number of operations at each time
step and the number of time steps required to achieve convergence increase
dramatically; the time of computation increases accordingly and the convergence
process becomes tediously slow. A larger computer would handle the computation
more quickly. Also, since vector processing can speed up the computation, a
vectorized computer program executed on a supercomputer or a machine with parallel

processors would signiﬁcantly reduce the time of computation.

ANALYSIS

The topics discussed in this chapter are related to the characteristics of the
numerical solution of the Navier-Stokes equations. In any numerical solution, it is
often useful to look into the stability and consistency of the solution, as well as the
error associated with it. The difference representation of the partial differential
equations of motion must meet the conditions of consistency and stability in order to
be considered acceptable; generally, this implies convergence. These issues will be
dealt with in the following sections. Furthermore, the ﬁnite-difference approximation
of the vorticity boundary condition is analyzed to determine its inﬂuence on the

numerical error; this will be addressed ﬁrst.

6.1 Error Estimate for Vorticity Boundary Conditions

It is often postulated that the overshoot in the velocity proﬁle near the boundary
is due to the order of approximation of the vorticity boundary condition. In order to
support the hypothesis that this overshoot is the outcome of an inadequate numerical

approximation, a simple error analysis is presented herein.

81

82
Let the velocity proﬁle be represented by the polynomial u(y) = a + by 4» cy2 +

dy3 where a, b, c, and d are arbitrary constants. At the centerline, the velocity is uCL
and its derivatives are (Bu/By)“ = O and (Eizu/Eiyz)CL as 0. Upon incorporating this

information into the deﬁnition of the velocity and solving the resulting system of

equations, we obtain the following values for the constants:
a =-- 0, b = 2uCL, c = -3uCL, d = “or; (83)

The error associated with the ﬁrst-order accurate vorticity boundary condition
(Thom, 1933) is considered ﬁrst. In the derivation leading to Eq. (79), the third-order

terms of the Taylor series expansion are neglected. The error’s leading term,

therefore, becomes proportional to ~(h3/3)(8CB/8y). From Eq. (12), we have CB = -

Bu/By since avlax = O at the no-slip boundary. It follows that at the boundary,

where y = 0, we obtain

 

3:: = -2c = on“. (84)

From the above analysis, it is obvious the that the error in the expression of CB

is ~ - 2 h uCL < 0. The injected vorticity is being overestimated by 2 h “c1.-

For the second-order accurate condition (Woods, 1954), however, the third-

order term in the Taylor series expansion is retained and the leading term of the error

is proportional to -(h2/8)(32CB/8y2). On the boundary, we require that

83

 

2
33528 _.._. '6d = '6“CL° (85)

The above approximation obviously underestimates the injected vorticity CB by

3/4 hzuCL since the error is ~ 3/4 h uCL > 0.

Indeed the overshoot can be thought of as being physically incorrect since its
presence cannot be substantiated by experimental data. Morihara and Chang (1973)
claimed that overshooﬁng is necessary to satisfy the continuity equation. Other
authors explained the phenomenon as a consequence of the accelerating inviscid core
or the large velocity gradients near the inlet. Others attributed these bulges to the
upstream diffusion of the vorticity which renders the use of a uniform inlet velocity as
incorrect; a cascade inlet condition was suggested instead. All previous analytical
analyses justifying the overshooting were either conducted for vanishing Reynolds

numbers or used in association with the boundary-layer approximation.

6.2 Consistency and Errors of a Numerical Solution

Consistency of a numerical solution is a measure of the extent to which the
ﬁnite-difference equations approximate the partial differential equations. A ﬁnite-
difference scheme is consistent if the difference between the partial differential
equation and its difference representation vanishes under grid reﬁnement. This
difference is deﬁned as the truncation error. The truncation error associated with any
numerical approximation is the sum of all terms not included in the Taylor series
expansions used in the discretization of the derivatives. The truncation error is

considered acceptable if the difference approximation is consistent and stable.

84
The difference between the exact and numerical solutions of any partial
differential equation is attributed to the error associated with the solution of the finite-
difference scheme. This error is the sum of the discretization error of the ﬁnite-

difference approximation and round-oﬂr error of the calculations.

The discretization error is caused by replacing the continuous derivatives by a
discrete representation. This is deﬁned as the difference between the exact solution
(even though an exact solution may not be available) of the partial differential
equations and the solution to the difference equations. This deﬁnition includes the
truncation error of the difference representation of the partial differential equations
plus the error introduced form the ﬁnite-difference approximations of the boundary
conditions.

The round-off error results ﬁ'om rounding the number of digits in the arithmetic
operations. The round-off error is proportional to the number of arithmetic operations
and may depend on other factors such as the number of nodal points of the
computational domain; in such a case, while the truncation error may be reduced from
grid reﬁnement, the round-off error will increase. The round-off error also depends on
the machine accuracy and the precision of the computational algorithm; these are

related to the number of signiﬁcant digits carried during the calculations.

One additional source of error is the tolerance error imposed by convergence
criteria for the termination of the computation. In all numerical solutions, certain
preassigned error criteria are used to decide whether additional computation
(iterations for steady problems or time steps for marching methods) is necessary. In
this study, these criteria are represented by Eqs. (65) and (82).

85

6.3 Numerical Stability of the Differenced Navier-Stokes Equations

Numerical stability is a subtle issue closely related to the convergence of the
numerical solution. It is related to the error of the ﬁnite-difference approximation. A
numerical scheme is considered stable if the errors do not grow as the solution
marches in time. It is not easy, however, to formally establish whether a certain
numerical scheme is stable or not. Restrictions limiting the size of the spatial
increments and/or time step are often required to avoid instability. If these
restrictions are not met, the numerical results exhibit oscillations that grow very
rapidly, and after a few time steps their amplitude becomes inﬁnite so that an
"overﬂow" is registered by the computer; this phenomenon is characteristic of
instability.

The convergence of a numerical solution is understood to mean that the ﬁnite-
difference solution tends to the exact solution of the boundary-value problem as the
increments in space and the time step are reduced in size. Although a certain ﬁnite-
difference approximation may be consistent, the solution must be stable so that it has
a chance of converging. A proof of this principle applied to a well-posed, linear, initial-
value problem is given by the Equivalence Theorem due to Lax (Richtmyer and
Morton, 1967). Lax’s theorem states that when the consistency condition is
satisﬁed, stability is the only necessary and sufﬁcient condition for convergence. On
the other hand, the theorem has never been proven true for a nonlinear partial

differential equation such as the Navier-Stokes equation.

The present numerical scheme requires the simultaneous solution of the

vorticity and stream function equations. It is common knowledge, however, that

86
Poisson’s equation possesses excellent convergence properties when solved alone.
Therefore, to address the stability problem of the numerical solution, it is reasonable
to assume that the more likely source of any instability is the nonlinear vorticity

equation.

Each of Eqs. (52), (53), or (64), represents a general system of N linear

algebraic equations in N unknowns. In matrix form, they are written as
[allitl = [bl - (86)

where [a] is the square coefﬁcient matrix of order N xN, [q] is the solution column

vector (containing the unknowns C or v), and [b] is the column matrix representing

the known quantities on the right-hand side.

In the case of the elliptic stream function equation, let us assume that the
system of equations, given by Eq. (64), is irreducible (cannot be arranged such that

some of the unknowns can be determined by solving less than n equations). If we

further assume that [a] is nonsingular then the solution [\v] exists and is unique. By
ordering the system of equations such that the coefﬁcient of the unknown at each row
(or equation) is largest in magnitude, then the leading diagonal elements of [a]

become largest in magnitude. It turns out that when the elements aij of the coefficient

matrix satisfy the condition

i=N

laiil 2 2 Iaijl (87)

i=1
jti

87

for all i, in general, except at points adjacent to the boundaries where it applies to a ij

not multiplying boundary values of v, then the solution converges as the matrix is
assured to be diagonally dominant (V arga, 1962). The strict inequality (the greater
than sign) must be true for at least one value of i and will normally be satisﬁed at

points adjacent to boundaries.

The interpretation of this sufﬁcient condition is clear in light of the required
diagonal dominance of the coefﬁcient matrix. For each equation, it requires that the
magnitude of the diagonal elements be greater than the sum of magnitudes of all other
non-diagonal coefﬁcients with the inequality applied to at least one equation usually
near the boundary. The diagonal dominance imposed by the above condition is forcing
the error, which is deﬁned as the difference between the calculated value and the
exact solution, to become smaller and smaller as the iteration is repeated cyclically.
Since the above is only a sufﬁcient condition, convergence may sometimes be attained
even though the above condition is not met. A necessary condition shall not be
stated since it is impractical to evaluate (Anderson et al., 1984). The above stability
condition also holds for all individual equations resulting from Eq. (63) for the SOR by

points method.

The stability requirements for the system of parabolic linear difference
equations, resulting from Eqs. (52) and (53), can be established from the von
Neumann (or Fomier') analysis (Anderson, et al., 1984) for marching methods applied
to a linearized system of equations. For stable ﬁnite-difference calculations, the
maximum eigenvalue of the ampliﬁcation matrix associated with the system has a

magnitude less or equal to one. This leads to the requirement that

Ilsa-955- $1 (88)

where A.“ is the largest eigenvalue of the Jacobian matrix of the system. The above

condition ensures that the error continue to be bounded as the solution progresses in

time.

It should be noted that the von Neumann stability analysis presented above
does not include the effect of boundary conditions even though a matrix notation for
the system is used. The inﬂuence of the boundary conditions is not easily included in
the analysis. However, the matrix method can be used for one-dimensional linear
equations to demonstrate the effect of boundary conditions on the stability (Anderson
et al., 1984). In reality, the vorticity boundary condition puts restrictions on the time
step, in contrast to the unconditional stability inferred from von Neumann’s condition

(Roach, 1982).

Situations can occur, practically at large Reynolds numbers, when the
convergence condition is violated at some points. This has led to an alternative
formulation in which the diagonal dominance is maintained by approximating the
convective terms with schemes of lower-order accuracy. This leads to coefﬁcients
that satisfy Eq. (88) under all circumstances (Dennis and Chang, 1969). The higher-
order accuracy can be preserved by the correction technique of Eqs. (56) and (57).

A comprehensive analysis of stability, including other methods of analysis with
many theorems and proofs, is contained in Richtmyer and Morton (1967).

DISCUSSION OF RESULTS AND CONCLUSIONS

The subject matter of this chapter is divided into two major topics. The results
obtained from the present numerical study are presented. The flow problems, in
addition to the role of the accelerating parameters, are discussed separately. The
related ﬁgures are displayed at the end of the chapter. Based on these results,

conclusions are then drawn.

Convergence is attained by advancing the solution in time until it no longer
varies signiﬁcantly. An error criterion is used to decide whether additional time steps
are necessary; a tolerance error e = 10" is selected as a condition for convergence.
The transient solution is of no interest (time steps are ﬁctitious for the steady
solution). The present numerical solution is second-order accurate since both the

numerical scheme and boundary conditions have second-order accuracy.

7.1 Entry Flow in a Plane Channel

The starting point for this analysis is the consideration of the full Navier-Stokes
equations in terms of the vorticity and Poisson’s equation for the stream function,
describing the viscous, incompressible ﬂow in the entry region of a two-dimensional

89

90

channel. The vorticity C, stream function \v, streamwise velocity u, and normal
velocity v are calculated at each nodal point for a uniform inlet velocity proﬁle and
Reynolds numbers ranging between 10 and 5000. A uniform grid with a spacing size
of 0.05 is used. The number of meshes in the x-direction is equal to Re/ (20 H). Run
CPU times on our available computer ranged from 2.5 minutes for Re = 20 to
approximately 12 hours for Re = 2000 using the full channel height and optimum
values of the accelerating parameters. Due to symmetry, the number of nodal points
of the computational domain is reduced by considering one half the height of the
channel. Consequently, the required CPU time is reduced to 0.5 minutes for Re = 20

and nearly 2 hours for Re = 2000.

The two regions making up the entry ﬂow (as shown in Fig. 2) are identiﬁed as
the inviscid-core region, in which a viscous layer is assumed to exist near the wall,
and the proﬁle-development region, in which viscous effects inﬂuence the entire
channel ﬂow. The inviscid-core region ends when the thickness of the viscous layer
becomes equal to half of the channel height. This is determined numerically to occur
when the velocity at the centerline becomes larger than the velocity at the ﬁrst node

below the centerline. The entrance length L,= is deﬁned as the distance from the inlet

to the point where the centerline velocity reaches 99% of the parabolic centerline
velocity. These regions and their ratios are obtained for various Reynolds numbers

and are presented in Table l.

The results show that the inviscid-core region is approximately one-ﬁfth of the
entrance length and increases with a decrease in the Reynolds number below 100.
This is somewhat shorter than one-fourth the entrance length reported by Mohanty
and Asthana (1978) for ﬂow in a circular pipe. No signiﬁcant difference was noted

91

 

 

 

 

 

 

 

Table l. The inviscid-core length Li, the proﬁle-development length L d, and the
entrance length L a for various Reynolds numbers.
Re _L_i_ .113. Le Li Ld Le Li
H H H H-Re H-Re H-Re Le
10 0.15 0.35 0.50 0.0150 0.0350 0.0500 0.3000
20 0.25 0.60 0.85 0.0125 0.0300 0.0425 0.2941
50 0.50 1.60 2.10 0.0100 0.0320 0.0420 0.2381
100 0.80 3.25 4.05 0.0080 0.0325 0.0405 0.1975
200 1.55 6.65 8.20 0.0078 0.0332 0.0410 0.1899
500 4.10 16.50 20.60 0.0082 0.0330 0.0413 0.2000
1000 7.85 33.35 41.20 0.0079 0.0334 0.0412 0.1897
2000 15.55 66.70 82.25 0.0077 0.0334 0.041 1 0.1892
5000 39.10 167.40 206.50 0.0078 0.0335 0.0413 0.1893

 

 

 

 

when the entrance length was compared with that of other researchers as well as the

well known empirical relationship: L,/H = 0.04 Re. Moreover, when L e is calculated

using different inlet conditions, the entrance length was found to be marginally

inﬂuenced by the inlet velocity distribution. A summary of L .3 values for a number of

channel studies is shown in Table 2.

The streamwise velocity proﬁle undergoes a change from an assumed inlet

proﬁle (U = 1.0) to that of a fully-developed, parabolic proﬁle at a location

92

Table 2. Comparison of the ratio L,/(H-Re) as reported in different channel
studies for Reynolds number = 200.

 

Analytical Schlichting (1934) 0.0400
Experimental Feliss et al. (1977) 0.0400
Numerical Bodoia and Osterle (1961) 0.0440
Roidt and Cess (1962) 0.0454

Hwang and Fan (1963) 0.0422

Gillis and Brandt (1964) 0.0456

Morihara and Cheng (197 3) 0.0452

Abdulla (1987) 0.0442

Present study 0.0410

with boundary-layer equations 0.0405

with cascade inlet condition 0.0415

with experimental inlet proﬁle 0.0417

 

 

 

 

downstream with centerline velocity uCL = 1.5. For a uniform inlet velocity, the

streamwise velocity proﬁles for various Reynolds numbers are shown in Fig. 8.
Several of these proﬁles are enlarged in Fig. 9. They demonstrate that for x-locations
in the inviscid-core region, the velocity proﬁles include a minimum on the axis and
symmetrically located maxima on either side of the centerline. The maximum velocity

in these proﬁles is at most 4.9% higher than that at the centerline.

This is in contradiction to reported results (Gillis and Brandt, 1964; Brandt and
Gillis, 1966; Morihara and Cheng, 1973) in which these local maxima are much more
pronounced over larger downstream distances and at low Reynolds numbers; the
overshoot diminishes for Reynolds numbers > 200. The consistency in the velocity

93
proﬁles over the range of Reynolds numbers considered here suggests that the
overshoot is a result of the numerical approximations. By solving the full Navier-
Stokes equations, it is found that these local maxima are essentially nonexistent.
This point is also evident in the work by Wang and Longwell (1964) who used a
ﬁnite-difference formulation that is less accurate in the representation of the vorticity
boundary condition than that used in this study. Their boundary condition (Thom and
Apelt, 1961) is only ﬁrst-order accurate thereby overestimating the boundary
vorticity that leads to the overshoot; our method retains the second-order terms in
the Taylor series expansion used to approximate the boundary vorticity. In addition,
it can be shown that in the case of Li and Ludford (1980), the use of boundary-layer
equations near the entrance leads to an incorrect representation of the initial

vorticity. These results are supported by the numerical ﬁndings of Abdulla (1987).

The effects of several factors on the develoth of the streamwise velocity
proﬁles are investigated. The proﬁles obtained from the use of boundary-layer
equations are shown in Fig. 10 for Re = 200. The proﬁles due to cascade inlet
condition for Re = 200 are displayed in Fig. 11. An experimental proﬁle for the inlet
velocity, obtained by ﬁtting a polynomial near the wall to a uniform central section
using measurements from a hot-wire anemometer (Backus, 1989), was also used.
The streamwise velocity distributions along the channel for this experimental proﬁle
at the inlet are shown in Fig. 12 for Re = 200. The difference due to the above three
ﬂow conditions, as compared to Fig. 8e, is clearly insigniﬁcant.

The streamwise velocity as a function of the location across the channel is
shown in Fig. 13 for Re = 200 and 2000. Also, the streamwise velocity at the
centerline is compared in Fig. 14 for Re = 20, 50, 100, 200, and 2000. The shape of the

94
centerline velocity is identical to that of other studies where, with the increase of the
Reynolds number, higher velocities are reached at shorter distances from the inlet.
Proﬁles of normal velocity components along the channel are shown in Fig. 15 for Re
= 200, 1000, and 2000. The changes in v are most apparent at distances near the
inlet. Although not completely fully-developed, the ﬂow in the proﬁle-development

region becomes parallel at a short distance after the inviscid—core region.

Proﬁles of the vorticity along the channel for Re = 200, 1000, and 2000 are
shown in Fig. 16. Figure 17 shows the values of the vorticity at the upper boundary
forRe-20,200,and2000. Thelasttwoﬁguresindicatethatthevalueofboundary
vorticity levels out to reach the fully-developed value (with magnitude of 6.0) as the
ﬂow progresses downstream; the vorticity on the lower boundary is equal but
opposite in sign. The changes, as functions of x/Re, are more rapid for higher
Reynolds numbers. A comparison of boundary vorticity using different formulas is
presented in Fig. 18. It further supports, as indicated in Section 6.1, the hypothesis
that the second-order Woods’ condition underestimates the boundary vorticity in
contrast to the ﬁrst-order Thom’s condition.

Vorticity contours and streamlines of the ﬂow are displayed in Figs. 19 and 20,
respectively. It is clear that most changes take place in the immediate vicinity of the
inlet section. The normalized pressure gradients along the channel at the wall and
centerline are shown in Figs. 21, 22, and 23, for different Reynolds numbers. Both

pressure gradients converge to the constant value (- Re/ 12) of a fully-developed ﬂow.

The details of the computational procedure for determining the proper time steps

and mesh sizes must be given some special attention. Concerning the procedure for

95
determining the mesh sizes, there are two problems involved. One is to determine
whether the size of the computational grid is large enough to impose the downstream
boundary condition and match the required computations with the available computer
capabilities. The other is to determine whether the chosen mesh sizes are small
enough for accurate results. These two problems are obviously interrelated.
Conversely, the size of the time step is important in determining the computation time

and the numerical stability of the solution.

To check whether the solution is independent of the mesh size, mesh sizes in
the range of 0.02 to 0.1 are used at selected Reynolds numbers. The effect of the
mesh size on the accuracy of the results is found to be minimal; for example, when
compared to a grid size of 0.025, our standard grid (0.05) produces a maximum
difference of 4.6% in the values of the streamwise velocity. The variation of the
streamwise velocity with the mesh size at selected locations is shown in Fig. 24 for
constant At and h/ At. The relationship between u and h2 is nearly linear, therefore,

confirming the second-order accuracy of the solution in the space variables.

7.2 Channel Flow with a Finite Step

Channel ﬂow in the presence of a ﬁxed ﬁnite step on the lower boundary is
identiﬁed with separation. Two distinct separated regions with pronounced eddies of
recirculating ﬂuid are developed on both sides of the step. The size of these
separated regions depends on the Reynolds number, the height (B), and length (D) of
the step; the separated region upstream of the step is considerably smaller in size
thanthataftofthestep. AchannelwithauniformgridandspacingofsizeOﬂS is
used in the numerical solutions for Reynolds numbers 20 and 200. A ﬁnite step of

variable height and length is placed at different locations Xl along the channel, as

96

shown in Fig. 5. A summary of the cases considered is presented in Table 3.

The effect of the length of the step on the separation is demonstrated in Fig. 25
for Re = 20. Steps, located at X1 = 1.0, with length D = 0.0, 0.05, 0.1, 0.25, and 0.3

andheightB=0.25arestudied. ThestepoflengthD=0.0isaspecialcasewhere

the step is represented by the grid line between two adjacent meshes.

 

 

Table 3. Length of the aft separated region L2 for the channel ﬂow with a ﬁnite
step.
Re Xl B D L2
20 1.00 0.20 0.00 0.368
1.00 0. 25 0. 00 0.441
1.00 0. 25 0. 05 0.468
1.00 0. 25 0.10 0.456
1.00 0. 25 0. 25 0.431
0. 85 0.25 0. 30 0.412
1.00 0.10 0. 05 0.190
1.00 0.20 0.05 0.386
1.00 0.30 0.05 0.552
1.00 0.40 0.05 0.663
1.00 0.50 0.05 0.821
200 1.00 0.20 0.05 1.794
1.00 0.25 0.05 2.175
2.00 0.20 0.05 1.733
2.00 0.20 0.25 1.506
2.00 0.20 1.00 1.116
2.00 0.25 0. 05 2. 309
8.00 0.20 0.05 1.702

 

 

 

For the cases

 

97
considered,thelengthLloftheseparatedregioninfrontofthestepappearstobe
uninﬂuenced by an increase in D while the length L2 of the separated region aft of the
step is slightly decreased with an increase in D. Conversely, the role of the height of
the step in determining the size of the separated region is shown in Fig. 26 for Re =
20. The step’s height B varies between 0.1 and 0.5 with a constant length D = 0.05.
Increasing the height of the step causes the size of the aft separated region to
increase, as expected. In all of the above cases, the length of the front separated

region is observed to be somewhat ﬁxed with value L1 ~ 0.10.

TheseparationduetoastepofheightB=0.2andRe=200isshowninFig. 27
for steps of length D = 0.05, 0.25, and 1.0 located at X1 = 2.0. The separated region
becomes smaller with the increase in D. For Re = 200, Fig. 28 shows an increase in
thesizeoftheseparatedregionduetoanincreaseintheheightofthesteme=
0.25; compared to Fig. 27a, the length of the separated region in the aft of the step is
nearly doubled.

Figure 29 shows the separation due to a step placed further upstream (at X1 =
1.0) for Re = 200, D = 0.05, and B = 0.2 and 0.25, respectively. The aft separated
region is again larger for B = 0.25. Evidently, the size of the separated region is
unaffected by the streamwise location of the step. This result is also confirmed by

Fig. 30 where the step is placed at X1 = 8.0 where the ﬂow is fully-developed; the

size of the separated region remains the same regardless of the location of the step.

Velocity vectors for the ﬂow in the presence of a step with B = 0.2, D = 0.05,
and Re a 20 are shown in Fig. 31. The velocity gradients are high in the vicinity of the
step as exhibited by the direction of the vectors. Streamwise velocity proﬁles for

98
several combinations of height B and length D are shown in Fig. 32 for Re = 20 and
X1 = 1.0, and in Figs. 33 and 34 for Re = 200 and X1 = 2.0 and 1.0, respectively. The

above ﬁgures reveal that the presence of a step has an upstream effect on the shape
of the velocity proﬁle and that the velocity distribution becomes parabolic at a

sufﬁcient distance downstream.

Inaddition,acomparisonofthecasesstudiedshowsthattheincreaseinthe
Reynolds number is accompanied by an increase in the length of the separated region
aft of the step. The size of the separated region upstream of the step, however,
seems to be relatively independent of the Reynolds number. The vorticity contours,
shown in Fig. 35, emphasize that the changes in the ﬂow ﬁeld are most pronounced in
the vicinity around the step.

7.3 Channel Flow with a Moving Boundary

For this flow, a section of the lower boundary is allowed to move with a
constant velocity in a direction the same as or opposite to the entry ﬂow. The channel
studied here has a grid with a uniform spacing of size 0.05. The vorticity C, stream
function \y, and streamwise velocity u, are calculated for ﬂows with boundary velocity

UB ranging from -10.0 to 5.0 and Reynolds numbers between 20 and 1000. The

moving boundary is deﬁned to exist between x = XI and x = X2 (as shown in Fig. 1).

Results demonstrated by plots of streamlines, velocity profiles, and vorticity
contours, show the formation of a separated region when the boundary section is
moving in the upstream direction. In general, the size of the separated region
depends on the Reynolds number and boundary velocity. Streamlines of the ﬂow for

99
Re a 20 with the lower boundary section (between X1 = 1.0 and X2 = 2.0) moving
opposite to the inlet velocity at UB = -10, -5, -2, and -1, are shown in Fig. 36. The
separated region becomes larger with an increase in the magnitude of boundary
velocity. Velocity vectors for the case Us = -10 are shown in Fig. 37 indicating the
recirculating pattern of the ﬂow inside the large separated region. Figure 38 shows
the flow streamlines when the boundary section is moving in the direction of the inlet
velocity at UB = 1, 2, and 5 for Re = 20. Observe that no separation occurs since the

boundary is moving in the direction of the ﬂow. The ﬂuid in the near vicinity of the
moving boundary section is accelerated due to the motion of the boundary.

The above observations are revealed in Figs. 39 and 40 for Re = 200.
Streamlines of the ﬂow with boundary velocity UB = -2.0, -1.0, -0.5 and 1.0 are
displayed. The inﬂuence of the Reynolds number on the size of the separated region
is demonstrated by comparing Fig. 36 and Fig. 39 for U}, = -1.0 and -2.0. It is clear
thatthesizeoftheseparatedregiondecreasesasReincreases. Thenumberof
iterations needed to attain convergence to a steady-state solution increases with an

increase in Re but decreases with a decrease in the magnitude of U3.

The size of the separated region is also related to the length and location of the
moving boundary section. As shown in Fig. 41, the separated region extends along
the length of the moving boundary with no apparent increase in its height due to the
increase in the length of the moving section of the boundary. Figures 42 and 43 show
that placing the moving section of the boundary further downstream (between X1 = 4
and X2 = 5) produces no signiﬁcant change in the size of the separated region; in this

case, however, a larger number of iterations is required due to the extended

1w
computational domain. Further examples regarding the effects of the length and
position of the moving boundary section are shown in Figs. 44 and 45, for Re = 200
and Us 2 -1.0, -0.5, and 1.0. In addition, placing the moving boundary section in the
fully-developed region of the flow showed no noticeable difference in the streamlines

of the ﬂow or the shape of the separated region, as shown in Fig. 46.

For Re a 1000, streamlines are shown in ﬁg. 47 when the lower boundary
section, positioned between X1 = 10 and X2 = 14, is moving opposite to the inlet
velocity at UB = -0.5; at Re = 1000 the separated region is reduced in size since it is
apparently suppressed by the relatively high velocity incoming ﬂow. Flow
streamlines for the Couette-Poiseuille ﬂow, when the entire lower boundary is
moving in the direction ofthe ﬂow at UB =1.0, for Re = 500 are shown in Fig. 48.

Development of streamwise velocity proﬁles along the channel are also shown

in Figs. 49 and 50 for UB = -1.0 , -0.5, and 1.0 and Re = 200, and in Fig. 51 for UB
= -0.5 and Re = 1000. The velocity at the boundary between X1 and X2 is equal to
the boundary velocity as required by the no-slip boundary condition. The parabolic
velocity proﬁle is attained at some distance beyond the inﬂuence of the moving
boundary. Velocity proﬁles for the Couette-Poiseuille ﬂow case are presented in Fig.
52 for Re = 200 and 500.

Vorticity contours with and without the presence of a separated region are
shown in Figs. 53 and 54 for Re = 20 and Re == 200, respectively. Evidently, changes
in the vorticity are concentrated in the immediate vicinity of the moving boundary
section. For the Couette-Poiseuille ﬂow case, vorticity contours are shown in Fig.

55, as well.

101

7.4 Channel Flow with an Oscillating Boundary

An extension of the previous problem is the introduction of an oscillatory motion
to the moving section of the lower boundary, which it is permitted to oscillate
sinusoidally in its own plane. The vorticity C, stream function w, and velocity
component u are calculated at each nodal point for a ﬂow with a uniform inlet velocity
proﬁle and boundary velocity U3 = A Cos (t) where A is the amplitude of oscillation of
boundary velocity. A mesh size of 0.05 is used with Reynolds numbers of 20, 200,
and 1000. Again, the oscillating boundary is deﬁned to exist between x = XI and x =

X2 (as shown in Fig. 1).

As predicted, plots of streamlines, as well as velocity proﬁles and vorticity
contours, demonstrate the existence of a separated region generated in the vicinity of
the oscillating boundary section during its upstream motion. Generally, the size of
the separated region depends on the Reynolds number, Strouhal number, the
amplitude of boundary velocity, and time; its size is a maximum at t = (2n+1)rt where
n is an integer representing the number of cycles of oscillation before the solution
converges. Streamlines for the flow with Re = 20, St = 1.0, and A = 1.0 are shown in
Fig. 56 at time t = 2m: and (2n+1)lt. For Re = 200, and different St and A,
streamline patterns are shown in Figs 57 and 58 at t = 2n1c and (2n+1)x. Note that
no separated region appears at the instant t = 2nrt since the boundary is moving with
the ﬂow, as in Fig. 56a, and that the formation of a separated region is clear at t =
(2n+ 1):: when the boundary is moving at maximum velocity in the direction opposite
to the flow, as in Fig. 56b. The inﬂuence of the oscillating boundary on the

progression of the ﬂow is further illustrated in Fig. 59 where a full cycle of oscillation

we
is shown; for this example, Re = 200, St = 1.0, and A = 1.0. The size of the separated
region increases with time but starts to decrease for t > (2n+l)1t. At t = (2n+l)rt,
the height of the separated region is equal to 0.098. As the velocity of the boundary

becomes > 0, the separated region ceases to exist.

The inﬂuence of the Strouhal number on the size of the separated region is
demonstrated by Figs. 60-62 for St = 0.5, 1.0, and 2.0, Re = 200, and various A and
length of the oscillating boundary. In all three cases, a separated region exists since
themotionoftheoscillatingboundaryisoppositetotheﬂowatthatpointintime(t=
(2n+l)rt). Clearly, the size of the separated region is inversely proportional with St
which is directly related to the ﬁequency of oscillation. The effect of amplitude of
oscillation A of boundary velocity is shown in Figs. 63—66 for A between 0.5 and 2.0
with different Re, St, and length of the oscillating boundary. Increasing the amplitude

of oscillation causes the separated region to increase in size.

Figure 67 displays stream function contours of the ﬂow at t = (2n+-1): for Re =
1000 and different St and A. The separated region barely exists since it is apparently
swept away by the relatively high velocity incoming ﬂow. The effect of the Reynolds
number on the size of the separated region can be established from the previous
ﬁgures, in particular, Figs. 56b, 60b, and 67a, for Re = 20, 200, and 1000, respectively,
where for all three ﬂows St = 1.0, A = 1.0, X1 = 1.0 and X2 = 2.0. It is evident that,
for ﬁxed St and A, the size of the separated region decreases as Re increases.
Comparing Figs. 63a and 60c for St = 2.0, A = 1.0, X1 = 1.0, X2 = 2.0, and Re = 20 and
200, respectively, also leads to the same conclusion. The number of iteration cycles it

necessary to obtain convergence increases with an increase in Re and St.

103

In addition, the size of the separated region may depend on the location and
length of the oscillating boundary section. Placing the oscillating section of the
boundary further downstream, as shown in Fig. 68, has no signiﬁcant effect on the
size of the separated region (as compared to Fig. 60b); in this case, however, a larger
number of iteration cycles is required due to the extended computational domain. On
the other hand, Fig. 69 shows that the separated region extends along the length of
the oscillating boundary with no apparent increase in its height due to the increase in
the length of the oscillating section of the boundary. The streamlines of the ﬂows
shown in Fig. 70 combine the effects of both location and length. The oscillating
boundary section is placed between X1 = 8.0 and X2 = 10.0; this is equivalent to a
parabolic inlet proﬁle. When compared to Figs. 61b and 69, for the ﬂow with Re =
200, St = 1.0 and A = 1.0, Fig. 70b illustrates that placing the oscillating boundary
section further downstream has no inﬂuence on the separated region. Moreover, from
Figs. 65a and 70c it is clear that the above is true for the case Re = 200, St = 2.0, and
A = 0.5.

Developments of streamwise velocity proﬁles along the channel, for several ﬂow

conditions, are also shown in Fig. 71 at t = 2th and in Fig. 72 at t = (2n+l)rt.
Variations in the velocity proﬁles due to the presence of the oscillating boundary
appear to be local. The velocity at the boundary between X1 and X2 is equal to the

boundary velocity as required by the no-slip boundary condition. A parabolic velocity

proﬁle prevails at the downstream end.

Vorticity contours for the flow with St = 1.0, A = 1.0 XI: 1.0, and X2 = 2.0 are

shown in Figs. 73 and 74 for Re = 20 and 200, respectively; these contours are at t =

104
(2n+1)rt with UB = -1.0 and t = 2n: with UB = 1.0. Vorticity contours are also
shown in Fig. 75 for the ﬂow with the oscillating boundary placed in the fully-
developed region. The above ﬁgures demonstrate that changes in the vorticity are
substantial only in the near vicinity of the oscillating section of the boundary. It is
interesting to notice that the vorticity at the centerline is almost zero, which suggests
that the oscillating boundary has a negligible effect on the ﬂow in the upper half of the
channel as long as the distance between the lower and upper boundaries is large

relative to the height of the separated region.

7.5 Accelerating Parameters

For the present numerical solution, the rate of convergence of the difference
equations is signiﬁcantly enhanced through the use of accelerating parameters.
Optimum values of the time step (At) for the vorticity equation, over-relaxation factor
(F8) for Poisson’s equation, and the optimum values of the weighting factors (a and

y) for the vorticity and stream function, respectively, are presented below.

In order to determine the optimum values of the accelerating parameters,
numerical experimentation with a large number of combinations of these parameters
is attempted for a channel with a uniform mesh of size 0.05. First, optimum values for
At are obtained using the optimum FS given by Eq. (61). The optimum values of F8
are then reﬁned by trying F8 in the range 1.0 to 1.9; FS = 1.0 corresponds to the
Gauss-Seidel iterative scheme. Weighting factors for the vorticity and stream
function are then applied, as prescribed by Eqs. (80) and (81), respectively. Values

of a and y are assumed to vary between 0.0 and 0.9.

H5
Optimum values of the accelerating parameters as functions of the Reynolds
number, based on the numerical solution over one half the height of the channel, are
given in Table 4 for Re = 20, 50, 100, 200, 500, 1000, and 2000. The number of outer
iterations and computing time for the accelerated solution along with the number of

iterations for the Gauss-Seidel solution, considered to be the unaccelerated case (F S

= 1.0, o = 0.0, and y = 0.0), are also given in Table 4. The same information is

Table 4. Optimum accelerating parameters, CPU time, and number of outer
iterations for the channel ﬂow with a uniform mesh of size 0.05 at various
Reynolds numbers; the solution is based on one half the height of the

channel.

 

 

rmmMrMFPmmurd'

Re At FS 0 7 CPU time am Wed
20 0.03 1.70 0.00 0.05 38.5 58 l 15
50 0.09 1.68 0.00 0.05 44.9 54 147
100 0.16 1.67 0.00 0.05 1:267 55 160
200 0.33 1.65 0.00 0.05 3:084 57 170
500 0.63 1.64 0.00 0.10 9:415 74 202
1000 0.67 1.64 0.00 0.15 31:59.9 122 228
2000 0.67 1.64 0.00 0.25 2:10:02.1 247 "'

 

 

 

 

l.(TUﬁmekgwmmhmmee

1L DwanwmmnammmmkuwmwdhwqﬁmmnmakmmmpmmMms

3. TheunacceleratedcaserefasmﬂleGauss-SeidelitemdcnwidropﬁmumALF821.0.6
a 0.0, andyc 0.0.

t ‘hﬂhmxmadnamﬁaubwnmumwnpmemmnmmmmm

1%

Table 5. Optimum accelerating parameters, CPU time, and number of outer
iterations for the channel ﬂow with a uniform mesh of size 0.05 at various
Reynolds numbers; the solution is based on the full height of the channel.

 

 

 

 

 

 

Number of Number of
Re At FS 0 7 CPU time Iterations; Iterations;
accelerated umccelerated
20 0.01 1.78 0.00 0.05 2:05.9 113 145
200 0.09 1.80 0.00 0.05 9:05.8 92 214
1000 0.24 1.80 0.00 0.10 2:18:29.0 266 406
2000 0.24 1.80 0.00 0.20 11:51:08.0 687 913
Notes: 1. CPUtirneisgiveninhrzminzsec.
2. Theaccelaatedeasehlcltxleﬂlemeofﬂleopﬁmumacceleraﬁngpamneters.
3. ThemacceleratedeasenfaswdreGamSeiddiwraﬁmwithopﬁmumALFS=lno

- 0.0, and y :- 0.0.

presented in Table 5 for Re 2 20, 200, 1000, and 2000 for the solution based on the full
height of the channel.

It is found that the optimum value of At increases with an increase of Re until it
achieves a constant value at Re > 500. For the solution of one-half the height of the
channel, this constant value is equal to 0.67, whereas for the full height of the channel
it is equal to 0.24. For all Reynolds numbers, the number of iterations required for

convergence initially decreases as At increases until the optimum value of At is

attained. Increasing At above the optimum value causes the number of iterations, and

hence the time of computation, to increase considerably. The solution, however,

107

becomes divergent as At continues to increase beyond a certain critical value.

The numerical results also demonstrate that the optimum value of FS is related
to the number of nodal points MxN of the computational domain, which is directly
proportional to the Reynolds number and mesh size. The computing time decreases
as FS increases from 1.0 until a minimum number of outer iterations is reached at the
optimum value of FS. Continued increases in the value of FS result in an increase in
the number of iterations. It is clear that the optimum value of FS decreases as the
Reynolds number increases up to Re = 500. However, for higher Reynolds numbers
the optimum value of FS approaches constant values equal to 1.64 and 1.8,
respectively, for solutions based on one half and the full height of the channel. It is
also noted that the results for Re < 500 converge most rapidly when the value of FS
is slightly higher than that given by the analytical expression of Eq. (61); this value is

approached asymptotically at higher Reynolds numbers.

The weighting factors 0 and 7 also have a signiﬁcant effect on the convergence
of the numerical solution. Again, they depend on the total number of nodal points of
the computational domain. The computing time decreases as 0‘ decreases to zero for
all Reynolds numbers. This result can be explained as being the consequence of the
reﬁnement of the optimum value of At; an optimum At can produce an optimum
convergence rate of the vorticity equation, thus implying that the rate of convergence
of successive iterations needs not be slowed down. Unlike a, a minimum solution
speed is obtained for values of 7 above zero. The optimum value of 7 has a constant
value of 0.05 for the lowest Reynolds numbers and increases steadily for Re > 200 to

reach a value as high as 0.25 for Re = 2000. Increasing the value of 7 above the

108

optimum value increased the number of iterations.

With the use of the optimum values of the accelerating parameters, the required
number of outer iterations is reduced approximately by a factor of 2 to 3 as compared
to that required by the unaccelerated case. This is evident from the number of
iterations displayed in Tables 4 and 5. Furthermore, for the range of Reynolds
numbers considered, the time of computation for this ADI/SOR scheme is better than
2.5 times less than that required by the iterative solution employed by Abdulla (1987).

Finally, the numerical results show that the optimum values of the accelerating
parameters which are used to reduce the computing time for the plane channel ﬂow,
also represent the optimum values for the channel ﬂow with a step or a moving

boundary.

7.6 Conclusiom

In closing, several conclusions are supported by the results of this numerical
study. These are summarized below.

1. For the two-dimensional, viscous, incompressible entry ﬂow, the entrance
length was found to be compatible with known experimental results, i.e., L e =

0.04 Re. The length of the inviscid-core region is approximately one-ﬁfth of the
entrance length. The inﬂuence of a cascade inlet condition or an experimental
inlet velocity proﬁle (with velocity gradients near the walls) on the development
of streamwise velocity proﬁles is insigniﬁcant. The local maxima resulting from
the overshoot in the streamwise velocity proﬁle near the walls, as predicted by

other researchers, are essentially negligible.

109
For the ﬂow in a channel constricted with a ﬁnite step, separation occurs on both
sides of the step. The size of the upstream separated region is smaller and
does not depends on the Reynolds number. The size of the aft separated region
is directly influenced by the Reynolds number and the height of the step but is
unaffected by the length of the step. Also, the size of these separated regions is
independent of the streamwise location of the step.

For the channel ﬂow with a boundary section moving at a constant velocity, the
upstream motion of the boundary clearly generates a separated region. The ﬂow
in the vicinity of the moving boundary is accelerated when the boundary is
moving in the direction of the ﬂow. The size of the separated region is a function
of the Reynolds number and boundary velocity. It decreases as Re increases

but increases with the magnitude of velocity Us of the moving boundary section.
The increase in Re and the magnitude of UB causes the number of iterations

needed to attain convergence to a steady-state solution to increase. Changes

are concentrated in the immediate vicinity of the moving boundary.

For the channel ﬂow with an oscillating boundary section, the formation of a
separated region is clear when the boundary is moving opposite to the ﬂow; no
separated region appears when the boundary is moving in the direction of the
ﬂow. The size of the separated region is inversely proportional to Re and St but
directly proportional to the amplitude of boundary oscillation. The number of
iteration cycles (i.e., time) needed to attain convergence increases with the

increase in Re and St. The effect of the streamwise location of the oscillating

boundary is insigniﬁcant.

110
The incorporation of accelerating parameters extends the domain of convergence
and improves the stability of the numerical solution, as well as minimizes the
time of computation by as much as a factor of 3 for certain Reynolds numbers.
Their optimal values are functions of the number of nodal points in the
computational domain; for each Reynolds number and mesh size there is an
optimum combination of these parameters. The optimum values of the four
accelerating parameters determined in this study for a plane channel ﬂow, are
also applicable for a channel ﬂow with a step or a moving boundary section;
these parameters should serve as a guide to reduce the computing time for other
ﬂow situations, which use this system of equations and have the same

computational domain.

The ADI/SOR scheme, utilizing optimum accelerating parameters, is
numerically stable for high Reynolds numbers with the mesh size considered
and appears to be a suitable technique for treating ﬂow problems that contain
time-dependent boundary conditions. It also provides the advantage of being
more rapidly convergent than the usual iterative schemes. Results obtained
using a ﬁner mesh produce a maximum difference of 4.9% in the values of the

streamwise velocity.

111

1.00 .
0.75:
y 0.50;
0.25-1
0.00 l , . , . T . , . , .

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0.0 0.1 0.2 0.3 0.4 0.5 0.6
x
a) Re=10.
1.m m V I
0.75-4
y 0350'1
0.25-
anol «
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2
x
b) Re=20.
1.00‘ N
0.75-
y 0.50“
0.25~
0.00 /
0.0 1.0 2.0 ISIVO 4.0 50. 60. W70 8.0 10.0 11.0 12.0
x
c) Re=200.

Figure 8. Development of streamwise velocity proﬁles through the entry region of
a two-dimensional channel at various Reynolds numbers.

112

1.00 m '1
0.75-:
y0.50j
0.25-j
.,.,.1r1r,.T.r1j.,

 

 

 

 

 

 

 

 

 

 

0.00
0.0 0.3 0.6 0.9 1.2 1.5 1.8 2.1 2.4 2.7 3.0
x
d) Re=100.
1.00 '\ 1 I r
0.75-
y 0.50-J
0.25«
0.00 / , eﬁ . , .
0.0 0.5 .5 W25? 3? 3.5 4.04.5 5.0 5.5 6.0
x
e) Re=200.

 

r

1.00. N V l
075*
y 0.50:
0.25;
0.00‘ . r i

0.0 2.5 #570 7.5 10.0 12.5 15.0 17.5 20.0 22.5 25.0 27.5 30.0

 

 

 

 

X

f) Re=500.

Figure 8. Development of streamwise velocity proﬁles through the entry region of
a two-dimensional channel at various Reynolds numbers.

113

1.00 .\ 1.1
0.75;
y 0.50;
0.25-1
'TﬁT' 'rrrﬁflrt'lrﬁ1

 

 

 

0.00 ﬁ/ 1 I
0.0 5.0 10.0 15.0 20.0 25.0 30.0 35.0 40.0 45.0 50.0 55.0 50.0

 

 

 

 

 

 

 

 

x
g) Re=1000.
1.00 .\ r1 e
0.75:
)1 0.50“:
0.25-
0.00. /
0.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0100.0110.0120.0
x
h) Re=2000.
1.00 ~\
0.75:
y 0.50“:
0.25:
0.001 /

 

 

 

 

0.0 25.0 50.0 75.0100012501500175020002250250027503000

X

i) Re = 5000.

Figure 8. Development of streamwise velocity proﬁles though the entry region of
a two-dimensional channel at various Reynolds numbers.

114

 

 

 
 
  
   
 
   
 

.G-OR.=20,x=O.1
.MRO=200,x=U.1
HR.-2W,x-U.5
HR.-200.x-1.0
G-EIROSZCKKJ,x=5.D

 

 

T I V T

'H ,, .. .fﬁﬁt ﬁr i
0.00 0.25 0.50 0.75 1.00 1 .25 1 .50

 

 

Figure 9. Streamwise velocity proﬁles at locations near the inlet section for
various Reynolds numbers.

F1

115

IoOOJ \
0.75e
y 0.50‘
0.25..
. M
T'l‘l'l'l‘l'l'
4.0 5.0 6.0 7.0 8.0

0.00 v v
0.0 1.0 2.0 3.0

 

>>

1
I Y I ' ﬁ

l
r
9.0 10.0 11.0 12.0

 

 

 

 

X

Figure 10. Development of streamwise velocity proﬁles through the entry region of
the channel using the boundary-layer equations at Re = 200.

1 .00 v
0.75;

y 0.50-1
0.25:

CHOW—F":

-1.0 0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0

 

 

 

 

 

 

X

Figure 11. Development of streamwise velocity proﬁles through the entry region of
the channel using a cascade inlet condition at Re = 200.

171'11I'
4050

 

 

1.00
0.75j
y 0.501
0.25;
0.001 . I .

r I j’ I ' I t T ' I I '
0.0 1.0 2.0 3.0 . . 6.0 7.0 8.0 9.0 10.0 11.0 12.0

 

 

X

Figure 12. Development of streamwise velocity proﬁles through the entry region of
the channel using an actual experimental inlet velocity proﬁle at Re = 200.

116

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1.5tttrIrtvII1v1_1_ITIII:Jxﬁ
1.2-
0.9—
3 d
3 1
0.6A
0.1%
00.. ,--..,.f-ﬁ. fﬁﬁ--.
0 00 O 01 O 02 0.05 0.04 0 05
x/Re
a) Re=200.
1’5 rTfﬁl"'ffﬁ'*'1 w
1.2“
d
0.9—
3 .
3 1
0.6-
0.3-4
0.0 ....,....T....,...-1....
0.00 0.01 0.02 0.03 0.04 0.05
x/Re
b) Re=2000.

Figure 13. Streamwise velocity at different y-locations across the channel for Re =
200 and 2000.

117

 

 

 

 

 

x/Re

Figure 14. Centerline velocity along the channel for various Reynolds numbers.

118

 

 

 

 

 

 

 

 

 

 

 

1.00
0.75-1
y 0.50-1
0.25:
0.00. v I T
0 0.1

0.

 

T r ' l ' I I
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 .0

a) Re=200.

 

1.00
035‘ l ‘
y 0.5m

0.254

 

 

 

 

 

 

 

 

 

 

 

 

0.00 I ' l ' l
0.0 1 .0 2.0 3.0

 

5.0 6.0 7.0 8.0 9.0 1 0.0

*4
o

b) Re=1000.

1.00 ’ 7’

0.75; ‘
y 0.503

0.25: l L |

0.000 r16.o'26.o'3or.o' r fr '

0. 40.0 50.0 60.0 70.0 80.0 90.0 100.0

 

 

 

 

 

 

 

 

 

 

 

 

 

X

c) Re = 2000.

Figure 15. Development of normal velocity proﬁles through the entry region of a
two-dimensional channel at various Reynolds numbers.

119

1°00 ' I I l ' l ' l I r I I 7 I ' I I l
.1

0.75-

y 0.50“

I .

0.25-j

D-ﬁi ' I r I I F V I 1 I 1 I I I r I I I I

0.1 0.2 0.3

 

 

 

 

 

 

f
0.0 0.4 0.5 0.6 0.7 0.8 0.9 1 .0
X
a) Re = 200.
1.m T I Y I I I T I 1 1 I I I I I I I I I
0.75-
y 0.50~
0.215q
0.“) T 1 1 I f r T I I I I I I I

 

 

 

r I U I T I
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 1 0.0

 

x
b) Re=1000.
1-00'T'I'r'rTIv1rII. r
0.75-j
y 0.501
0.25:
0.001

 

 

 

I T . I I I I I ' f I I I I I r ' I
0.0 1 0.0 20. 0 30.0 40 .0 50.0 60.0 70.0 80.0 90.0 100.0
X

c) Re = 2000.

Figure 16. Development of vorticity proﬁles through the entry region of a two-
dimensional channel at various Reynolds numbers.

120

 

 

 

 

 

 

 

 

 

 

35.0
0 Re - 20
A Re - 200
30.0 o D R. = 2000
25.0
20.0
15.0
10.0
5.0~
0.0 T W T T—I v I I v I I I t r I I I T wiI T T T 1
0.00 0.01 0.02 0.03 0.04 0.05
x/Re
Figure 17. Vorticity at the upper boundary along the channel at various Reynolds
numbers.
35.0
are Woods'
“ H Jenson!
300.1 H Thom'l
1
25.0—
20.0“
c, .
15.04
10.0-
5.0-
0.0 I j T rﬁ I r I I W T r f I I I I I T I I V I
0.00 0.00 0.00 0.01 0.01 0.01
x/Re

Figure 18. Vorticity at the upper boundary along the channel calculated from
different formulas for Re = 200.

121

 

 

 

 

 

 

 

 

 

 

I IjleIIIIIIIIIIFIITIIIIIﬁITIIITiI
-
n—q
-
‘
-q
ﬂ
q

2. f... 3'

ﬂ
c—I
q
q
—
q
—

1 llllllllllllllllLllLlllllllLJJlll

a)Re=200.

 

 

 

a) Re = 2000.

Figure 19. Vorticity contours of the flow in the entry region of a two-dimensional
channel at Re = 200.

 

 

 

 

 

 

 

 

 

Figure 20. Streamlines of the ﬂow in the entry region of a two-dimensional channel
at Re = 200.

122

 

 

 

 

 

 

 

 

35.0
4 o O R. - 20
A Re - 200
30°C? c: R: = 2000
=3 25.0:
A
é 20.o-« .
\ J
3 15.0-
53 10.0-
\ J .
cg 5.0-r
T 1. faint-IIIIIIII-I'I-.--qu.nu.luau-ul-ut
0.0 canon. o n 0
—5.0
—10-0 I ‘% I T I I I I I T I I T I T I I I T I T
0.00 0.01 0.02 0.03 0.04
x / Re
Figure 21. Normalized pressure gradients at the wall along the channel for various
Reynolds numbers.
12.5
3 0 R0 - 20
. A Re = 200
a U R. = 2000
8 10.0-‘
3 3
5 1
3? 7.5: I
v I1
\ 4
o. 4 .
3 5.0-j
I"\ '1 g
(1:. 1% A
E 2.5: ”ﬁg“: - .0. ' o o . - . .
g : "lIIII-Iu-ollI:-I.Ial.-Po-lu-Ouloll
.1
I 0.0:
-2.5 1} I. T T r I r 1 I I r ﬁ I I I I
0.00 0.01 0.02 0.03 0.04-
x/Re
Figure 22. Normalized pressure gradients at the centerline along the channel for

various Reynolds numbers.

123

 

 

 

 

30.0 c
I 0 Vol
I A Centerline
x 20.0—
U .
\ .
a d
U
.l
RT 10.0—g
\ j ‘
Q) a
9:2 ‘ o “a
I 4 1“...."“‘0llIt.no.0.0000000000000000
0.0— .o
-10.0 .2 1 I I I I I I I I f I I I l' I I I I I
0.00 0.01 0.02 0.03 0.04

x/Re

Figure 23. Normalized pressure gradients at the wall and centerline along the
channel for Re = 200.

124

 

HOW)
8
‘0
(7|

J [1

I o x=D.2, y=0.2, h/M=I
d =
1.080 El x-o.z, y-0.2, At 0.05

IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

0.000 0.002 0.004 0.605 0.008 0.010

 

 

 

a) x = 0.2, y = 0.2.

 

1.445

1

11

—-L
A
4..
N
l
U
D

lllLllllllllJJlJllllLll
DO

1 .439

U(x.y)

2?:

OD

1 .433

o x=D.5, y=0.5, h/M=l
1.430 E’ﬁ’i'Pi' 2"9‘5’: 933.035.

0.000 0.002 0004' ' ' ' 0606 0.008 0.010
hZ

 

 

 

I IIIII'IIIIII'TIIT‘FIII

b) x = 0.5, y = 0.5.

Figure 24. Streamwise velocity as a function of the mesh size at selected locations;
Re = 200.

125

 

a) D = 0.0.

 

b) D = 0.05.

 

c) D = 0.1.

Figure 25. Flow separation due to a ﬁnite step of various width; X1 = 1.0, B = 0.25,
Re = 20.

176

 

d) D = 0.25.

 

e) D = 0.30, (x1 = 0.85).

Figure 25. Flow separation due to a ﬁnite step of various width; X1 = 1.0, B = 0.25,
Re = 20.

 

a) B =0.1.

 

b) B = 0.2.

 

c) B = 0.3.

Figure 26. Flow separation due to a finite step of various height; X1 = 1.0, D = 0.05,
Re = 20.

 

a) D = 0.05.

 

b) D = 0.25.

 

c) D = 1.0.

Figure 27. Flow separation due to a ﬁnite step of various width; X1 = 2.0, B = 0.20,
Re = 200.

130

 

Figure28. Flow separation due to a ﬁnite step; Re = 200, X1: 2.0, B = 0.25, D =
0.05.

 

a) B = 0.20.

 

b) B = 0.25.

Figure 29. Flow separation due to a ﬁnite step of various height; X1 = 1.0, D = 0.05,
Re = 200.

l3]

 

Figure 30. Effect of step location on the ﬂow separation due to a ﬁnite step; X 1 =
8.0, B = 0.20, D = 0.05, Re = 200.

 

---—---------_---_-----— -------------

 

—.....—.—

 

.—.———._._.

 

“.4“...

 

_._-_._._..

—.—_———_—__—..._._.—

.4-.- — ......................

 

 

 

Figure 31. Velocity vectors of the ﬂow in the presence of ﬂow separation due to a
ﬁnite step; X1 = 1.0, B = 0.20, D = 0.05, Re = 20.

132

 

 

 

 

1.00 m T j I
0
y 0.50;
0.25:
m:_...=J

 

 

0.0 0.1 0.2 0.3 0.4 0-5 0.6 0.7 0.8 0.9 1-0 1.1 1.2

a) B = 0.25, D = 0.0.

 

1.00

 

 

 

 

 

1 NW ' I r
0.75~
, 0.503
0.25;
0.00‘ 1“")

 

0.0 0.1 0.2 0.3 0.4 0-5 0.6 0.7 0.8 0.9 1-0 1.1 1.2

b) B = 0.25, D = 0.05.

1 .00 .w
. \
D.75~
y 0.50-1
0.25- J

0’00 I F I W
0.0 0.1 0"2 0.3 W04'05106'107T0‘Br01‘9 1.0 1.1 1.2 1.3

 

 

 

 

 

 

 

X

c) B = 0.5, D = 0.05.

Figure 32. Development of streamwise velocity proﬁles along the channel in the
presence ofa ﬁnite step; Xl =1.0, Re = 20.

133

1.00 .\
0.75:
y 0.501
0.251
‘ /

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0.00 I
0.0 1.0 20 31:0 40- 60. T70 80 10.011.012.0
x
a) D=0.05.
1.00 N
0.75~
1
y 0.50“
0.25-
0.00 ./
0.0 1.0 I20 .30 40. 60. 7'0. 8.0 10.011.012.0
x
b) D=0.25.
1.00 I\ ' l
0.75%
y 0.50“
0.25~
0.00 ./ . , . . .
0.0 1.02.0 .0 60. 7.0 «80 10.0 11.0 12.0
x
c) D=1.0.

Figure 33. Development of streamwise velocity proﬁles along the channel in the
presence of a ﬁnite step; X1 = 2.0, B = 0.2, Re = 200.

134

1 .00 N
0.75 1
4
y 0.50 r
.
0.25 ~
0.00 /

0.0 1.0 2.0 ':30 4.0 60. r{7078'10 10.011.012.0

 

 

 

 

 

Figure 34. Development of streamwise velocity proﬁles along the channel in the
presence of a ﬁnite step; X1 = 1.0, B = 0.25, D = 0.05, Re = 200.

 

IIIIIIIIIITITIIIIIIIII‘IITIIIWT II

  
 
 
  

llll

 

 

 

E
I
L
1
l1

 

 

llJl

 

ll LlJlllLlLlLL

a) B = 0.20, D = 0.05.

 

 

 

 

 

 

 

lJlJllllli/llLlllll

 

1 141 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

b) B = 0.25, D = 0.25.

Figure 35. Vorticity contours in the presence of ﬂow separation due to a f'mite step;
Xl =1.0, Re = 20.

135

 

_lll|lllllllllllllllllrITlllllIllIIIITIIIIIIIIIIIIIIIlllllll

 

 

 

a) UB=-10.0.

 

b) UB =-5.0.

 

c) UB=-2.0.

Figure 36. Streamlines of the ﬂow with the lower boundary section moving opposite
to the inlet velocity at UB;X1=1.0, X2 = 2.0, Re = 20.

 

 

d) UB=-l.0.

Figure 36. Streamlines of the ﬂow with the lower boundary section moving opposite
to the inlet velocity at UB; X1 = 1.0, X2 = 2.0, Re = 20.

 

‘.----—_-—-—~—~‘--~~na

 

 

 

 

Figure 37. Velocity vectors of the ﬂow with the lower boundary section moving
opposite to the inlet velocity at UB = -10.0; X1 = 1.0, X2 = 2.0, Re = 20.

137

 

 

a) Un=1-0-

 

b) 0,, = 2.0.

 

TI—TIIIIITllllllllllIlllllllllllllllllllllllIllllllllilllill

xx. 960

    

 

 

 

c) UB = 5.0.

Figure38. Streamlines of the flow with the lower boundary section moving in the
direction ofthe inlet velocity at UB;X1=1.0, X2 = 2.0, Re = 20.

 

a) UB=-2.0.

 

 

 

c) UB=-0.5.

Figure 39. Streamlines of the ﬂow with the lower boundary section moving opposite
to the inlet velocity atUB;X1=1.0, X2 = 2.0, Re = 200.

139

 

 

 

 

 

 

 

 

 

Figure 40. Streamlines of the ﬂow with the lower boundary section moving in the
direction of the inlet velocity at UB = 1.0; X1 = 1.0, X2 = 2.0, Re = 200.

 

 

 

 

 

 

 

 

 

 

Figure 41. Streamlines of the ﬂow with the lower boundary section moving opposite
to the inlet velocity at UB = -l.0; X1 =1.0, X2 = 5.0, Re = 200.

140

 

 

 

 

 

 

b) UB = -1.0.

Figure 42. Streamlines of the flow with the lower boundary section moving opposite
to the inlet velocity at UB; X1 = 4.0, X2 = 5.0, Re = 200.

141

 

 

 

a) U3 3 -2.0.

 

 

 

b) UB=-l.0.

Figure 43. Streamlines of the flow with the lower boundary section moving opposite
to the inlet velocity at UB; X1 = 4.0, X2 = 5.0, Re = 200. (A close up of

Fig. 42.)

142

 

 

 

a) U3 3 -l.0.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

c) UB=1.0.

Figure 44. Streamlines of the ﬂow with the lower boundary section moving at UB;
X1 = 2.0, X2 = 4.0, Re = 200.

143

 

 

 

 

 

 

 

 

 

a) U3 = -l.o.

 

 

 

 

 

 

 

 

b) UB 2-0.5.

 

 

 

 

 

 

 

c) UB = 1.0.

Figure 45. Streamlines of the ﬂow with the lower boundary section moving at UB;
X1: 2.0, X2 = 6.0, Re = 200.

144

 

a) U3=-1.0.

 

b) UB=-0.5.

 

 

c) UB = 1.0.

Figure 46. Streamlines of the ﬂow with the lower boundary section moving at UB;
X1 = 8.0, X2 =10.0, Re = 200.

145

 

 

 

 

 

 

 

 

Figure 47. Streamlines of the flow with the lower boundary section moving opposite
to the inlet velocity at UB = -0.5; X1 = 10.0, X2 = 14.0, Re = 1000.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

.960 .960 .960———:
~ :1
ha.“- . 729- .72! _.729—_-:
H“ -‘
Rh; :
H __
M J“ 480 Jae—g
:-
g. —. 240 .240 .240—-—-;l
1-

 

 

 

 

 

 

Figure 48. Streamlines of the ﬂow with the entire lower boundary moving in the
direction of the inlet velocity (Couette-Poiseuille ﬂow) at UB = 1.0; Re =

500.

146

 

1.00 \ . , .
0.75;

y 0.50-:
0.25:
0.00‘ ./

 

 

 

 

0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 12.0

X

 

 

 

 

 

 

a) Una-1.0.
1.00 .\ . r 4
0.75%
y 0.50:
0.25;
0.001 v r
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 12.0
x
b) Un=-0.5.
1.00 .\ . r .
0.75;
y 0.50:
0.253
0.00. / 4

 

 

 

 

0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 12.0

X

C) U3 3 1.0.

Figure 49. Development of streamwise velocity proﬁles along the channel with the
lower boundary section moving at UB; X1 = 2.0, X2 = 4.0, Re = 200.

147

1&1 .\ . , .
0.751
y 0.50-1
0.25:
0"" / . - r , . . r 1 . .
0 3.0 4.0 5.0 6.0

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

.w I l , 1
0.0 1. 2.0 7.0 8.0 9.0 10-0 11.0 12.0
X
a) UB='1.0.
.30 '\ II 7' I
4
0.75-
y 0.50-
0.25m
onn‘ /
0W ' '
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10-0 11.0120
X
b) UB=-0.5.
1.00 .\ . u .
0.75—
y 0.50“
0.25—
0.00 '/ ' I 1 I I I I I I I I I
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 12.0

 

X

c) UB=1.0.

Figure 50. Development of streamwise velocity proﬁles along the channel with the
lower boundary section moving at UB; X1 = 2.0, X2 = 6.0, Re = 200.

148

 

 

 

1.00 .\ . 1
0.75;

y 0.50j
0.25:
0.001 .1”

 

 

0.0 5.0 10.0 15-0 20.0 25-0 30.0 35.0 40.0 45.0 50.0 55.0 60.0
X

Figure 51. Development of streamwise velocity proﬁles along the channel with the
lower boundary section moving opposite to the inlet velocity at UB

= -05; x1 = 10.0, X2 = 14.0, Re = 1000.

 

1.m '\ r r I

0.754

J

y 0.504

0.25-1

1

01134
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 12.0

 

 

 

X

a) UB = 0.5.11: = 200.

 

1.00 I\ ' I ‘

075‘ \

y 0.50 4

0.25 -
1

 

 

 

 

0.00
0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 22.5 25.0 27.5 30.0

X
b) UB = 1.0.11: = 500.

Figure 52. Development of streamwise velocity proﬁles along the channel with the
entire lower boundary moving in the direction of the inlet velocity
(Couette-Poiseuille ﬂow) at UB .

149

  
   

 

 

 

 

 

 

 

 

 

J.
a
l
6
j.-
in
11111111111111111

 

 

 

 

 

-
-
1— —1
-
t: -1
1— u-
p- -1
1— —1
- H
1— —I
p— -
F- -. ' —1
p—
.11
.I ’ q
. .. - -1
m‘ 4 I. " ..---. .

 

 

Figure 53. Vorticity contours with the lower boundary section moving at UB; X1 =
1.0, X2 = 2.0, Re = 20.

{é

 

 

 

 

’1

111‘111[IIIII
II
l
i
a
a
’1111 1

1

1111111;

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

a) U3 = -0050
.r- _
F- -4
C 2
1- -I
*- -1
_ d——-——.
f ' -' '_.___ ﬁ.‘ _‘
1— .1
1- .—
E a
L.
I -

 
     

Figure 54. Vorticity contours with the lower boundary section moving at UB; X1 =
1.0, X2 = 2.0, Re = 200.

 

 

—4-W

 

 

 

 

 

 

 

l I
1111 1111111 JLJ

 

 

 

 

Figure 55. Vorticity contours for the entire lower boundary moving in the direction of
the inlet velocity (Couette-Poiseuille ﬂow) at U B = 1.0; Re = 500.

151

 

a) t = 2n, Us = 1.0 (at the instant when the boundary is moving at
maximum velocity in the direction of the flow).

 

b) t = (2n+l)11:, UB = -1.0 (at the instant when the boundary is
moving at maximum velocity in the direction opposite to the ﬂow).

Figure 56. Streamlines of the ﬂow at different time with the lower boundary section
oscillating at UB = Acos(t); X1 = 1.0, X2 = 2.0, Re = 20, St = 1.0, A =

1.0, n = 8.

152

 

 

 

 

 

 

 

 

 

a) t = 201:,UB = 2.0.

 

 

 

 

 

 

 

 

b) t=(20+1)rt,UB==-2.0.

Figure 57. Streamlines of the flow at different time with the lower boundary section
oscillating at UB = Acos(t); X1 = 2.0, X2 = 4.0, Re 2 200, St = 1.0, A =

2.0, n = 25.

153

 

 

 

 

 

 

 

a) t= 2111:,0B = 2.0.

 

 

 

 

 

 

 

 

b) t=(2n+1)1t, UB = -2.0.

Figure 58. Streamlines of the flow at different time with the lower boundary section
oscillating at UB = Acos(t); X1 = 2.0, X2 = 6.0, Re = 200, St = 0.5, A =

2.0, n =- 14.

154

 

 

   

a) t = 2n1t, UB =1.0. e) t=(2n+l)11:, UB = -1.0.

 

 

   

b) t= (2n-1~l/4)1t,UB = 0.707. f) t= (2n~1~5/4)1t',UB = -0.707.

 

   

c) t = (2n+1/2)1t, UB = 0.0. g) t = (2n+3/2)1t, UB = 0.0.

 

 

d) t= (20+3/4)1c, UB = -0.707. h) t= (2n+7/4)1t, UB = 0.707.

Figure 59. Streamlines of the flow at different time with the lower boundary section
oscillating at UB = Acos(t); X1 = 1.0, X2 = 2.0, Re = 200, St = 1.0, A =

1.0, n = 21.

155

 

 

 

 

 

 

 

 

 

a) St = 0.5, n = 13.

 

 

 

 

 

 

 

 

 

b) St = 1.0, n = 21.

 

 

 

 

 

 

 

 

 

c) St = 2.0, n = 27.

Figure 60. Size of the separated region as a function of the Strouhal number; X1 =
1.0, X2 = 2.0, Re = 200, A =1.0, t = (20+1)1t, UB = -1.0.

156

 

 

 

 

 

 

 

 

 

 

a) St = 0.5, 0 =15.

 

 

 

 

 

 

 

 

 

b) St = 1.0, n = 23.

 

 

 

 

 

 

 

 

 

c) St = 2.0, n = 42.

Figure 61. Size of the separated region as a function of the Strouhal number; X1 =
2.0, X2 = 4.0, Re = 200, A = 1.0, t = (20+ l)1t, UB = -1.0.

157

 

 

 

 

 

 

 

 

 

 

a) St = 0.5, n =16.

 

 

 

 

 

 

 

 

 

 

b) St = 1.0, n = 26.

 

 

 

 

 

 

 

 

 

 

c) St = 2.0, n = 45.

Figure 62. Size of the separated region as a function of the Strouhal number, X1 =
2.0, X2 = 6.0, Re = 200, A = 0.5, t = (2n+1)1t, UB = -0.5.

 

 

a) A =1.0, n = 8, UB = -1.0.

 

b) A = 2.0, n = 8, UB = -2.0.

Figure 63. Flow streamlines showing the size of the separated region as a function
of the amplitude of oscillation of the boundary; X1 = 1.0, X2 = 2.0, Re =

20, St = 2.0, t = (2n+1)1t.

159

 

 

a) A = 0.5, n = 25, UB = -0.5.

 

 

b) A = 2.0, n = 20, 0,, = -2.0.

Figure 64. Size of the separated region as a function of the amplitude of oscillation
of the boundary; X1 = 1.0, X2 = 2.0, Re = 200, St = 1.0, t = (20+1)1c.

160

 

 

 

 

 

 

 

 

 

a) A = 0.5, n = 46, UB = -0.5.

 

 

 

 

 

 

 

 

b) A = 2.0, n = 25, U13 = -2.0.

Figure 65. Size of the separated region as a function of the amplitude of oscillation
of the boundary; X1 = 2.0, X2 = 4.0, Re = 200, St = 2.0, t = (2n+1)1t.

161

 

 

 

 

 

 

 

 

 

a) A = 1.0.11 = 24, UB = -1.0.

 

 

 

 

 

 

 

b) A = 2.0.11 = 25, 0B = -2.0.

Figure 66. Size of the separated region as a function of the amplitude of oscillation
of the boundary; X1 = 2.0, X2 = 6.0, Re = 200, St = 1.0, t = (2n+1)11:.

162

 

 

 

 

 

 

 

 

 

a) X1=1.0,X2=2.0,St=1.0,n = 68.

 

 

 

 

 

 

 

 

 

b) x1 = 10.0, x2 =14.0,s:= 1.0.11 =101.

 

 

 

 

 

 

 

 

 

0) x1 = 10.0, x2 = 14.0.31: 0.5, n = 56.

Figure 67. Streamlines of the flow with the lower boundary section oscillating at UB
= Acos(t); Re =1000, A = 1.0, t = (20+1)rt, U3 = -1.0.

163

 

 

 

 

 

 

 

 

 

Figure 68. Effect of the location of the oscillating section of the boundary on the size
of the separated region; X1 = 5.0, X2 = 6.0, Re = 200, St = 1.0, A = 1.0, t

= (2n+l)1r, n = 28, UB = -1.0.

 

 

 

 

 

 

 

Figure 69. Effect of the length of the oscillating section of the boundary on the size
of the separated region; X1 = 1.0, X2 = 3.0, Re = 200, St = 1.0, A = 1.0, t

= (20+1)x, n = 26, U], = -1.0.

 

a) St = 0.5, A = 2.0, n = 17, U3 = -2.0.

 

b) St = 1.0, A =1.0, n = 24, U3 = -l.0.

 

c) St = 2.0, A = 0.5, n = 60, U3 = -0.5.

Figure 70. Effect of the location and length of the oscillating section of the boundary
on the size of the separated region; X1 = 8.0, X2 = 10.0, Re = 200, t =

(2n+ Du.

.
’1

We 7.

165

 

 

 

1.00 .\ . . , .
0.75:
y 0.50:
0.sz
0.001 ./ -

 

 

0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10-0 11.0 12.0

X

a) XI: 2.0, x2 = 4.0, s: = 1.0, A =1.0,n = 23, UB =1.o.

 

 

1.00 .\ ' 1
0.753
y 0.50:
0.25;
0.00‘ ./ - -='—~-

 

 

 

0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 12.0

X

b) x1 = 2.0, x2 = 4.0, s: = 0.5, A = 2.0, n = 13, UB = 2.0.

 

 

 

1.00 .\ . r .
0.75:

y 0.50:
0.251
0.001 ./ . . -

 

 

0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10-0 11.0 12.0

X

0) x1: 2.0,X2=6.0, St= 1.0, A = 1.0, n =24, UB =10.

Figure 71. Development of streamwise velocity proﬁles along the channel at the
instant when the boundary is moving at maximum velocity in the

direction of the flow; Re = 200, t = 2n1r.

166

 

 

1.00 .\
0.75;
y 0.50:
0.25:
‘ Z

0.00 i
0.0 1.0 2.0

 

 

X

a) xl =2.0,x2 =4.0, St= l.0,A = 1.0.11 =23, UB = -1.0.

 

tT70 rBrio 10-0 11.0 12.0

 

1.00 .\

1

0.75-1
i
y 0.50 '-

0.25 -

4

0.00-
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0

 

 

X

 

10.0 11.0 12.0

b) x1 = 2.0, x2 = 4.0, s: = 0.5, A = 2.0, n =13,UB = .20.

 

1.00 N 12 T
0.751
y 0.50;
0.25-j
0.004 ./ '

 

 

 

0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0

X

 

10.0 11.0 12.0

0) XI: 2.0, x2 = 6.0, s: = 1.0, A =1.0,n = 24, UB = -1.0.

Figure 72. Development of streamwise velocity proﬁles along the channel at the

instant when the boundary is moving at maximum
direction opposite to the ﬂow; Re = 200, t = (2n+1)1c.

velocity in the

   

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

11

«0.1

0.0 g_g_____-_._— 0.! :

.0 Z

I _.
1

N. .- .a. m a .m.....---____-

a) t= 2nrt, UB =1.0.

0 0 0 D —0.0 7

1', 12 0 3

"M... ______ , m.- .. a. '4

b) t= (2n+1)1t, UB = -1.0.

Figure 73. Vorticity contours for the case of the lower boundary section oscillating
at UB = Acos(t);X1=1.0, X2 = 2.0, Re = 20, St = 1.0, A =1.0, n = 8.

168

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

‘_ I

1— n-

P —I

— -

h _

I— .1

o a.-—_ 4.0 g g

In- —

1— a—

E —

r- _ ..
F A ‘—

a) t= 2nx’ UB = 100.

ﬁ

14

—1

I

-1

w

.4

-1

J x '—

' A —

"I§.‘s u _-

b) t=(2n+1)1t, UB = -1.0.

Figure 74. Vorticity contours for the case of the lower boundary section oscillating
at UB = Acos(t);Xl=1.0, X2 = 2.0, Re = 200, St = 1.0, A =1.0, n = 21.

169

 

 

 

 

 

 

 

    
   

 

 

 

     

P- I-
_ -
- -
1- — d
1— _.
1— .4
1- —1
_ I
- -1
- d
- i -1
E ' l ‘ 1a :
I ’7?,73.35:"I."'"‘,'.:—.‘-‘Ei‘:‘?::.::r:::;:::42553;.“1- Tull-n- ........ -

 

a) t= 2n1r, UB =10.

 

 

 

 

 

 

 

 

 
     

 

 

 

 

P -
1- -1
- -1
F
_ _
_ ‘
0.0 0.0 90.0-----
b -
— ‘
E _ -
1- .1
P -
w
b I-

b) t=(20-1-1)1c,UB =-1.0.

Figure 75. Vorticity contours for the case of the lower boundary section oscillating
at U13 = Acos(t); X1 = 8.0, X2 =10.0, Re = 200, St = 1.0, A = 1.0, n = 24.

APPENDIX

00000000000000000000OOOOOOOOOOOOOOOOOOOOOOO

A

 

COMPUTER PROGRAM
PROGRAM ENTRY
File Name: ENTRY.FOR
Written By: Bashar S. AbdulNour - December 1988

 

A computer program for the numerical solution of the vorticity
transport equation and stream function equation describing the time-
dependent, viscous, incompressible flow in the entry region of a
two-dimensional channel with part of the lower boundary oscillating.
A finite-difference scheme, using the Alternating-Direction Implicit
(ADI) method to solve the vorticity transport equation for the
vorticity and the Successive Over-Relaxation (SOR) method to solve
Poisson's equation for the stream function, is employed.

 

DESCRIPTION OF VARIABLES

M .
M1 -
M2 -
N a

I, J =

number of grid points in streamwise direction (X-dir.);
along channel's length

number of grid points in X-dir. to the left of the
oscillating boundary

number of grid points in X-dir. to the right of the
oscillating boundary

number of grid points in normal direction (Y—dir.);
along channel’s height

subscripts corresponding, respectively, to nodal points
in the x- and Y-directions

I=1 is the left boundary

12M is the right boundary

J=1 is the lower boundary

J=N is the upper boundary

J-NC is the centerline

[domain of interior reigon : I=2 to M-l, J-2 to N-l]
mesh size; equal in x- and Y—directions

Reynolds number based on channel height and inlet
velocity

Strouhal number based on channel height and inlet
velocity [- WL/U, where W is frequency of oscillation]

170

000000OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO00000000000000000000000

INC '
KT -
KK,KKA(T)'
T .
TA(T)
TKK
KTP
KINC
AA
PHI

PHIA(T)
X(I)

Y(J)
XR(I)
F(I,J)
Q(I.J)
FS(I,J)
QS(I,J)
FC(I,J,T)

QC(IIJIT)-

FT(I,J,T)-

QT‘IIJIT)-

RFS
AKS
AKV
BKV
ITMAX

ITER

ITERA(T)
ITMXF -

ITMXQ

ITERF
ITERQ
EPSILON
EPSILONF
EPSILONQ
ERFO,

ERQO
ERF, ERQ

DERFO,
DERQO

KOUNT

NCYCLE
U(I,J)
V(I,J)
U0(I,T)
ILE

LE

171

number of time increments in one time period

time step used in the ADI method

counters for time steps

time

array to store time (when convergence begins)

time from start of the 2*PI period

time fraction of the 2*PI period

lable for time steps included in the convergence period
amplitude of oscillation of boundary velocity

input forcing function (i.e., velocity of the
oscillating boundary which is a given function of time)
array to store values of PHI at successive time steps
streamwise (i.e., axial) coordinate [X=0 to (M-1)*H]
normal (i.e., transverse) coordinate [Y=0 to 1]

X(I)/RE

calculated values of stream function

calculated values of vorticity

stored values of stream function

stored values of vorticity

values of stream function at selected locations used to
check the convergence

values of vorticity at selected locations used to check
the convergence

values of stream function at specified time during the
last period before convergence; values are used in
output

values of vorticity at specified time during the last
period before convergence; values are used in output
over—relaxation factor for Poisson's equation
weighting factor for stream function

weighting factor for vorticity in interior region
weighting factor for vorticity boundary conditions
preselected maximum number of outer iterations (or time
steps)

number of outer iterations for the stream function and
vorticity

array to store ITER (when convergence begins)

maximum permitted number of inner iterations for the
stream function

maximum permitted number of inner iterations for the
vorticity

number of inner iterations for the stream function
number of inner iterations for the vorticity

tolerance error (infinitesimally small quantity)

stream function tolerance error

vorticity tolerance error

maximum error in stream function and vorticity,
respectively, for the outer iteration

maximum error in stream function and vorticity,
respectively, for inner iterations

denominators of the relative error terms for the stream
function and vorticity, respectively, used in the outer
iteration to avoid division by zero

counter to insure that the outer iteration is
convergant for consecutive iterations over one full
period

number of cycles needed for the solution to converge
streamwise velocity (in the X-direction); u-U(X,Y)
normal velocity (in the Y-direction); V=V(X,Y)

velocity of the boundary (a function of time)

grid point along the centerline representing the
entrance length

entrance length defined when the centerline velocity
reaches 99% of its non-dimensional value for a fully-

000 0000OOOOOOOOOOOOOOOOOOOOOOO

GOO GOO

GOO

1T2

developed parabolic velocity profile

LERE - LE/RE ~

MSKIP - counter to skip X-location data points in order to
reduce printing of output at high RE

MMl, NMl, PI, ETA, H2, TH, C1, C2, GI, 62:
are constants defined throughout the program to

facilitate the computation. They are defined as they
occur in the program.

DESCRIPTION OF SUBROUTINES

Subroutine Function

 

Q_ADI Solving the vorticity equation using the ADI method

F_COL Solving the stream function equation using the SOR by
columns method

TRI( ) Solving a system of equation with tri-diagonal matrix
using Thomas algorithm (it is used within other
routines)

MORE( ) Writing the complete results at different time into the

output file

 

Declaration of Variables (double precision for all real variables):

 

INCLUDE 'INCLUDE.FOR’ ! include file INCLUDE.FOR to define
arrays dimensions IDIM, JDIM, and KDIM;
declare variables; and define common blocks

REAL*8 AKS, AKV, ETA

REAL*8 FC(IDIM10,JDIM,KDIM), QC(IDIM10,JDIM,KDIM)
REAL*8 FT(IDIM,JDIM,KDIM8), QT(IDIM,JDIM,KDIM8)
REAL*8 T, AA

REAL*8 UO(IDIM10,KDIM)

REAL*8 ERFO, DERFO

REAL*8 ERQO, DERQO

REAL*8 LE, LERE

INTEGER INC, KK, KINC

INTEGER KOUNT, KSKIP, KNDEX, K

INTEGER NCYCLE

INTEGER ILE

INTEGER LUNITR, LUNITS, LUNITW

Define Logical Unit Numbers:

 

LUNITR - 12 ! logical unit number of input data file
LUNITS - 6 ! logical unit number of screen (default output)
LUNITW - 14 ! logical unit number of output file

Read Input Data from Datafile:

 

OPEN (LUNITR,FILE-'ENTRY.DAT',STATUS='UNKNOWN')

READ (LUNITR,*) M,N,M1,M2,RE,ST,INC,RFS,AKS,AKV,BKV,

& EPSILON,EPSILONF,EPSILONQ,ITMAX,ITMXF,ITMXQ
CLOSE (UNIT-LUNITR)

Compute Grid Points and Coordinates:

 

MMl - M-l

NMl - N-l

NC (N+1)/2 ! label centerline nodes
H - l.0/FLOAT(N-1)

X(1)

173

- 0.

XR(1) - 0.
DO 9110 I-2,M

X(I)

- X(I-1) + H

9110 XR(I) - X(I)/RE

Y(1)

- 0

no 9120 3-2,n

9120 Y(J)

for the

0000

- Y(J-1) + H

Calculate the Optimum Value of the Over—Relaxation Factor

Stream Function:

 

PI -
ETA -
RFS -

000

4.*DATAN(1.D0)
0.5*(COS(PI/MM1) + COS(PI/NM1))
2.*(l.-SQRT(1.-ETA**2))/ETA**2

Introduce and Calculate Some Constants and Parameters:

 

0
H
IIIIIII

COMPUTE

2.*PI/INC
H*H

2.*H
RPS/4.

1. - RFS
KT/(RE*H)
2.*Gl

BOUNDARY CONDITIONS

***********i***************
Stream Function Boundary Conditions:

 

l. Upstream - Uniform Profile at Inlet

 

0000000

DO 9215 J-2,NM1

F(1,J)

‘0
N
H
U‘

- Y(J)

FS(1,J) - Y(J)

C 2. Downstream

 

DO 9220 J-2,NM1

 

 

 

 

F(M,J) - 3.*(Y(J)**2) - 2.*(Y(J)**3)
9220 FS(M,J) - F(M,J)
C
C 3. Lower and Upper Boundaries
C
DO 9230 I-1,M
F‘I’l) - O.
FS(I,1) - 0.
F(I,N) - 1.0
9230 FS(I,N) - 1.0
C
C Vorticity Boundary Conditions:
C
C l. Upstream - Uniform Profile at Inlet
C
DO 9255 J-1,N
Q‘lrJ) - 0.
9255 QS(1,J) - 0.
C
C 2. Downstream
C

DO 9260 J-1,N

Q(M,J)

- 12.*Y(J) - 6.

9260 QS(M,J) 3 Q(M,J)

174

3. Lower and Upper Boundaries [to be evaluated later since they
depend on Q(I,2) and Q(I,N-l), respectively]

 

Initialization of Values on Grid of Interior Region:

 

0()0(30(3C)

DO 9400 I-Z,MM1
DO 9400 J-2,NM1
F(I,J) - Y(J)
9400 Q(I,J) - 0.
C

C Initialization of Values of F and Q Selected to Test Convergence:
C

 

DO 9420 I-11,M-10,10
DO 9420 J-2,NM1
DO 9420 KK-1,INC
II - (I - 1)/10
FC(II,J,KK) - 0.
9420 QC(II,J,KK) - 0.
C
C Initialization of Values of F and 0 Selected for Output:
C

 

DO 9422 I-2,MM1
DO 9422 J-2,NM1
D0 9422 K=1,KDIM8
FT(I,J,K) - o.

9422 QT(I,J,K) - o.

c

C BEGIN OUTER ITERATION FOR THE STREAM FUNCTION AND VORTICITY

C ***t*******************************************************
KR - O
KOUNT 0
KSKIP INC/8 ! number of time steps skipped
KNDEX ! time step index for FT and QT
K ! same as KNDEX
ITERF
ITERQ

(DOCDO

DO 9590 ITER=1,ITMAX

IF (KK.EQ.INC) KK - 0
RR - KK + 1

T - ITER*KT

PHI - AA*DCOS(T)

Storing All Old Values of the Stream Function and Vorticity
(Values from the Previous Time Step):

0(3C)0

 

DO 9450 I=Z,MM1

DO 9450 J82,NM1

FS(I,J) - F(I,J)
9450 QS(I,J) - Q(I,J)
C

SOLVING THE VORTICITY TRANSPORT EQUATION FOR THE VORTICITY
*******t**************************************************

3. (Continued) Calculate New Values of Lower and Upper
Vorticity Boundary Condition at Begining of Vorticity Inner
Iteration [at new time step] Using the Second-Order Woods’
Formula with Values of the Stream Function and Vorticity
Known from Previous Iteration

 

C)0(1()0(7()0

D0 9512 I'Z,M1-1
QS‘Ill) - Q‘Irl)

9512
9514

9516
9520
C
C
C
C
C
C
C
C
C

9530
9540

550
560

0000090 0000000 00000

175

Q‘Irl) - 3-*(F(I11)’F(Ilz))/Hz " Q(IIZ)/2o
Q(I,1) I BKV*QS(I,1) + (1.-BKV)*Q(I,1)

QS(M1I1)- Q‘er 1)

Q(Ml,1) I 3.*F(M1,1)/HZ - 3.*(F(M1-l,2)+F(M1+1,2))/(2.*H2)
- Q(Ml,2)/2. + 3.*PHI/TH

Q(Ml,1) I BKV*QS(M1,1) + (1.-BKV)*Q(M1,1)

DO 9514 IIM1+1,M2-1

QS(II1) ' Q‘Irl)

Q(I,l) - 3.*(F(I,l)-F(I,2))/H2 - Q(I,2)/2. + 3.*PHI/H
Q‘Ill’ ' BKV*QS(II1) + (1.-BKV)*Q(I,1)

QS(M2.1)- 0(M2.1)

Q(M2,l) - 3.*F(M2,l)/H2 - 3.*(F(M2-1,2)+F(M2+1,2))/(2.*H2)
- Q(M2,2)/2. + 3.*PHI/TH

Q(M2.1) - BKV*QS(M2.1) + (1.-BKV)*Q(M2.1)

DO 9516 IIM2+1,MM1

QS(I.1) - 0(I.1)

0(I.1) - 3.*(F(I,1)-F(I.2))/H2 - 0(I,2)/2.
Q(I,1) - BKV*QS(I,1) + (l.-BKV)*Q(I,1)

DO 9520 II2,MM1

QS(IIN) ' Q(IIN)

Q‘IIN) ' 3-*(F(I:N)'F(I:NM1))/HZ - Q(IINM1)/20
C(IrN) ' BKV*QS(I:N) + (1.-BKV)*Q(I,N)

Calling Vorticity Equation Solver:

 

CALL Q_ADI

End of inner iteration for the vorticity.

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA

Recalculate Q(I,J) Using Weighting Factor:

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

DO 9540 II2,MM1
DO 9540 JI2,NM1
Q(I,J) I AKV*QS(I,J) + (1.-AKV)*Q(I,J)

SOLVING POISSON EQUATION FOR THE STREAM FUNCTION
t***********t***********************************

Calling Poisson's Equation Solver:

 

CALL F_COL

End of inner iteration for stream function.

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA

Recalculate F(I,J) Using Weighting Factor:

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

DO 9560 II2,MM1
DO 9560 JI2,NM1
F(I,J) I AKS*FS(I,J) + (1.-AKS)*F(I,J)

Test the Values of the Stream Function and Vorticity
Recalculated from the Outer Iteration for Convergence
(by comparing selected values over one time period):

 

ERFO I 0.7

ERQO I 0.

DO 9580 II11,M-10,10
DO 9580 JI2,NM1

176

II I (I - 1)/10
DERFO I F(I,J) + EPSILON

DERQO I Q(I,J) + EPSILON

ERFO I DMAX1(ERFO,DABS((F(I,J)-FC(II,J,KK))/DERFO))
9580 ERQO I DMAX1(ERQO,DABS((Q(I,J)-QC(II,J,KK))/DERQO))
C

IF ((ERFO.LE.EPSILONF).AND.(ERQO.LE.EPSILONQ)) THEN
IF ((KK-l).EQ.0) KK I INC + 1
IF ((KOUNT.GT.0).AND.(ITER-ITERA(KK-1)).NE.1) KOUNT I 0
IF (KK.GT.INC) KK I 1
9581 KOUNT I KOUNT + 1

ITERA(KK) I ITER
KKA(KK) I KK
TA(KK) I T
PHIA(KK) I PHI

DO 9584 IIM1,M2,10
II - (I - 1)/1o
9584 U0(II,KK) - pa:

IF (KOUNT.GE.KK) THEN
KNDEX I INT(KK/KSKIP)
IF (KNDEX.EQ.(K+1)) THEN
K I KNDEX
DO 9586 I-1,M ' storing values of the stream
DO 9586 JIl,N ! function and vorticity of
FT(I,J,K) - F(I,J) ' the last full period before
9586 QT(I,J,K) I Q(I,J) ' exiting
ENDIF
ENDIF

IF ((KOUNT.LT.INC).OR.(KK.LT.INC)) GO TO 9588
GO TO 9600
ENDIF

Storing Selected Values of the Stream Function and Vorticity at
the End of Each Time Step (to be Used to Check for Convergence):

 

900000

588 D0 9589 I=11,M-10,10
DO 9589 J=l,N
II - (I - 1)/1o
FC(II,J,KK) I F(I,J)
9589 QC(II,J,KK) I Q(I,J)
C
9590 CONTINUE
C
C Error Message if Number of Outer Iterations Exceed ITMAX:
C

 

IF(ITER.GE.ITMAX) THEN
9980 WRITE (LUNITS, 9981) ITMAX
9981 FORMAT (/, 2X,'ERROR ,/, 2x, ----- ',
5 /, 2X,’JOB was ABORTED DUE To OUTER ITERATION PROBLEM.
& /, 2X,'NUMBER OF OUTER ITERATIONS ',
& 'HAS EXCEEDED ITMAX - ',I4,/)
ENDIF

Calculate the Number of Cycles for Convergence:

 

600 NCYCLE I ITER/INC

End of outer iteration.

AAAAAAAAAAAAAAAAAAAAAAA

0000\0000

177

Inlet and Downstream End:

000

 

U(1,1) I U(l,NC)
U(l,N) I U(l,NC)

Open Output File:

 

WRITE (LUNITW,9801)
801 FORMAT (/.1x,'F11e Name: ENTRY.OUT',/)

Write Data and Calculated Results into Output File:

000“) 0 000

 

WRITE (LUNITW,9810)
9810 FORMAT (/,1x,'Data:',

Compute the Streamwise Velocity on the Lower Boundary at the

a /,1x,'---- ',/)
C
WRITE (LUNITW,9811) RE,ST
9811 FORMAT (2X,’REYNOLDS NUMBER, RE - ’,F10.1,/,
& 2X,’STROUHAL NUMBER, ST - ',F12.3,/)
C
WRITE (LUNITW,9815) M,N,H
9815 FORMAT (2X,’TOTAL x-DIR. GRID POINTS, M - ',I8 ,/,
& 2x,'TOTAL Y-DIR. GRID POINTS, N - ’,I8 ,/,
& 2X,'MESH SIZE IN x- AND Y-DIRECTION, H - ',F12.3 )
C
WRITE (LUNITW,9817) KT
9817 FORMAT (2x,'TIME STEP, KT - ',F12.3,/)
C

WRITE (LUNITW,9822) RFS

OPEN (LUNITW,FILEI’ENTRY.OUT’,STATUSI'NEW',FORM='FORMATTED')

9822 FORMAT (2X,'STREAM FN OVER-RELAXATION FACTOR, RFS - ',F7.3)
WRITE (LUNITW,9824) AKS,AKV,BKV

9824 FORMAT (2X,'STREAM FN WEIGHTING FACTOR, AKS - ',F7.3,/,
& 2X,'VORTICITY WEIGHTING FACTOR, AKV - ',F7.3,/,
& 2X,'VORTICITY B.C. WEIGHTING FACTOR, BKV - ’,F7.3,/)

C
WRITE (LUNITW,9831) EPSILONF,EPSILONQ

9831 FORMAT (2x,'STREAM FN PERMITTED ERROR, EPSILONF - ',E14.6,/,
& 2X,'VORTICITY PERMITTED ERROR, EPSILONQ - ',El4.6,/)

C
WRITE (LUNITW,9833) ERF,ERQ,ERFO,ERQO

9833 FORMAT (2X,’STREAM FN ERROR (INNER ITERATION), ERF - ',El4.6,/,
& 2X,’VORTICITY ERROR (INNER ITERATION), ERQ - ',E14.6,/,
& 2X,'STREAM FN ERROR (OUTER ITERATION), ERFO - ’,E14.6,/,
& 2X,'VORTICITY ERROR (OUTER ITERATION), ERQO - ',El4.6,/)

WRITE (LUNITW,9842) ITMAX,ITMXF,ITMXQ,ITERF,ITERQ,ITER
9842 FORMAT (
2X,’PERMITTED MAXIMUM NO. OF ITERATIONS, ITMAX
2X,’PERMITTED MAX. NO. OF STREAM FN ITERATIONS, ITMXF
2X,’PERMITTED MAX. NO. OF VORTICITY ITERATIONS, ITMXQ
2X,'FINAL ACTUAL NO. OF STREAM FN ITERATIONS, ITERF
2X,’FINAL ACTUAL NO. OF VORTICITY ITERATIONS, ITERQ
2X,’ACTUAL NO. OF OUTER ITERATIONS, ITER

hhhhhh

WRITE (LUNITW,9846) NCYCLE
9846 FORMAT (
& 2X,'NUMBER OF CYCLES FOR CONVERGENCE, NCYCLE

WRITE (LUNITW,9350)
9850 FORMAT (1X,'Results:’,/,
8 1x" _______ "/)

\D

9!

986

966

9872

9877

9875
9879

(Itvnnnn
lﬂl'l‘ﬂlf"

Q
~~

 

9855

9859

9870

9867

9861

9872

9877

9875
9879

(UC)0<1()0

9886
C
9888
C

C
C
C

178

WRITE (LUNITW,9855) X(M)
FORMAT (2X,’TOTAT LENGTH OF CHANNEL - ’,F8.4 ,/)

WRITE (LUNITW,9859)
FORMAT (’1',/)

WRITE (LUNITW,9870)
FORMAT (2X,'VORTICITY AT LOWER AND UPPER BOUNDARIES FOR ',
'ONE FULL PERIOD BEFORE ',/,
2X,'EXITING ',
'(ONE FULL PERIOD AFTER START OF CONVERGENCE) ',//)

bl!”

DO 9879 KINCI1,INC

TKK - KKA(KINC)*KT

KTP - TKK/(2.*PI)

WRITE (LUNITW,9867)

FORMAT (2x,' ',
I I)

WRITE (LUNITW,9861) KKA(KINC),ITERA(KINC),TA(KINC),TKK,KTP

FORMAT (/,2x,'KK - ',16 ,/,

 

 

&

& 2X,'ITER - ',16 ,/,
& 2x,'TIME, T - ',F9.4,/,
& 2X,'TKK - ',F9.4,/,
& 2X,'KTP - ',F9.4,/)
WRITE (LUNITW,9872)
FORMAT (2x,' x ' ,
5 5X:' U(XIOIT) 'IZXI' Q(X:0:T) ',2X,' Q(X:H,T) ',/)

WRITE (LUNITW,9875) X(1),U(l,1),Q(l,l),Q(l,N)

DO 9877 II11,M-10,10

II - (I - 1)/1o

WRITE (LUNITW,9875) X(I),U0(II,KINC),QC(II,1,KINC),QC(II,N,KINC)
WRITE (LUNITW,9875) X(M),U(M,1),Q(M,1),Q(M,N)

FORMAT (2x,F7.3,2x,3(2x,F10.6))

CONTINUE

Call Subroutine MORE to Continue Writing Results (X, Y, U, V, F,
and Q) at Selected Time Steps (Within One Full Period Before
Exiting) into Output File; Values of Stream Function and Vorticity
are Reassigned First:

 

DO 9838 K=1,KDIM8
KINC ' K*KSKIP
DO 9886 II1,M

DO 9886 J=I,N
F(I,J) I FT(I,J,K)
Q(IIJ) ' QT(IIJIK)

CALL MORE (ILE,KINC,LUNITW)

Close Output File:

 

CLOSE (UNIT-LUNITW)

STOP
END

 

**

********************************************************************

 

00000000

0 0000000000000000000000000000000000000000000000000

0000

179

SUBROUTINE Q_ADI

 

A subroutine for solving vorticity transport equation for the
vorticity using the Alternating-Direction Implicit (ADI) method.
The ADI method splits each time step into two halves. For the
first half, the discrete form of the vorticity equation is
implicit in X. For the second half the difference equation is
implicit in Y. Each equation yields a tri-diagonal system of
simultaneous linear algebraic equations. The resulting systems
are solved using Thomas algorithm which is a special adaptation of
the Gaussian elimination. This process is repeated over
successive time steps.

 

DESCRIPTION OF VARIABLES EXECLUSIVE TO THIS ROUTINE

old values of the vorticity from the previous inner

iteration

intermediate values of the vorticity at the end of

first half of time step

U(I,J)*K/2

V(I,J)*K/2

finite-difference multiplier used to convert the

central difference approximation of the first order

derivatives of the vorticity with respect to

x, contained in the convective terms, to forward/

backward difference approximation

ALPHAl-l give backward difference approximation

ALPHAlI—l give forward difference approximation

ALPHA2( ) same as above but for derivatives with respect to Y

Al(I), I one-dimensional arrays containing elements of the

A2(I), upper, main, and lower diagonals, respectively, of

A3(I) the tri-diagonal matrix of the vorticity equation;
it applies to both the X-equation (first half of
time step) and the Y-equation (second half of time
step)

BB(I) I one-dimensional array containing elements of the
right-hand-side column matrix of the vorticity
equation

WW(I) I column matrix of the solution of the vortivity
equation

DERQ I denominator of the relative error term for the

vorticity; used in inner iteration to avoid division

by zero

QOLD(I:J7

QINT(I,J)

UDT(I,J)
VDT(I,J)
ALPHA1( )

 

INCLUDE 'INCLUDE.FOR’ ! include file

REAL*8 QOLD(IDIM,JDIM), QINT(IDIM,JDIM)
REAL*8 UDT(IDIM,JDIM), VDT(IDIM,JDIM)
REAL*8 ALPHA1(IDIM,JDIM), ALPHA2(IDIM,JDIM)
REAL*8 DERQ

Calculate the Velocities and difference operators at All
Interior Nodes:

 

DO 8200 I=2,MM1

DO 8200 J-2,NM1

U(I,J) - (F(I,J+1)-F(I,J-1))/TH
IF (U(I,J).GT.0.0) THEN

3
C
CC
GUCCCC
C
C
(c

84C
55C

Ab
ECCC
Ch»

CC

8200

0 ()0(3

250

0(36)0(3()m

0

8300

C
C
8400
C
8500

0(1()0

8550

0(0630

180

ALPHA1(I,J) I 1.0
ELSEIF (U(I,J).LT.0.0) THEN
ALPHA1(I,J) I -1.0
ELSE
ALPHA1(I,J) ' 0.0
ENDIF
UDT(I,J) I U(I,J)*KT/2.

V(I,J) - (F(I-1,J)—F(I+1,J))/TH
IF (V(I,J).GT.0.0) THEN
ALPHA2(I,J) - -1.o
ELSEIF (V(I,J).LT.0.0) THEN
ALPHA2(I,J) - 1.0
ELSE
ALPHA2(I,J) - 0.0
ENDIF
VDT(I,J) - V(I,J)*KT/2.

Begin ADI Iteration to Solve the Vorticity Equation:

 

DO 8900 ITERQI1,ITMXQ
DO 8250 II2,MM1

DO 8250 J=2,NM1
QOLD(I,J) ' Q(IIJ)

COMPUTE THE VORTICITY FOR INTERIOR REGION:

 

l. X-Equation for the First Half of Time Step;
Solve for the Intermediate Value of the Vorticity

 

DO 8500 JI2,NM1
DO 8300 I=2,MM1

A1(I) I -(1+ALPHA1(I,J))*UDT(I,J) - Gl

A2(I) I (2*ALPHA1(I,J))*UDT(I,J) + 62 + TH*ST

A3(I) I (l-ALPHA1(I,J))*UDT(I,J) - 61

38(1) I (( (1+ALPHA2(I,J))*VDT(I,J) + 61) *Q(I,J-l)
+(—(2*ALPHA2(I,J))*VDT(I,J) - GZ + TH*ST)*Q(I,J)
+(-(1-ALPHA2(I,J))*VDT(I,J) + 61) *Q(I,J+l))

CONTINUE

BB(2) I BB(2) - Al(2)*Q(1,J)

A1(2) I 0.

BB(MMI) I BB(MMl) - A3(MM1)*Q(M,J)

A3(MM1) I 0.

CALL TRI (2,MM1) ! tri-diagonal system solver;

solves for Q at fixed J

DO 8400 II2,MM1
QINT(I,J) I WW(I)

CONTINUE

2. Rename the Intermediate Values of the Vorticity X-Boundary
Conditions at End of First Half of Time Step

 

DO 8550 J31,N
QINT(1,J) ' Q(1aJ)
QINT (“I J) - Q (“I J)

3. Compute the Stream Function and recalculate the Velocities
and difference operators at All Interior Nodes:

 

8580
C

C
C
C

00

8600
C

C
C
8700
C

8800
C

0 000

8850

8900

&
&

181

CALL F_COL

DO 8580 II2,MM1
DO 8580 JI2,NM1
U(I,J) I (F(I:J+1)-F(I,J-1))/TH
IF (U(I,J).GT.0.0) THEN
ALPHA1(I,J) - —1.0
ELSEIF (U(I,J).LT.0.0) THEN
ALPHAl(I,J) - 1.0
ELSE
ALPHA1(I,J) - 0.0
ENDIF
UDT(I,J) - U(I,J)*KT/2.

V(I,J) - (F(I-1,J)-F(I+1,J))/TH
IF (V(I,J).GT.0.0) THEN
ALPHA2(I,J) - 1.0
ELSEIF (V(I,J).LT.0.0) THEN
ALPHA2(I,J) - -1.o
ELSE
ALPHA2(I,J) - 0.0
ENDIF
VDT(I,J) - V(I,J)*xT/2.

4. Y-Equation for the Second Half of Time Step;

Compute a New Value for the Vorticity at End of Time Step

 

DO 8800 II2,MM1
DO 8600 JI2,NM1

A1(J) I -(1+ALPHA2(I,J))*VDT(I,J) - G1

A2(J) I (2*ALPHA2(I,J))*VDT(I,J) + GZ + TH*ST

A3(J) - (l-ALPHA2(I,J))*VDT(I,J) - Gl

BB(J) I (( (1+ALPHA1(I,J))*UDT(I,J) + G1) *QINT(I-1,J)
+(-(2*ALPHA1(I,J))*UDT(I,J) - 62 + TH*ST)*QINT(I,J)
+(-(1-ALPHA1(I,J))*UDT(I,J) + G1) *QINT(I+1,J))

CONTINUE

BB(2) I BB(2) — A1(2)*Q(I,1)

A1(2) I 0.

BB(NM1) - BB(NM1) - A3(NM1)*Q(I,N)

A3(NM1) - 0.

CALL TRI (2,NM1) ! tri-diagonal system solver;

solves for Q at fixed I
DO 8700 JI2,NM1
Q(I,J) - WW(J)
CONTINUE

IF (ITMXQ.EQ.1) RETURN

Test the Vorticity for Convergence:

 

ERQIO.

DO 8850 II2,MM1

DO 8850 JI2,NM1

DERQ I Q(I,J) + EPSILON

BRO " DMAX1(ERQ,DABS( (Q(I:J)-QOLD(I,J) ) /DERQ) )
IF (ERQ.LE.EPSILONQ) GO TO 8950

CONTINUE

182

\O
U
C)

RETURN
END

 

**********************************************************************

 

SUBROUTINE F_COL

 

A subroutine for solving Poisson's equation for the stream
function using the Successive Over-Relaxation (SOR) method. The
difference form of the stream function equation yields a system
of simultaneous linear algebraic equations. The SOR by blocks
(or lines) method is used to solve the resulting system of
equations. Blocks are chosen to be columns of the computational
domain.

 

DESCRIPTION OF VARIABLES EXECLUSIVE TO THIS ROUTINE

FOLD(I,J) I stored values of the stream function (from previous

iteration)
AI(J), I one-dimensional arrays containing elements of the
BI(J), upper, main, and lower diagonals, respectively, of
CI(J) the tri-diagonal matrix of Poisson's equation
(representing nodes J in column I)
DI(J) I one-dimensional array containing elements of the
right-hand-side column matrix of Poisson’s equation
WI(J) I column matrix of the solution of Poisson's equation
DERF I denominator of the relative error term for the

stream function; used in inner iteration to avoid
division by zero

 

INCLUDE ’INCLUDE.FOR' ! include file
REAL*8 DERF

Begin SOR Iteration to Solve Poisson's Equation:

 

DO 5700 ITERF-1,ITMXF

0 000 0 000000000000000000000000000000 00000 (:30

DO 5300 II2,MM1
DO 5300 JI2,NM1
5300 FOLD(I,J) I F(I,J)

 

C
C COMPUTE THE STREAM FUNCTION FOR INTERIOR REGION:
C Start Solution of Block-Iteration by Columns and SOR:
C
DO 5500 II2,MM1
C
AI(2) I 0.
81(2) I -4.
CI(2) I l.
DI(2) - -H2*Q(I,2) - F(I-1,2) - F(I+1,2) - F(I,1)

DO 5400 JI3,NM1-1
AI(J) I 1.
BI(J) I I4.
CI(J) I l.

5400

5450
C
5500

0

0 000

5600

5700
5950

183

DI(J) -H2*Q(I,J) - F(I—1,J) - F(I+1,J)
AI(NMl) - 1.
BI(NM1) I -4.
CI(NM1) - 0

DI(NM1) - IH2*Q(I,NM1) - F(I-1,NM1) — F(I+1,NM1) - F(I,N)

CALL TRI (2,NM1) ! tri-diagonal system.solver;
solves for F at fixed I

DO 5450 JI2,NM1

F(I,J) I WI(J)

F(I,J) I RFS*F(I,J) + C2*FOLD(I,J)
CONTINUE

IF (ITMXF.EQ.1) RETURN

Test the Stream Function for Convergence:

 

ERF I 0.

DO 5600 II2,MM1
DO 5600 JI2,NM1
DERF I F(I,J) + EPSILON

ERF - DMAXl(ERF,DABS((F(I,J)-FOLD(I,J))/DERF))

IF (ERF.LE.EPSILONF) GO TO 5950
CONTINUE

RETURN .
END

 

*‘k'k****i**************************************************************

 

00000000000000000000000000 00000

SUBROUTINE TRI (K,L)

 

A subroutine for solving a system of simultaneous linear algebraic
equations with a tri-diagonal coefficient matrix. It uses the
Thomas algorithim method which is a special adaptation of the

Gaussian elemination method in the form of a
The equations are numberd from K through L.
diagonal, and superdiagonal coefficients are
A, B, and C, respectively, each of dimension
of length (L) containing the right—hand side
The computed solution is the vector W(L).

recursion formula.
Their subdiagonal,
stored in the arrays
(L). D is a vector
of the linear system.

 

DESCRIPTION OF OTHER VARIABLES EXECLUSIVE TO

THIS ROUTINE

K, L I subscripts corresponding, respectively, to
the first and last equations to be solved
[usually KIl for a system of L equations]

BETA(L), GAMA(L) I vectors of the intermediate coefficients used
in the solution of the system as demanded by

the recursion formula

 

INCLUDE ’INCLUDE.FOR’

! include file

0(30

2100

0 0 0(30

2200

0

184

REAL*8 BETA(IDIM), GAMA(IDIM)
INTEGER K,L
INTEGER KP1,LMK

Computing intermediate arrays:

 

BETA(K) I B(K)
GAMA(X) I D(K)/BETA(K)

KPl - x + 1

DO 2100 IIKP1,L

BETA(I) - B(I) - A(I)*C(I-1)/BETA(I-1)
GAMA(I) - (D(I) - A(I)*GAMA(I-1))/BETA(I)

Computing the solution vector:

 

W(L) I GAMA(L)

LMK I L I R

DO 2200 JI1,LMK

I I L I J

W(I) I GAMA(I) I C(I)*W(I+1)/BETA(I)

RETURN
END

 

**********************************************************************

 

0 0()0(00()0(30 0(30030

0()0 0

9620

9630
9635

SUBROUTINE MORE (ILE,KINC,LUNITW)

 

Subroutine to compute the velocities and write data and calculated
results at any time step during the one full period before the
final time step (at convergence) into the output file.

 

INCLUDE 'INCLUDE.FOR' ! include file

INTEGER KINC, ILE, MSKIP
INTEGER LUNITW

PHI I PHIA(KINC)

Compute the Streamwise Velocity:

 

DO 9630 I-1,MM1

U(I,l) - o.

U(I,2) - (I25.*F(I,2)+48.*F(I,3)-36.*F(I,4)+16.*F(1,5)
I3.*F(I,6))/(12.*H)

Do 9620 JI3,N-2

U(I,J) - (-F(I,J+2)+8.*F(I,J+1)-8.*F(I,J-1)+F(I,J-2))/(12.*H)

JINMl

U(I,J) - ( 25.*F(I,J)-48.*F(I,J-1)+36.*F(I,J-2)-16.*F(I,J-3)
+3.*F(I,J-4))/(12.*H)

U(I,N) - 0.

Do 9635 I-M1,M2

U(I,1) - PHI

O

96

9E

000w

D ‘0
no

981

986

185

 

DO 9640 JI1,N

9640 U(M,J) I 6*Y(J) I 6*(Y(J)**2)
U(1,1) I U(l,NC)
U(l,N) I U(l,NC)

C Compute the Normal Velocity:

 

DO 9650 II1,M
9650 V(I,1) I 0.
DO 9660 JI2,N
V(1,J) I o.
V(2,J) - I(I25.*F(2,J)+48.*F(3,J)I36.*F(4,J)+16.*F(5,J)
& I3.*F(6,J))/(12.*H)
DO 9655 II3,M-2
9655 V(I,J) I I(-F(I+2,J)+8.*F(I+1,J)I8.*F(I-l,J)+F(I-2,J))/(12.*H)
V(MM1,J) I I(25.*F(MM1,J)I48.*F(MI2,J)+36.*F(MI3,J)I16.*F(MI4,J)
& +3.*F(MIS,J))/(12.*H)
9660 V(M,J) I 0.
C

C Continue Writing Data and Calculated Results into Output File:

 

WRITE (LUNITW,9859)
9859 FORMAT ('1',/)

C

WRITE (LUNITW,9860)
9860 FORMAT (2X,’OUTPUT FROM THE LAST PERIOD BEFORE EXITING',/)
C

TKK I KKA(KINC)*KT

RTP - TKK/(2.*PI)

WRITE (LUNITW,9861) KKA(KINC),ITERA(KINC),TA(KINC),TKK,KTP
9861 FORMAT (/,2X,'KK I ',16 ,/,

& 2X,'ITER I ',16 ,/,
& 2X,'TIME, T - ',F9.4,/,
& 2X,'TKK I ',F9.4,/,
& 2X,'KTP - ’,F9.4,/)

WRITE (LUNITW,9862)
9862 FORMAT ( I.

 

 

& 22x,' VELOCITY’,2X,’ VELOCITY',
& 4x,' STREAM FN',2X,' VORTICITY’,/,
& 2x,' x ' ,2x,' Y ' ,
a 4x,' U ',2x,' v ',
8 4X,’ F ',2X,' Q ',/)
C
I I 0
9863 IF (I.LT.M) THEN
IF (I.LE.100) THEN
I I I + 1
ELSE
MSKIP I 2*(INT(I/100)+1)
I I I + MSKIP
IF (I.GT.(100*MSKIP/2 + 1)) I I 100*INT(I/100) + 1
ENDIF
DO 9866 JI1,N
WRITE (LUNITW,9864) X(I),Y(J),U(I,J),V(I,J),F(I,J),Q(I,J)
9864 FORMAT (2(2X,F7.3),2X,2(2X,F10.6),2X,2(2X,F10.6))
9866 CONTINUE
WRITE (LUNITW,9867)
9867 FORMAT (2X,' ',
& l I)
GO TO 9863
ENDIF

9881 MSKIP I INT(MM1/50)
IF (MSKIP.LT.1) MSKIP I 1

0

0 00000000000000000000000000000000000000

186

RETURN
END

 

*******************t****‘k*********************************************

 

**

File Name: INCLUDE.FOR

 

A FORTRAN file to include the following FORTRAN statements which

are shared between the main program and all subroutines:

1. PARAMETER statements defining arrays dimensions IDIM, JDIM and

KDIM; these are the maximum allowable number of array elements

in the x and Y directions and in time, respectively.

2. DECLERATION statments for all global variables shared between
the main program and all subroutines.

3. Common statements representing common variables and arrays;
this will reserve the same memory location for these arrays.

To include this file in the FORTRAN code (main program and all
related subroutines) use the VAX/VMS FORTRAN statement:
INCLUDE 'INCLUDE.FOR'
The contents of this file will then be included in all FORTRAN
programs and subroutines; i.e., IDIM and JDIM need no longer be
be defined in a PARAMETER statement; variables need no longer be
be passed as a subroutine argument or in a COMMON statement
(i.e., common block) in the main programs or subroutines.

All variables and arrays in this file are defined in the main
program or in the subroutines where they occur.

 

Declaration of Variables; Double Precision for All Real Variables:

 

PARAMETER (IDIM-2001,JDIMI51) ! define space dimensions
of arrays

PARAMETER (IDIM10I(IDIMI1)/10)

PARAMETER (KDIMIIZO)

PARAMETER (KDIM8I8)

! define time dimension of arrays

REAL*8 RE, ST

REAL*8 RFS

REAL*8 va

REAL*8 H, H2, TH

REAL*8 X(IDIM), Y(JDIM), XR(IDIM)
REAL*8 F(IDIM,JDIM), Q(IDIM,JDIM)
REAL*8 FS(IDIM,JDIM), QS(IDIM,JDIM)
REAL*8 FOLD(IDIM,JDIM)

REAL*8 PI

REAL*8 KT

REAL*8 PHI

REAL*8 TA(KDIM), PHIA(KDIM)

REAL*8 TKK, KTP

REAL*8 C1, C2

REAL*8 G1, 62

REAL*8 EPSILON, EPSILONF, EPSILONQ
REAL*8 ERF, ERQ

000

0000

187

 

REAL*8 U(IDIM,JDIM), V(IDIM,JDIM)
REAL*8 A(IDIM), D(IDIM), C(IDIM), D(IDIM), W(IDIM)
REAL*8 A1(IDIM), A2(IDIM), A3(IDIM), BB(IDIM), WW(IDIM)
REAL*8 AI(IDIM), BI(IDIM), CI(IDIM), DI(IDIM), WI(IDIM)
INTEGER M, N, 1041, NM1, NC, M1, M2
INTEGER ITERA(KDIM), KKA(KDIM)
INTEGER ITMAx, ITER, ITMXF, ITERF, ITMXQ, ITERQ
EQUIVALENCE (A,A1,AI), (B,A2,BI), (C,A3,CI), (D,BB,DI),
& (W,WW,WI)
Common Blocks:
COMMON / NAME010 / RE
COMMON / NAME011 / ST
COMMON / NAMEozo / M, N
COMMON / NAME021 / MM1, NM1
COMMON / NAME023 / NC
COMMON / NAME025 / M1, M2
COMMON / NAME028 / M81, M52, NS
COMMON / NAMEO30 / H
COMMON / NAME031 / H2
COMMON / NAME032 / TH
COMMON / NAME039 / PI
COMMON / NAME040 / RFS
COMMON / NAME041 / BKV
COMMON / NAME045 / C1, C2
COMMON / NAME051 / KT
COMMON / NAME053 / 61, G2
COMMON / NAME054 / PHI
COMMON / NAME056 / ITERA, KKA
COMMON / NAME058 / TA, PHIA
COMMON / NAME058 / TKK, KTP
COMMON / NAME060 / ITMAx, ITER
COMMON / NAME061 / ITMXF, ITERF
COMMON / NAME062 / ITMxo, ITERQ
COMMON / NAME070 / EPSILON
COMMON / NAME071 / EPSILONF
COMMON / NAME072 / EPSILONQ
COMMON / NAME075 / ERF
COMMON / NAME076 / ERQ
COMMON / NAME080 / x, Y, XR
COMMON / NAME091 / F, Q
COMMON / NAME093 / Fs, QS
COMMON / NAME095 / FOLD
COMMON / NAMEloo / A, B, C, D, W
COMMON / NAME111 / U, v

 

******************i**********t****************************************

 

BIBLIOGRAPHY

BIBLIOGRAPHY

Abarbanel, S., Bennett, 8., Brandt, A., and Gillis, 1., 1970, “Velocity Proﬁles Of Flow at
Low Reynolds Number”, ASME Journal of Applied Mechanics, Vol. 37, pp. 2-4.

Abdulla, N. N., 1987, “Laminar Flow Separation in a Constricted Channel”, Ph.D.
Dissertation, Department Of Mechanical Engineering, Michigan State University, East Lansing, MI.

Allen, D., and Southwell, R. V., 1955, “Relaxation Methods Applied to Determine the
Motion, in Two Dimensions, Of a Viscous Fluid Past 11 Fixed Cylinder”, Quarterly Journal of
Mechanics and Applied Mathematics, Vol. 8, pp. 129-145.

Ames, W. F., 1977, Numerical Methods for Partial Diﬂ’erential Equations, 2nd Ed., Academic
Press, New York.

Anderson, C. R., 1989, “Observations on Vorticity Creation Boundary Conditions”,
Mathematics Aspects of Vortex Dynamics, E. C. Caﬂisch, ed., pp. 144-159, SIAM.

Anderson, D. A., Tannehill, J. C., and Fletcher, R. H., 1984, Computational Fluid Mechanics
and Heat Transfer, Hemisphere Publishing, New York.

Aref, H., 1986, “The Numerical Experiment in Fluid Mechanics”, Journal of Fluid
Mechanics, Vol. 173, pp. 15-41.

Armaly, B. F., Durst, R, Pereira, J. C. F., and Schdnung, B., 1983, “Experimental and
Theoretical Investigation of Backward-Facing Step Flow”, Journal of Fluid Mechanics, Vol. 127,
pp. 473-496.

Arpaci, V. 8., and Larsen, P. S., 1984, Convection Heat Transfer, Prentice-Hall, New Jersey.

Aziz, K., and Hellums, J. D., 1967, “Numerical Solution of the Three-Dimensional Equations
of Motion for Laminar Natural Convection”, Physics of Fluids, Vol. 10, pp. 314 -324.

Backus, J. W., 1989, Private Communications, Department of Mechanical Engineering,
Michigan State University, East Lansing, MI.

188

189

Barakat, H. 2., and Clark, J. A., 1966, “Analytical and Experimental Study of the Transient
laminar Natural Convection Flows in Partially Filled Liquid Containers”, Proceedings of the Third
International Heat Transfer Conference, Vol. 2, pp. 152- 162, AIChE, New York.

Benjamin, A. S., and Denny, V. B., 1979, “On the Convergence of Numerical Solutions for 2-D
Flows in a Cavity at Large Re”, Journal of Computational Physics, Vol. 33, pp. 340-358.

Birkhoff, G., and Varga, R. S., 1959, “Implicit Alternating Direction Methods”, Transactions
of the American Mathematical Society, Vol. 92, pp. 13-24.

Birkhoff, 6., Varga, R. S., and Young, D., 1962, “Alternating Direction Implicit Methods”,
Advanced in Computers, F. L. Alt et al., ed, Vol. 3, pp. 189-237, Academic Press, New York.

Bodoia, J. R., 1959, "Ihc Finite Difference Analysis of Conﬁned Viscous Flows”, Ph.D.
Thesis, Department of Mechanical Engineering, Carnegie Institute of Technology, Pittsburgh, PA.

Bodoia, J. R., and Osterle, J. F., 1961, “Finite Difference Analysis of Plane Poiseuille and
Couette Flow Developments”, Applied Scientiﬁc Research, See. A, Vol. 10, pp. 265-276.

Bontoux, P., Gilly, B., and Roux, B., 1980, “Analysis of the Effect of Boundary Conditions
on Numerical Stability of Solutions of Navier-Stokes Equations", Journal of Computational Physics,
VOL 36, pp. 417-427.

Boussinesq, J., 1891, “Sur la martiere dont les vitesses, dans un tube cylindrique dc section
circulaire, evasé 8 son entree, se distribuent depuis cette entree jusqu’aux endroits 00 se trouve établi
un régime uniftrme”, Comptes Rendus Hebdomadaires des Seances de l'Académie des Sciences,
Vol. 113, pp. 9-15; also pp. 49-51.

Bozeman, J. D., and Dalton, C., 1973, “Numerical Study of Viscous Flow in a Cavity”,
Journal of Computational Physics, Vol. 12, pp. 348-363.

Brandt, A., and Gillis, J., 1966, “Magnetohydrodynamics Flow in the Inlet Region of a
Straight Channel”, Physics of Fluids, Vol. 9, pp. 690-699.

Briley, W. R., 1970, “A Numerical Study of Laminar Separation Bubbles Using the Navier-
Stokes Equations”, Report 1110614-1, United Aircraft Research laboratories, East Hartford, CT.

Briley, W. R., 1971, “A Numerical Study of laminar Separation Bubbles Using the Navier-
Stokes Equations”, Journal of Fluid Mechanics, Vol. 47, pp. 713-736.

Burggraf, O. R., 1966, “Analytical and Numerical Studies of the Structure of Steady Separated
Flows”, Journal of Fluid Mechanics, Vol. 24, pp. 113-151.

Campbell, W. D., and Slattery, J. C., 1963, “Flow in the Entrance of 8 Tube”, ASME Journal
of Basic Engineering, Vol. 85, pp. 41-46.

Caughey, D. A., and Iyer, R. K., 1989, “Diagonal Implicit Multigrid Calculation of Inlet
Flowﬁelds”, AIAA Journal, Vol. 27, pp. 110-112.

Chow, C.-Y., 1979, An Introduction to Computational Fluid Mechanics, John Wiley, New
York.

Christiansen, E. B., and Lemmon, H. B., 1965, “Enhance Region Flow”, AIChE Journal, Vol.
ll,pp. 995-999.

190

Collins, M., and Schowalter, W. R., 1962, “laminar Flow in the Inlet Region of a Straight
Channel”, Physics of Fluids, Vol. 5, pp. 1122-1124.

Collins, M., and Schowalter, W. R., 1963, “Behavior of Non-Newtonian Fluids in the Inlet
Region of a Channel”, AIChE Journal, Vol. 9, pp. 98-102.

Daube, 0., and Ta Phoc Lee, 1979, “A Mixed Compact Hermitian Method for the Numerical
Study of Unsteady Viscous Flow Around an Oscillating Airfoil”, Proceeding of the Third GAMM
Conference on Numerical Methods in Fluid Mechanics, E. H. Hirschel, ed., Notes on Numerical
Fluid Mechanics, Vol. 2, pp. 56-65, Vieweg and Sohn, Braunschweig.

Davis, R. T., 1967, “Laminar Incompressible Flow Past a Semi-Infinite Flat Plate”, Journal of
Fluid Mechanics, Vol. 27, pp. 691-704.

Dennis, S. C. R., and Chang, G.-Z., 1969, “Numerical Integration of the Navier-Stokes
Equations for Steady Two-Dimensional Flow”, Physics of Fluids, VOL 12 (Supplement 11), pp. 11-
88-11-93.

Dennis, S. C. R., and Smith, F. '11, 1980, “Steady Flow Through a Channel with a Symmetrical
Constriction in the Form of a Step”, Proceedings of the Royal Society of London, Series. A, Vol.
372, pp. 393-414.

Deshpande, M. D., Giddens, D. P., and Mabon, R. F., 1976, “Steady laminar Flow Through
Modelled Vascular Stcnoses", Journal of Biomechanics, Vol. 9., pp. 165-174.

Douglas, J., 1955, “On the Numerical Integration of 32u/3x2 + azu/ayz . au/at by Implicit
Methods", Journal ofthe Society oflndustrial and Applied Mathematics, Vol. 3, pp. 42455.

Douglas, J., and Rachford, H. H., 1956, “On the Numerical Solution of Heat Conduction
Problems in Two and Three Space Variables”, Transactions of the American Mathematical Society,
Vol. 82, pp. 421-439.

Douglas, J., and Gunn, J. B., 1964, “A General Formulation of Alternating Direction
Methods: Part I. Parabolic and Hyperbolic Problems”, Numerische Mathematik, Vol. 6, pp. 428-453.

Durst, F., and Pereira, J. C. F., 1988, “lime-Dependent Lamina Backward-Facing Step Flow
in a Two-Dimensional Duct", ASME Journal of Fluids Engineering, Vol. 110, pp. 289-296.

Duck, P. W., 1982, “Oscillatory flow inside a square cavity”, Journal of Fluid Mechanics,
Vol. 122, pp. 215-234.

Fairweather, G., and Mitchell, A. R., 1967, “A New Computational Procedure for A.D.I.
Methods”, SIAM Journal of Numerical Analysis, Vol. 4, pp. 163-170.

Fargie, D., and Martin, B. W., 1971, “Developing laminar Flow in a Pipe of Circular Cross-
Section", Proceedings of the Royal Society of London, Series. A, Vol. 321, pp. 461-476.

Feliss, N. A., Smith, M. C., and Potter, M. C., 1977, “An Experimental Investigation of
Inoompressible Channel Flow Near Transition”, ASME Journal (1' Fluids Engineering, Vol. 99, pp.
693-698.

Fasythe, G. E., and Wasow, W. R., 1960, Finite-Diference Methods for Partial Diﬂ'erential
Equations, John Wiley, New York.

191

Friedman, M., 1970, “Flow in a Circular Pipe with Recessed Walls”, ASME Journal of
Applied Mechanics, Vol. 37, pp. 5-8.

Friedman, M., 1972, “Laminar Flow in a Channel with a Step”, Journal of Engineering
Mathematics, Vol. 6, pp. 285-290.

Friedmann, M., Gillis, J., and Liron, N., 1968, “lamina- Flow in a Pipe at Low and Moderate
Reynolds Numbers", Applied Scientiﬁc Research, Vol. 19, pp. 426- 438.

Frankel, S. P., 1950, “Convergence Rates of Iterative Treatments of Partial Differential
Equatiom”, Mathematical Tables and Other Aids to Computation, Vol. 4, pp. 65-75.

Fromm, J. B., 1963, “A Method frr Computing Nonsteady, Incompressible, Viscous Fluid
Flows”,ReportLA-2910,Los Alamos Scientiﬁc laboratory, Los Alamos, NM.

Fromm, J., 1964, "The Time Dependent Flow of an Incompressible Viscous Fluid”,
Fundamental Methods in Hydrodynamics, B. Alder et al., ed., Methods in Computational Physics,
Vol. 3 , pp. 345-382, Academic Press, New York.

Fromm, J. B., 1972, “Numerical Solution of the Navies Stokes Equations at High Reynolds
Numbers and the Problem of Discretization of Convective Derivatives”, Numerical Methods in Fluid
Dynamics, J. J. Smoldern, ed., AGARD Lecnlre Series No. 48, Advisory Group for Aerospace
Research and Development, North Atlantic Treaty Organization, pp. 6.1-6.44.

Fromm, J. B., and Harlow, F. H., 1963, “Numerical Solution of the Problem of Vortex
Street Development”, Physics of Fluids, Vol. 6, p. 975-982.

Ghia, U., Ghia, K. N., and Shin, C. T., 1982, “High-Re Solutions for Incompressible Flow
Using the Navier-Stokes Equations and a Multigrid Method”, Journal of Computational Physics,
Vol. 48, pp. 387-411.

Gillis,J.,andBrandt,A.,l964,“1heNurnericallntegrationoftheEquationsofMotionofa
Viscous Fluid”, Report SR-1,AFEOAR 63-73,AirForce EmopeanOffrceof Aerospace Research.

Goldstein, S., 1965, Modern Developments in Fluid Dynamics, Vol. 1, Dover Publications,
New York.

Gasman, A. D., Pun, W. M., Runchal, A. K., Spalding, D. B., and Wolfshtein, M., 1969, Heat
and Mass Transfer in Recirculating Flows, Academic Press, New York.

Greenspan, D., 1965, Introductory Numerical Analysis of Elliptic Boundary Value Problems,
Harper and Row, New York.

Greenspan, D., 1967, “Numerical Studies of Two Dimensional Steady State Navies-Stokes
Equations for Arbitrary Reynolds Number", Report No. 9, Computer Sciences Department,
University of Wisconsin, Madison, WI.

Greenspan, D., 1968, “Numerical Studies of Viscous, Incompressible Flow for Arbitrary
Reynolds Number”, Report No. 11, Computer Sciences Department, University of Wisconsin,
Madison, WI.

Greenspan, D. J., 1969a. “Numerical Studies of Steady, Viscous, Incompressible Flow in a
Channel with a Step”, Journal of Engineering Mathematics, Vol. 3, pp. 21-28.

192

Greenspan, D. J., 1969b, “Numerical Studies of Prototype Cavity Flow Problem”, The
Computer Journal, Vol. 12, pp. 88-93.

Greenspan, D. J., 1974, Discrete Numerical Methods in Physics and Engineering, Academic
Press, New York.

Gupta, M. M., and Manohar, R. P., 1979, “Bormdary Approximations and Accuracy in Viscous
Flow Computations”, Journal of Computational Physics, Vol. 31, pp. 265-288.

Harlow, F. H., and Fromm, J. E., 1965, “Computer Experiments in Fluid Dynamics”,
Scientiﬁc America, Vol. 212, pp. 104-110.

Hanbeck, R. W., 1964, “Laminar Flow in the Entrance Region of a Pipe”, Applied Scientific
Research, See. A, Vol. 13, pp. 224-232.

Hwang, C.-L., and Fan, L.-T., 1963, “A Finite Difference Analysis of Laminar Magneto-
Hydrodynamic Flow in the Entrance Region of a Flat Rectangular Duct”, Applied Scientﬁc
Research, Sec. B, Vol. 10, pp. 329-343.

Jenson, V. G., 1959, “Viscous Flow Round a Sphere at Low Reynolds Numbers (<40)”,
Proceedings of the Royal Society of London, Series A, Vol. 249, pp. 346-366.

Kawaguti, M., 1961, “Numerical Solution of the Navier-Stokes Equations for the Flow in a
Two-Dimensional Cavity”, Journal of the Physical Society of Japan, Vol. 16, pp. 2307-2315.

Khosla, P. K., and Rubin, S. G., 1974, “A Diagonally Dominant Second-Order Accruate
Implicit Scheme", Computers and Fluids, Vol. 2, pp. 207-209.

Krrmae, B., 1985, “Computational Fluid Dynamics: Its Present Status and Future Direction”,
Computers and Fluids, Vol. 13, pp. 239-269.

Km, A., and Yajnik, K. S., 1980, “Internal Separated Flows at Large Reynolrk Number”,
Journal of Fluid Mechanics, Vol. 97, pp. 27-51.

Kuskova, T. V., 1967, “Difference Method for Calculating the Flow of a Viscous Fluid” (in
Russian), Sbornik Robot Vychislitel’nogo Tsentra Moskovshogo Gosudarstvennogo Universiteta, Vol.
7, pp. 16-27; also 1969. APPlied Mechanics Reviews (Review No. 380), Vol. 22, p. 61.

Kwon, 0., Pletcher, R. H., and Delaney, R. A., 1988, “Solution Procedure for Unsteady Two-
Dimensional Boundary Layers”, ASME Journal of Fluids Engineering, Vol. 110, pp. 69-75.

Lsnghmr, H. 1... 1942, “Steady Flow in the Transition Length of a Straight Tube”, ASME
Journal of Applied Mechanics, Vol. 9, pp. ASS-A58.

Leal, L. G., and Acrivos, A., 1969, “Structme of Steady Closed Streamline Flows within a
Bormdary Layer”, Physics of Fluids, Vol. 12 (Supplement 11), pp. II-105 -II-l l3.

Lew, H. 8., and Fung, Y. C., 1970, “Entry Flow into Blood Vessels at Arbitrary Reynolds
Number", Journal of Biomechanics, Vol. 3, pp. 23-38.

Li, L. C., and Ludford, G. S. S., 1980, “The Overshoot in Entry Flow”, Archives of
Mechanics, Vol. 32, pp. 741-746.

193

Lilly, D. IL, 1965, “On the Computational Stability of Numerical Solutions of Time-
Dependent Non-Linear Geophysical Fluid Dynamics Problems", Monthly Weather Review, Vol. 93,
pp. 11-26.

erdgren, T. S., Atabek, B. H., and Chang, C. C., 1961, “Transient Magnetohydrodynarnic
Duct Flow”, Physics of Fluids, Vol. 4, pp. 1m6-1011.

erdgren, T. S., Sparrow, E. M., and Starr, J. B.,1964, “Pressure Drop Due to the Entrance
Region in Ducts of Arbitrary Cross Section”, ASME Journal of Basic Engineering, Vol. 86, pp. 620-
626.

Lynch, R. E., and Rice, J. R., 1968, “Convergence Rates of ADI Methods with Smooth Initial
Error”, Mathematics of Computation, Vol. 22, pp. 311-335.

Mallinson, G. D., and de Vahl Davis, 6., 1973, "The Method of False Transient for the
Solution of Coupled Elliptic Equations”, Journal of Computational Physics, Vol. 12, pp. 435-461.

McDonald, J. W., Denny, V. E., and Mills, A. F., 1972, “Numerical Solutions of the Navier-
Stokes Equations in Inlet Regions", ASME Journal of Applied Mechanics, Vol. 39, pp. 873-878.

Mei, R., and Plotkin, A., 1986, “A Finite-Difference Scheme for the Solution of the Steady
Navier-Stokes Equations", Computers and Fluids, Vol. 14, pp. 239-251.

Mills, R. D., 1965, “Numerical Solutions of the Viscous Flow Equations for a Class of
Closed Flows”, Journal of the Royal Aeronautical Society, Vol. 69, pp. 714-718; also correction p.
880.

Mitchell, A. R., 1969, Computational Methods in Pam'al Dﬁ‘erential Equations, John Wiley,
New York.

Mitchell,A.R.,andFairweather,G., 1964,“Im}rovedFormsoftheAlternatiugDirection
Methods of Douglas, Peaceman, and Rachford for Solving Parabolic and Elliptic Equations”,
Numerische Mathematik, Vol. 6, pp. 285-292.

Mohanty, A. K., and Asthana, S. B. L.,1978, “laminar Flow in the Entrance Region of a
Smooth Pipe, Journal of Fluid Mechanics, Vol. 90, pp. 433-447.

Maihara, H., and Cheng, R. T., 1973, “Numerical Solution of the Viscom Flow in the
Enhance Region of Parallel Plates”, Journal of Computational Physics, Vol.11, pp. 550-572.

Morris, D. J., 1975, “Solution of the Incompressible Driven Cavity Problem by the
Alternating-Direction Implicit Method”, Numerical Studies of Incompressible Viscous Flow in a
Driven Cavity, NASA SP-378, National Aeronautics and Space Administration, washinsm D.C.,
pp. 47-60.

Nallasamy, M., 1986, “Numerical Solution of the Separating Flow Due to an Obstruction”,
Computers and Fluids, Vol. 14, pp. 59-68.

Orszag, S. A., and Israeli, M., 1974, “Numerical Simulation of Viscous lncornpressible
Flows”, Annual Review of Fluid Mechanics, Vol. 6, pp. 281-318.

Pan, F., and Acrivos, A., 1967, “Steady Flows in Rectangula' Cavities”, Journal of Fluid
Mechanics, Vol. 28, pp. 643-655.

194

Pao, Y.-I-1., and Daugherty, R. J ., 1969, “Tune-Dependent Viscous Incompessible Flow Past a
Finite Flat Plate”, Report D1-82-0822, Boeing Scientiﬁc Research Iaborataies, Seattle, WA.

Paris, J., and Whitaker, S., 1965, “Conﬁned Wakes' A Numerical Solution of the Navier-
Stokes Equations”, AIChE Journal, Vol. 11, pp. 1033-1041.

Peaceman, D. W., and Rachford, H. H., Jr., 1955, "The Numerical Solution of Parabolic and
Elliptic Differential Equations”, Journal of the Society of Industrial and Applied Mathematics, Vol.
3, pp. 28-41.

Pearson, C. E., 1964, “A Coruputational Method for Time-Dependent Two-Dimmsional
Incompressible Viscous Flow Problems”, Report SRRC-RR-64-17, Sperry Rand Research Center,

Subdtu'y, MA.

Pearson, C. E., 1965a, “A Computational Method Fa Viscous Flow Problems", Journal of
Fluid Mechanics, Vol. 21, pp. 611-622.

Pearson, C. E., 1965b, “Numerical Solutions for the Time-Dependent Viscous Flow Between
Two RotatingCoaxialDisks”, Journal of Fluid Mechanics, Vol. 21,pp. 623-633.

Peyret, R., and Taylor, T. D., 1983, Computational Methods for Fluid Flow, Springer-Verlag,
New York.

Peyret, R., 1971, “Sun 1e calcul numerique des 6coulements stationnaires d’un ﬂuide visqueux
incompressible" Comptes Rendus Hebdomadaires des Seances de l'Academie des Sciences, Série A,
Vol. 272, pp. 1274-1277.

Prandtl, I... and Tietjens, O. G., 1934, Applied Hydro- and Aeromechanics, Dover
Publications, New York.

Richtmyer, R. D., lid Morton, K. W., 1967, Dw'erence Methods for Initial-Value Problems.
2nd Ed., Intaacience Publishers, New York.

Roache, P. J, 1975, “The LAD, NOS and Split NOS Methods for the Steady-State Navier-
Stokes Equations", Computers and Fluids, Vol. 3, pp. 179-195.

Roache, P. J, 1982, Computational Fluid Dynamics, Revised Ed., Hermosa Publishes,
Albuquerque, NM.

Roidt, M, and Cess, R. D., 1962, “An Approximate Analysis of laminar
Magnetohymodynamic Flow in the Entrance Region of a Flat Duct”, ASME Journal ofth
Mechanics, Vol. 29, pp. 171-176.

Rubin, S. G., and Harris, J. E., 1975 “Incompressible Viscous Flow in a Driven Cavity”,
Numerical Studies of Incompressible Viscous Flow in a Driven Cavity, NASA SP-378, National
AeronauticsandSpace Administration, Washington, D.C., pp. 1-6.

Rubin, S. G., and Khosla, P. K., 1977, “Polynomial Interpolation Methods for Viscous Flow
Calculations”, Journal of Computational Physics, Vol. 24, pp. 217-244.

Runchal, A. K., Spalding, D. B., and Wolfshtein, M., 1969, “Numerical Solution of Elliptic
Equations for Transport of Vorticity, Heat, and Matter in Two-Dimensional Flow”, Physics (1‘
Fluids, Vol. 12 (Supplement 11), pp. II-21-II-28.

195

Samuels, M. R., and Chmchill, S. W., 1967, “Stability of a Fluid in a Rectangular Region
Heated from Below”, AIChE Journal, Vol. 13, pp. 77-85.

Schiller, L., 1922, “Die Entwicklrmg der laminaren Geschwindigkeitsverteilung und ilue
Bedeutung fur Zhhigkeitsmessungen” Zeitschriﬂ ﬁlr angewandte Mathematik und Mechanih, Vol. 2,
pp. 96-106.

Schlichting, H., 1934, “Lmninare Kanaleinlaufstrolnung", Zeitschrift ﬁZr angewandte
Mathematik und Mechanik, Vol. 14, pp. 368-373.

Schlichting, H., 1979, Boundary-Layer Theory, 7th Ed., McGraw-Hill, New York.

Schmidt, F. W., and Zeldin, B., 1969, "Laminar Flows in Inlet Sections of Tubes and Ducts”,
AIChE Journal, Vol. 15, pp. 612-614.

Son,J.S.,andHanratty,T.J.,1969,“NtunericaISolution£tn'theFlowArormdaCylinderat
Reynolds Numbers of 40, 200 and 500”, Journal of Fluid Mechanics, Vol. 35, pp.369-386.

Smith, F. T., 1977, “Upstream Interaction in Channel Flows”, Journal of Fluid Mechanics,
Vol. 79, pp.631-655.

Smith, F. T., 1979, "The Separating Flow Through a Severely Constricted Syrrunetric Tube”,
Journal of Fluid Mechanics, Vol. 90, pp.725-754.

Sparrow, E. M., Anderson, C. B., 1977, “Effect of Upstream Flow Processes on
Hydrodynamic Development in a Duct”, ASME Journal of Fluids Engineering, Vol. 99, pp. 556-560.

Sparrow, E. M., Hixon, C. W., and Shavit, G., 1967, “Experiments on Laminar Flow
Development in Rectangular Ducts”, ASME Journal (1 Basic Engineering, Vol. 89, pp. 116124.

Sparrow, E. M., Lin, S. H., and Lundgren, T. S., 1964, “Flow Development in the
Hydrodynamic Enhance Region of Tubes and Ducts”, Physics of Fluids, Vol. 7, pp. 333-347.

Ta Phoc Loc, and Daube, 0., 1981, “Higher Order Accurate Numerical Solution of Unsteady
Viscous Flow Generated by a Transversely Oscillating Ellipt'u: Cylinder”, Vortex Flows, W. L.
Swift et al., ed., pp. 155-171, ASME, New York.

Taylor, P. J., 1970, “The Linear Stability of No-Slip Boundary Conditiom in the Numerical
Solution of Nonsteady Fluid Flows”, Journal of Computational Physics, Vol. 6, pp. 268-387.

Taylor, T. D., and Ndefo, E., 1971, “Computation of Viscous Flow in a Charurel by the
Method of Splitting”, Proceedings of the Second International Conference on Numerical Methods in
Fluid Dynamics, M. Holt et al., ed., Lecture Notes in Physics, Vol. 8, pp. 356-364, Stringer-
Verlag, Berlin.

Taylor, T. D., and Ndefo, B., and Masson, B. S., 1972, “A study of Numerical Methods for
Solving Viscousandlnviscid Flow Problems”, Journal of Computational Physics, Vol. 9, pp. 99-119.

Thom, A., 1933, “The Flow Past Circula' Cylinders at Low Speeds”, Proceedings of the
Royal Society of London, Series. A, Vol. 141, pp. 651-669.

Thom, A., and Apelt. C. J., 1961, Field Computation in Engineering and Physics, Van
Nostrand, London.

196

Thompson, J. F., 1969, “Optimized Acceleration of Convergence of an Implicit Numerical
Solution of the Time-Dependent Navier-Stokes Equations”, AlAA Journal, Vol. 7, pp. 2186-2188.

Torrance, K. E., 1968, “Comparison of Entire-Difference Computations of Natmal
Convection”, Journal of Research of the National Bureau of Standards, Vol. 72B, pp. 281- 301.

Tuann, S.-Y., and M. D. Olson, 1978, “Review of Computing Methods for Recirculating
Flows”, Journal of Computational Physics, Vol. 29, pp. 1-19.

Van Dyke, M., 1970, “Entry Flow in a Cinnuel", Journal of Fluid Mechanics, Vol. 44, pp.
813-823.

Varga, R. S., 1962, Matrix Iterative Analysis, Prentice-Hall, New Jersey.

Vrentas, J. S., Duds, J. L., and Bargeron, K. G., 1966, “Effect of Axial Diffusion of Va'ticity
on Flow Development in Circular Conduits: Part 1. Numerical Solutions”, AIChE Journal, Vol. 12,
pp. 837-844.

Wachspress, E. L., 1966, Iterative Solution of Elliptic Systems, Prentice-Hall, New Jersey.

Wang, Y. L., and Longwell, P. A., 1964, “Laminar Flow in the Inlet Section of Parallel
Plates”, AIChE Journal, Vol. 10, pp. 323-329.

Wiginton, C. L., and Dalton, C., 1970, “Incompresm’ble laminar Flow in the Entrance Region
of a Rectangular Duct”, ASME Journal of Applied Mechanics, Vol. 37, pp. 854-856.

Wiginton, C. L., and Wendt, R. L., 1969, “Flow in the Entrance Region of Ducts”, Physics of
Fluids, Vol. 12, pp. 465-466.

Wilkes, J. 0., and Chmchill, S. W., 1966, “The Finite Difference Computation of Natural
Convection in a Rectangular Enclosrue”, AIChE Journal, Vol. 12, pp. 161-166.

Wilson, 8., 1969, “The Development of Poiseuille Flow”, Journal of Fluid Mechanics, Vol.
38, pp. 793-806.

Wilson, 8. D. R., 1971, “Entry Flow in a Channel. Part 2”, Journal of Fluid Mechanics, VOL
46, pp. 787-799.

Woods, L. C., 1954, “A Note on the Numerical Solution of Foruth Order Differential
Equations”, The Aeronautical Quarterly, Vol. 5, pp. 176-184.

Yanenko, N. N., 1971, The Method of Fractional Steps: The Solution of Problems of
Mathematical Physics in Several Variables, Springer-Verlag, New Yak.

Young, D., 1954, “Iterative Methods for Solving Partial Difference Equations of Elliptic
Type”, Transactions of the American Mathematical Society, Vol. 76, pp. 92-111.