HYDRODYNAMIC LIFT IN SEDIMENT
TRANSPORT
By
BARRY ARDEN BENEDICT
A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF
THE UNIVERSITY OF FLORIDA
IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE
DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1968
ACKNOWLEDGMENTS
The author wishes to express his sincere appreciation to
Dr. B. A. Christensen, Dr. J. H. Schmertmann, Dr. E. A. Farber,
and Dr. T. O. Moore for their service on his supervisory committee
and their interest in his work. The author is especially indebted
to the committee chairman, Dr. B. A. Christensen, whose guidance,
encouragement, and personal interest and enthusiasm were truly
valuable.
Appreciation is extended to the National Science Foundation,
under whose traineeship the author has been working.
The author also wishes to express his gratitude to the
University of Florida Computing Center for the computing facilities,
services, and aid extended by the Center.
Finally, the author wishes to extend thanks for the aid and
encouragement provided by his wife, who has persevered through many
trying times.
TABLE OF CONTENTS
Page
ACKNOWLEDGMENTS . . . . . . . . . . .... ii
LIST OF TABLES . . . . . . . . . . . . vi
LIST OF FIGURES . .. ......... . . . . . .. vii
LIST OF SYMBOLS . . . . ... .. . . . . .. x
ABSTRACT ...... .. ........... .. ........ xiii
CHAPTER
I INTRODUCTION. .. . . . . . . . . ..1
II RELATED BACKGROUND ON LIFT FORCE STUDIES . . . . 4
2.1 Historical interest . . . . . . . 4
2.2 Neglect of lift force .. . . . . 5
2.3 Indications of significance of lift . . . 6
2.4 Use of potential theory in lift studies on
single bodies . . . . . . . . 7
2.4.1 General . . . . . . . 7
2.4.2 Applicability of potential flow . . 8
2.4.3 Work of Kutta and Joukowsky . . . 9
2.4.4 Fuhrmann's work and other studies . 11
2.4.5 Jeffreys' analysis . . . . . 12
2.4.6 Flow around a single sphere .. . . 14
2.5 Multiparticle studies .. . . . . . 19
2.5.1 Einstein and El Samni . . . . 19
2.5.2 Chepil . . . . . . . . 20
2.5.3 Chao and Sandborn . . . . . 22
2.6 Shapes of bodies studied . . . . . .. 22
2.7 Relation to work of dissertation . . . .. 23
III LOGARITHMIC VELOCITY DISTRIBUTION. . . . . . 25
3.1 Development of logarithmic velocity
distribution . . . . .. ... .. 25
3.2 Special problems of present expressions . 26
TABLE OF CONTENTS (Continued)
CHAPTER Page
3.2.1 Theoretical bed . . . . . ... 27
3.2.2 Conditions near the bed . . . ... 29
3.3 Use of proposed adjusted velocity distribution 31
3.3.1 Comparison of distributions at wall . 31
3.3.2 Comparison of distribution with
increasing y . . . . ..... 32
3.4 Determination of k and vf from experimental
data . . .. . . . . . . 33
3.5 Effect of sidewalls . . . . . ... 35
IV TWODIMENSIONAL WORK: EARLIER RESULTS . . . ... 41
4.1 Shapes studied . . . . . . .. 41
4.2 General methodsvelocity and pressure results 41
4.3 Lift integration . . . . . . ... 44
4.4 Calculation of predicted lift by use of
measured logarithmic velocity profile . . 47
4.5 Application to experimental results . . .. 49
V THREEDIMENSIONAL NUMERICAL SOLUTIONS . . . ... 52
5.1 General . . . . . . . .... . 52
5.2 Problem formulation . . . . . .... 53
5.2.1 Choice of solution method ...... 53
5.2.2, Depth of flow space ... . . . . 54
5.2.3 Boundary conditions . . . . ... 56
5.3 General finite differences approach . . .. 61
5.4 Finite differences equations: interior space 64
5.4.1 General lattice point ... .. . . . 65
5.4.2 Object point on planar noflow boundary 68
5.4.3 Object point on "foldedsymmetry"
boundary . . . . . . ... 70
5.4.4 Adjacent point on hemispherical surface 70
5.4.5 Graded lattice . . . . . ... 72
5.5 Finite difference equations: hemispherical
boundary . . . . . . . .. 78
5.5.1 General . . . . . . . . 78
5.5.2 Xdirection derivative . . . .. 79
TABLE OF CONTENTS (Continued)
CHAPTER Page
5.5.3 Ydirection derivative . . . .. 81
5.5.4 Zdirection derivative . . . ... 87
5.5.5 Adjacent zpoints and ypoints subjected
to normal derivative condition . 88
5.5.6 Final boundary formulation . . ... 91
5.5.7 Singular points . . . . . . .. 93
5.6 Velocity and lift calculations . . . .. 94
5.7 Implementation of solution . .. . . .. 99
VI RESULTS AND COMPARISONS . . . . . . .. .. 101
6.1 General . . . . . . . .... ... 101
6.2 Numerical results for Chepil arrangement . .. 102
6.3 Comparisons with Chepil's observations . . 110
6.3.1 Details of Chepil's work . . . .. 110
6.3.2 Comparison of lift forces . . .. . 116
6.4 Numerical results for closely packed
hemispheres . . . . . . . . 125
6.5 Comparison with EinsteinEl Samni observations 138
6.5.1 Physical details of experiments . . 138
6.5.2 Values of lift for hemisphere bed . .. 141
6.5.3 Lift on a gravel bed. . . . . 142
VII CONCLUSIONS AND FUTURE WORK . . . . . .. 146
APPENDIX . . . . . . . . . .. . . 149
NOTES ON FORTRAN IV COMPUTER PROGRAM . . . .. . 150
NOTES ON EQUIVALENT GRAIN SIZE . . . . . .. 200
REFERENCES . . . . . . . . . . . . . 201
BIOGRAPHICAL SKETCH . . . ... . . .. . . . 207
LIST OF TABLES
Table Page
1 Comparison of Proposed and Former Distributions ... .. 32
2 Chepil's Experimental Data . . . . . . ... 112
3 CL for Chepil's Work .................. 114
4 RoughnessGrain Size Ratios . . . . . . ... 115
5 Comparison of Theoretical with Chepil's Work ...... 123
LIST OF FIGURES
Figure Page
1 Flow around a Joukowsky profile . . . . ... 10
2 Calculated and measured pressure distribution around
aJoukowsky profile ...... .........10
a Joukowsky profile . . . . . . . . 10
3 Pressure distribution around a Fuhrmann body . . .. 11
4 Jeffreys' cylinder . . . . ..... . . 12
5 Pressures on a single sphere . . . . . ... 17
6 Arrangement of spheres and theoretical bed in
EinsteinEl Samni work .............. .28
7 Arrangement of spheres and theoretical bed in
Chepil's work . . . . . . . .... . 30
8 Sketches for sidewall effect . .. .. . .... 36
9 Evaluation of constant for sidewall analysis .... 38
10 Grains placed in rough bed configuration . . ... 42
11 Pressure distribution on twodimensional grain ..... 43
12 Piezometric head distribution on twodimensional grain 43
13 Distributions of surface hydrodynamic pressure decreases. 45
14 Lift coefficient CL . . . . 46
15 Theoretical bed used in relating logarithmic and
potential velocity profiles ......... .. 48
16 Definition sketch for experimental application .... 50
17 Comparison of theoretical and measured values . . .. 51
18 Solution space for closely packed hemispheres . . 58
19 Solution space for Chepil's arrangement . . . ... 59
vii
LIST OF FIGURES (Continued)
Figure
20 Foldedsymmetry boundary . . . . .
21 Sevenpoint finite difference scheme . .
22 General lattice point . . . . . .
23 Examples for object point on planar boundary
24 Arc interpolation for surface value . . .
Grading the lattice . . . . .
Xderivative condition at boundary .
Yderivatives . . . . . .
Special points for yderivative .
Zpoint for normal derivative condition
Pressures on area of hemisphere .
Chepil's case: velocities on y = 0
Chepil's case: velocities on y = 0.2a .
Chepil's case: velocities on y = 0.4a .
Chepil's case: velocities on y = 0.6a .
Chepil's case: velocities on y = 0.8a .
Chepil's case: velocities on x = 0 .
* *
37 Chepil's case: trace of some equipotential surfaces
in plane y = 0 .. ... .. . . . ... .
38 Measured and theoretical pressure distributions:
Vf = 68 cm/sec . .. . . . . . .
39 Measured and theoretical pressure distributions:
f = 91 cm/sec . . . . . . . . .
40 Measured and theoretical pressure distributions:
V = 128 cm/sec . . . . . . . .
viii
rr
r
109
117
118
119
Page
S 60
S 62
S 65
S 69
S 71
S 73
S 79
S 82
S 86
S 89
S 97
S 103
S 104
S 105
. 106
S 107
S 108
. .
* .
.*
r
. .
LIST OF FIGURES (Continued)
Figure Page
41 Measured and theoretical pressure distributions:
Vf = 159 cm/sec . . . . . . . .... 120
42 Chepil's hemisphere . . . . . . . .... 122
43 Closely packed hemispheres: velocities on y = 0 .... .127
44 Closely packed hemispheres: velocities on y = 0.2a .. 128
45 Closely packed hemispheres: velocities on y = 0.4a . 129
46 Closely packed hemispheres: velocities on y = 0.5a . 130
47 Closely packed hemispheres: velocities on y = 0.6a . .131
48 Closely packed hemispheres: velocities on y = 0.8a . 132
49 Closely packed hemispheres: velocities on x = 0 . . 133
50 Closely packed hemispheres: trace of some equipotential
surfaces in plane y = 0 . . . . . .... 134
51 Closely packed hemispheres: trace of some equipotential
surfaces in plane y = 0.5a . . . . . .. .135
52 Closely packed hemispheres: flow pattern on surface
viewed toward xyplane ... .. ..... .. . 136
53 Closely packed hemispheres: flow pattern on surface
viewed toward yzplane . . .. .. . . . .... 137
54 Velocity suppression near sidewall. ...... .. 140
LIST OF SYMBOLS
a sphere (or circular cylinder) radius
B coefficient in logarithmic velocity profile
CL lift coefficient based on total bed area
CLu lift coefficient based on projected area of grain
d grain diameter
d equivalent grain diameter
dA projected elemental area in xyplane
dA1 elemental surface area on hemisphere
g acceleration due to gravity
h full lattice increment
I point index, xdirection
J point index, ydirection
K point index, zdirection
k equivalent sand roughness
K diameter of hemispheres
s
k1 sidewall equivalent sand roughness
k2 bottom equivalent sand roughness
L hydrodynamic lift force
p pressure
p reference pressure, defined as used
pl reference pressure, defined as used
PI field point in xdirection
PI1 field point in xdirection
PJ field point in ydirection
PJI field point in ydirection
PK field point in zdirection
PK1 field point in zdirection
q total velocity at a point
qb velocity along base of the grain
r spherical coordinate
r roughnessgrain size ratio based on large hemisphere
diameter in Chepil's tests
R Reynolds number
Rew wall Reynolds number (vfk/v)
u velocity
U free stream velocity at a
u mean velocity
m
u velocity at uppermost point of hemisphere, sphere,
or elliptic cylinder
utop same as ut
u35 velocity at 0.35 Ks above theoretical bed
v velocity
v friction velocity (J7 /)
w velocity in zdirection
x Cartesian coordinate
y Cartesian coordinate; also elevation above datum
yb location of theoretical bed
yo distance at which velocity equals zero
fraction finer than (in soil gradation curve)
Cartesian coordinate; also elevation above datum
angle through "corners" of isovels
ratio of given lattice leg length to full increment, h
unit weight of fluid
spherical coordinate; also angle for arc interpolation
von Karman's constant
lift per unit area (total bed area)
lift per unit area (only projected area of grain)
kinematic viscosity of fluid
mass density of fluid
density of particle in Jeffreys' analysis
bed shear stress
potential function
value of potential function on boundary
overrelaxation factor
xii
Abstract of Dissertation Presented to the Graduate Council
in Partial Fulfillment of the Requirements for the
Degree of Doctor of Philosophy
HYDRODYNAMIC LIFT IN SEDIMENT TRANSPORT
by
Barry Arden Benedict
December, 1968
Chairman: Dr. B. A. Christensen
Major Department: Civil Engineering
Hydrodynamic lift is a force often neglected in studies of
sediment movement, despite being of the same order of magnitude as
the drag force. The goal of this dissertation is to demonstrate that
potential flow theory can be used to predict hydrodynamic lift in
sediment transport. It is known that concentration of streamlines,
rather than viscous forces, contribute most to lift. This is rein
forced by many early airfoil works using potential flow theory, which
yield good values for lift, even though exhibiting zero drag.
The theoretical methods used treat mean, steady flows with no
free surface or sidewall effects. Since potential flow involves an
inviscid fluid, the method considered is only applicable to studying
cases of real flow in the hydrodynamically rough range, where viscous
forces are negligible. The two sets of experimental work studied
involve flows in the rough range over beds of hemispheres arranged in
hexagonal patterns. One work uses hemispheres three diameters apart
center to center; the other uses touching hemispheres.
xiii
The analysis begins with a solution of the potential flow over
the two sets of hemispheres, solving here the threedimensional prob
lems by use of finite difference methods. The solutions enabled cal
culation of velocities near the hemispherical surface, relation to
pressures by Bernoulli's equation valid for rotational flow, and sub
sequent integration to find the lift force. The velocity distribution
near the surface (from potential flow) is then linked to the velocity
distribution in the actual flows considered, yielding values for lift
corresponding to the experimental works. This procedure can be used
for any rough bed and velocity field desired.
For the widely spaced hemispheres, theory produces (for three
cases most closely related to the theoretical model studied) values of
lift differing from measured values by 19 per cent, 13 per cent,.and
8 per cent. For the closely packed hemispheres, results from theory
are 16 per cent above measured values. This discrepancy is reduced to
perhaps 10 per cent or less if allowance is made for the sidewall
effects of the narrow flume used in the experimental work. Good
results are also obtained for a natural gravel bed replaced by a bed
of equal hemispheres.
The results show quite good agreement between theory and
experiment. Hence, the goal of this dissertation is accomplished, and
a new analytical tool is found effective in studying lift forces on
a rough bed. It is hoped that the tool will be useful in future work
in the field.
xiv
CHAPTER I
INTRODUCTION
One problem of great practical concern to the hydraulic engineer
is the development of a full knowledge of the transport of material in a
stream. Understanding sediment transport and developing advanced meth
ods in this area are physically and economically important in such prob
lems as scour evaluation and prevention, construction of stable channels,
planning and design of reservoirs, and the maintenance of harbor chan
nels. The problem is of true importance when some of the sediment magni
tudes involved are viewed. Brown [1]1 estimated in 1950 that the aver
age gross sediment production of the United States was on the order of
four billion tons per year. He also indicated that deposition of river
borne sediment in impounding reservoirs was estimated to be equivalent
to a loss each year of sufficient storage capacity to hold the annual
water supply for a city of 250,000 people. On a smaller scale,
Langbein and Leopold [2] note that watersheds composed of fine wind
blown soil, such as in western Iowa, yield as much as 2,000 tons of
sediment per square mile per year. The magnitude of the sediment prob
lem causes a continuing need for better understanding in sediment trans
port. One facet of this understanding must include the mechanism of
Numbers in brackets refer to the References.
a particle being moved from the stream bed. A need for a better defin
ition of the fundamental forces acting on such particles has prompted
the author's study of hydrodynamic lift in sediment transport, with
observations of flow characteristics in the vicinity of a rough bed.
Lift has been an often neglected force in fluvial hydraulics
despite the fact that it has the same order of magnitude as the bed
shear stress. This is reflected in the many stable channel design
procedures which neglect lift. These procedures may result in a design
which is too costly. Hence, studies of the lift force take on further
economic significance.
Most work on lift has been on single particles, with no con
sideration for the particle interaction existent in the actual stream.
Those works which have included interaction have generally measured
lift forces over a given area rather than in the form of pressure
distributions on the individual particles. The author intends to study
interaction effects, trying to provide a base for analytical studies
of lift forces.
First, potential flow theory will be used as a guide to provide
a means of studying the flow near the surface of a grain or bed shape
and then relating that flow to other flow characteristics. The poten
tial flow velocity profile will be related, for this work, to the
logarithmic velocity profile, furnishing a means for predicting the
pressure distribution and total lift force on single particles within
a series. It should be noted that forms of velocity profile other than
the logarithmic might be used if they are characteristic of the flow.
This dissertation will concern itself with cases involving a
steady, mean flow with free surface effects (wave resistance) being
insignificant. A report of some earlier results for twodimensional
bed shapes will be made. Then a numerical solution will be made for
two cases of flow over a bed of hemispherical particles. These solu
tions will be related to experimental measurements obtained elsewhere,
as will some of the twodimensional results.
It is the hope of the author to present herein an analytical
method wherein the hydrodynamic lift, often neglected, can be computed
with a reasonable degree of accuracy. It is also intended to form a
basis for future analytical and experimental work in sediment transport.
CHAPTER II
RELATED BACKGROUND ON LIFT FORCE STUDIES
2.1 Historical interest
For centuries engineers have been interested in the movement of
sediment due to flowing water, which Rouse and Ince call "A class of
flow phenomena inherently hydraulic in nature ." [3, p. 246].
Concerned with problems of scour, deposition of material, stable
channels, and the like, engineers have for years attempted to increase
their knowledge of sediment motion. Both empirical and theoretical
means have been employed in these attempts.
Domenico Guglielmini [3, p. 70] was perhaps representative of
the whole Italian school interested in flow resistance in open channels.
His work in the seventeenth century made some qualitative observations
which were very accurate, though his analytical work was more faulty.
Later, Pierre Louis Georges Du Buat [3, p. 1291 collected a
vast array of experimental results, including extensive data on the
beginning of sediment movement. His eighteenth century works over
shadowed that of other hydraulicians for about a century.
Work continued through various periods until the work of Grove
Karl Gilbert in the years around 1910. His tests on initial sediment
movement and various phases of transport covered a wide range. It has
been noted that ". . the results he presented in U.S.G.S. Professional
Paper No. 86 in 1914 still continue to be those most often quoted of
any in the field" [3, p. 225].
From the time of Gilbert's work, laboratory facilities increased,
enabling further, broader studies. As various governmental agencies
began to attack the problems in the United States, more and more quan
titative information became available. Additionally, interest in analyt
ical approaches to sediment problems grew with the advent of more
theories on movement. Consideration of the forces actually acting on
particles exposed to flow received increasing emphasis.
2.2 Neglect of lift force
As Leliavsky [4] indicated, the drag component of forces acting
on particles received the bulk of work earlier in this century. Indic
ative of this was the work by C. M. White [5] in 1940. White gave
great care to study of the drag, while the lift was treated only briefly,
with a guarded conclusion that lift did not exist. Leliavsky noted
in 1955 that there still exists ". . an almost unexplored aspect of
the problem, viz., the vertical component of the resultant of the
hydraulic forces applied to the grain, i.e., 'the lift'" [4, pp. 6465].
Similarly, Young stated, "Although some attention has been given
to the effect of the drag components on the behavior of suspensions,
little work has been done in connection with the determination of the
lift component" [6, p. 47].
This disregard for the lift force has resulted in a lack of
understanding of this force, with probably attendant shortcomings in
full understanding of drag. One area of very practical interest where
lift is often neglected is in stable channel design, where the widely
used method of E. W. Lane [7], as well as other methods, does not
include a lift force. It should be noted, however, that other methods
are beginning to include lift effects. This neglect of the lift force
is on the safe side but is economically costly.
2.3 Indications of significance of lift
As Leliavsky noted, contrary to White's results, Jeffreys [8]
and Fage [9] produced results which definitely indicate lift as an
important factor. Fage's experimental work provided evidence, and
Jeffreys' gave a theoretical approach. Studying a cylinder resting on
the flat bed of a deep stream and applying the principles of classical
hydrodynamics, Jeffreys found a relation for the lift and then for
a scour criterion. Use of realistic sizes indicated that lift alone
should be capable of dislodging particles.
Chang [10], considering the lift on particles as due to the
pressure of a velocity head, worked from the simple fact that the
particle would tend to lift when the vertical lifting force (or hydro
dynamic lift) plus the buoyant force equalled the particle's weight.
He also did work on drag and made the following comparison. "Theo
retically, the force required to lift a particle from the bottom of
a stream is about 40 per cent greater than that required to move it
along the bed" [10, p. 1282]. This certainly indicates the same order
of magnitude for drag and lift.
Further validity is given to the value of studying lift by
Yalin [11] who says "A consideration of the paths of saltating particles
by R. A. Bagnold reveals that saltation begins with a motion directed
'upward'; . However, if this is true, then, as has already been
maintained by various authors, the lift force . must be the cause
of the detachment" [11, p. 229]. Support for Yalin's statement is
offered by the work of Einstein [12], Bagnold [13], and Velikanow [14].
Young, doing work in 1960 on spheres in a cylindrical tube,
found the lift to be of the order of onehalf the drag for experiments
with a Reynolds number based on the pipe diameter and the mean veloc
ity in the range 3601115 (in the laminar range). He summarized,
"It is thus apparent that the lift force should not be overlooked in
studies related to the incipient motion of particles resting on stream
beds or pipe walls" [6, p. 57].
It seems apparent that sufficient evidence exists to prompt
efforts to increase understanding of the lift phenomena. The author
will now consider earlier works on lift.
2.4 Use of potential theory in lift studies on single bodies
2.4.1 General.Much of the analytical work done on lift has
dealt with single bodies or particles, especially spheres and circular
cylinders, isolated from any other bodies. An example of this approach
is the work done by Jeffreys [8] mentioned earlier.
The overwhelming amount of work on lift has been done in the
realm of aerospace engineering. While the principles thus developed
are applicable to hydraulics, the work concerns single bodies only,
often airfoils, struts, and the like. One of many works giving dis
cussions of several of these approaches is the famous work edited by
Goldstein [15]. Another, viewing historical development, is by
von Karman [16]. In a number of these, lift is studied by means of
potential flow.
2.4.2 Applicability of potential flow.The idea that lift
could be predicted by use of potential theory, even perhaps in those
areas where it might seem out of place, was offered by Prandtl. He
stated that any explanation of drag requires a consideration of
viscosity, ". . whereas the lift can be explained entirely without
the concept of viscosity so that the wellknown methods of the clas
sical hydrodynamics of the ideal fluid are applicable" [17, p. 159].
Part of the reason for some reluctance to use such methods here is the
fact that within a boundary layer adjacent to a surface, irrotational
flow does not exist. However, in the cases of interest, the sublayer,
where viscous and inertial forces are of the same order, over the surface
is very thin. Outside this layer, up to the point of separation, the
equations of inviscid fluid flow are valid. Hence, the pressures on
the outer edge of the sublayer can be found from such equations, and
since the pressure difference from the outer edge of the layer to the
surface is assumed negligible, the normal pressures on the surface, and
thus lift, can be predicted by inviscid flow principles. Also, since
the sublayer is thin, it is possible to discuss velocities "on the
surface," though the velocities actually considered are those a small
distance away at the edge of the boundary layer. Of course, the
actual velocities at the wall in a real fluid would be zero.
Reasoning such as Prandtl's has through the years prompted
many potential flow solutions to the lift problem, probably the most
famous of which are those by Kutta and Joukowsky.
2.4.3 Work of Kutta and Joukowsky.Generally, an understand
ing of the work of Kutta and Joukowsky begins with flow around a cir
cular cylinder. It is well known that no lift or drag forces are
predicted by any transformation of the symmetrical flow around a cylin
der. The addition of a vortex in the center of the cylinder produces
a streamline pattern, the effect of which is to yield a lift force.
This force is due to the circulation's tendency to increase the velocity
above the cylinder and decrease it below, thus causing a pressure dif
ference and a consequent lifting action.
Kutta [18] first applied the methods of conformal transforma
tion to transform the cylinder with circulation into a line inclined
to the flow. This produces a force acting perpendicular to the veloc
ity at infinity. This force is, of course, the one originally exerted
on the cylinder.
Kutta's work and use of transformations prompted further work.
Joukowsky [19] wanted to avoid difficulties at the sharp leading edge
of Kutta's plane, and he employed a mapping function by which a curvi
linear profile very similar to actual airfoil shapes was developed.
Flow around a Joukowsky profile is shown in Figure 1.
Numerous investigations of Joukowsky profiles have been carried
out. Joukowsky [20] himself performed experiments in 1912, and
Blumenthal [21] calculated the pressure distribution from theory
in 1913. Betz [22] provided experimental comparisons for the lift
and pressure distributions predicted by theory. Figure 2 shows this
comparison, and the agreement with theory is seen to be satisfactory.
The lower lift than that predicted can be accounted for by the fric
tion which causes flow separation from the profile near the end; the
resultant failure to attain the full pressure difference predicted
yields a smaller lift. The overall agreement is, however, good.
Figure 1.
Flow around a Joukowsky profile.
Profile
Calculated
Measured
Figure 2. Calculated and measured pressure distribution
around a Joukowsky profile. (From reference 17,
p. 181.)
The classic work of Kutta and Joukowsky, with subsequent
experimental verification of the validity of their approach, has led
to many other lift solutions by means of potential flow. These works
have generally produced similarly satisfactory results.
2.4.4 Fuhrmann's work and other studies.Prandtl and Tietjens
[23, p. 137] mention the work of Fuhrmann [24] who calculated, by poten
tial theory, and then measured the pressure distribution on some slender
bodies whose shapes were derived from sourcesink combinations. As
found in reference [15] and pictured in Figure 3, Fuhrmann's experi
ments revealed good agreement with theory except near the body's end, a
result common to other works and also expected. Rodgers [25] notes
this result in reference to several bodies of revolution analytically
treated by Lamb [26] and MilneThomson [27]. He attributes the devia
tions from theory primarily to generation of vorticity along the body.
i.
0.8  Measured
_ Potential Flow
0.6_
PPo i I
2
p v /2 1
o 0.2
0
0.4 x
0 0.2 0.4 0.6 0.8 1.0
x/L
Figure 3. Pressure distribution around a Fuhrmann body.
Other works where pressure distributions have been predicted
and checked experimentally could also be shown. The important factor
is not to go into details of numerous cases, but rather to point out
that there is strong historical backing for use of potential theory
to study lift forces. It should be recalled that these earlier uses
of potential theory have employed it to describe the entire flow
pattern. Only two further cases will be studied, the first being the
one case where potential flow was used to describe the entire flow
pattern in an effort to solve the problem of sediment movement.
2.4.5 Jeffreys' analysis.Jeffreys [8] dealt with a single
long circular cylinder resting on a bed in a twodimensional study,
as shown in Figure 4.
Uniform Flow
U
a Long Circular
Cylinder
Figure 4. Jeffreys' cylinder.
Jeffreys developed his work based on the complex potential for
this case,
na
W = raU coth () (21)
z
He found the amount by which the lift exceeded the weight of the
enclosed liquid, and thence wrote the following condition for move
ment of the cylinder:
( + 2 )U2 > ga (22)
3 9 p
where a = density of particle
g = acceleration due to gravity
a = cylinder radius
p = mass density of water
This can be written as
U2 > 1 ( ) ga (23)
1.43 p
Jeffreys noted that J. S. Owens [28] had earlier measured
the velocity required to move pebbles, finding
U2 = 1.65 ga (24)
The motion observed by Owens was not, however, a jumping or
lifting motion, but rather a rolling motion, yielding a difference
from values predicted by theory. However, as Jeffreys says, "The
proportionality of U to the square root of the linear dimensions is in
agreement with theory" [8, p. 276].
Use of a value for 7 of 2.7 times p in Jeffreys' equation yields
the following:
U2 > 1.19 ga
(25)
Differences from theory are primarily due to two causes.
First, the motion measured was rolling rather than totally lifting,
which (25) is based on. Second, the measured particles were three
dimensional, contacting the bed at only a small number of points,
while Jeffreys' cylinder made contact over the whole length of its
body. These effects are compounded by the fact that the Uvalues used
by the two men are not the same. Jeffreys employed a potential flow
field, but Owens made his measurements in a flow field exhibiting a
logarithmic velocity distribution. The differences involved here will
be discussed in Chapters III and IV. Despite the discrepancies, the
theory seems to provide a much better starting point than might first
be thought possible.
2.4.6 Flow around a single sphere.The results from measure
ments on flow around a single sphere will be presented because the shape
relates to this study and the results reveal some factors influencing
the actual flow pattern.
The flow involved is that around a single sphere suspended in
an otherwise uniform flow. For this case the potential can be
expressed [29] as
Ua
S= cos e + Ur cos e (26)
2r
2r
where a = radius of sphere
r = radial distance
0 = angle measured from horizontal and
through sphere center
U = free stream velocity at 0
Note that this is written for one meridian plane, since the
case is axisymmetric. The tangential surface velocity on r = a
can be found as
1 3
= U sin 8 = q (27)
Using Bernoulli's principle, with a pressure of po at the
forward stagnation point, yields the following expression for the
surface pressure, p.
p p = yz + q (28)
where z = elevation
y = unit weight of fluid
What is of interest here is to study only the vertical force
component over the upper half of the hemisphere. This will involve
integrating (28) over the upper surface. Integration of the first
term, z, would yield a hydrostatic lift (buoyancy) which would, in
fact, not even be measured, as the difference in piezometric head is
the measured value. Hence, integration of the last term in (28)
will yield the desired force component. It should be noted that this
would actually be the theoretically predicted value for the case
of a single hemisphere on a flat bed. However, the measurements herein
used are for a suspended sphere. The total vertical force can be
found as follows:
F = U2 sin2e dA (29)
v 2 4
o
where dA = 2a2 sin28 dG
yielding
2
S= a p 2 (210)
Note that if this was expressed per unit area (based on the area
na ), with U replaced by ut = (3/2)U, the velocity at the top of
the sphere, theory would show
2
3 u
3 t
F =  (211)
p 4 2
where F = vertical force per unit area
Measurements were made by Flachsbart [30], as shown in
Schlichting [31]. The results were very similar to those by Fage [32]
shown below. The measured values are shown below with the theoretical
curve. The values are plotted with reference to the stagnation pres
sure pl + ~ pU where pl is the pressure at where the velocity
equals U. This stagnation pressure is the p being used here. Replac
ing pl by po yields
S o= B 1 (212)
1 2
2pU
where B is the value shown in the graph. The expression thus found
for (p po) can be integrated numerically.
In the case of the higher Reynolds number, the force can be
expressed as (210) with a coefficient of 1.64, only about 3 per cent
below the theoretical. The lower Reynolds number produces a coefficient
of 1.32, about 20 per cent below the theoretical. It therefore seems
that above some critical Reynolds number the potential theory becomes
more and more capable of predicting pressures on a body outside the
separation zone. The effect of the Reynolds number involves the change
of laminar boundary layer at the body to turbulent as a critical range
of Reynolds number is reached. The boundary layer, becoming turbulent,
is capable of proceeding further downstream, therefore causing separa
tion to be delayed to points further and further back on the sphere.
Ppl
 0
1 2
pU
0 60 120
180
0 (degrees)
Figure 5. Pressures on a single sphere.
This allows the surface pressure distribution to more nearly approach
the theoretical. For even lower Reynolds numbers than those indicated,
the viscous forces would play an even greater role, thereby causing
further deviation in the flow from theoretical. These ideas are of
importance in relation to the limitations of work in this dissertation.
The forces evaluated in this section were only those in the
vertical direction on the upper surface of the sphere, with an eye to
noting agreement with theory.
For later comparison a value of CLu will be determined here
for the theoretical case of a single hemisphere on a flat bed. This
will entail integrating the pressure at the hemispherical base, found
from (28) and subtracting it from the force of (211) to give a
resultant vertical lift. Integration here occurs in a direction normal
to that in (29). Using 81 for this integration yields a relation
cos 0 = sin 0. Therefore, the ratio of the basal velocity (ub) to the
top velocity (u t) equals cos 81. The pressure decrease along the
bottom can be found as below.
2 2
t 22 2 1 t
Fv = P T (2a2) sin29 cos2 dO = 0 P T a2 (213)
The lift per unit area is therefore
F F u
S vvv vb t1 t
S= ,a = p  (2P.)
u 2 2
Tra
This lift.coefficient of 0.50 should form an upper bound for the
work to be done later.
At this point it seems appropriate to define certain lift
coefficients to be employed in this dissertation. The difference
lies in the area over which the force is considered. The subscript,
u, will denote those cases where the area considered is only that
directly beneath the body being considered, the projected area of the
grain. Thus, \ denotes lift per unit area based on the total bed
area, while X is based on the area of the grain projected onto the
U
bed. Similarly, the coefficients CL and CLu are used with the corre
sponding X's.
As indicated, most work done on lift has dealt with single
bodies. Attention will now be turned to systems with more than one
particle.
2.5 Multiparticle studies
2.5.1 Einstein and El Samni.It was natural that multi
particle studies should arise, as these begin to approach the sediment
conditions found in nature. Unfortunately, however, work in this area
has been limited. Some values for the lift force came from the work
by Einstein and El Samni [33,34]. Using the upper onehalf of
plastic spherical balls 0.225 feet in diameter placed in a hexagonal
pattern, they measured the lift force as a pressure difference. They
made the following statement.
The procedure in making such measurements was as follows:
if a lift force is exerted on the top layer of a stream bed,
the solid support of the sediment particles is relieved of part
of their load and this load is transmitted hydrostatically to
the fluid between the solid bed particles. Thus, it must be
possible to detect and measure this lift as a general pressure
increase of the pore fluid in the bed. [33, p. 521]
Their results enabled them to write
2
p = CL P (215)
where Ap = pressure difference
CL = lift coefficient
p = fluid density
u = velocity
They found CL a constant 0.178 if u was taken as the velocity 0.35
sphere diameters above the theoretical bed, determined by experiment
as 0.20 sphere diameters below the sphere tops. Further studies
which they made on natural gravel yielded the same expression for Ap
with some redefinition of u along lines consistent with Einstein's
earlier work [12]. This work forms the only example of a lift force
essentially integrated over a number of particles, though the distri
bution over individual particles was not ascertained. Support for the
rationale of measuring lift as a pressure difference is given by
Engelund and Hansen [35, p. 19] in discussing variations from hydro
static pressure due to streamline curvature.
2.5.2 Chepil.Chepil [36] performed experiments in a wind
tunnel on hemispherical elements placed on a plane bed in a hexagonal
pattern three diameters apart. He chose this spacing based on work
by Zingg [37], which indicated this is the average spacing between
particles erodible from a sand bed.' The processes of erosion of sand
by wind and by water involve essentially the same factors, as indicated
by, among others, Kadib [38].
Chepil measured pressures over the surface of one metal
hemisphere. This was done by placing pressure taps, starting at the
base of the hemisphere, 30 degrees apart along one line running parallel
with and another normal to the wind direction. The negative pressure
end of the manometer was connected to a tap on top of the hemisphere.
The remainder of the hemisphere pattern consisted of gravel hemispheres.
The hemispheres occupied 11 per cent of the total floor area. The lift
and drag forces on the hemisphere were determined by integrating the
measured pressure distributions and also, as a check, by means of two
torsion balances measuring the forces directly. Hemispheres of three
different sizes were used, and also some measurements were made on
relatively small sand and gravel mounds. Measurements were also made
at different points downstream in the tunnel to note the effects as the
air boundary layer developed to its full extent.
Chepil found that increasing the depth of the fluid boundary
layer, after a certain limiting depth is reached, has little effect on
the lift to drag ratio, though the depth of boundary layer had a
profound effect on the magnitude of both lift and drag. For the study,
it was found that
LIFT = 0.85 Drag (216)
For this study it was also found that the pressure difference between
the top and the bottom of the hemisphere is about 2.85 times the lift
per unit bed area directly under the hemispheres. Using the latter
finding, Chepil surmised that the CL from Einstein and El Samni [33]
should equal 0.178/2.85, or 0.0624. He then used equation (216)
with his drag and velocity measurements to attempt a correlation with
the case of closely packed hemispheres and found a CL of 0.0680.
Chepil concludes by noting that "This study shows that lift on
hemispherical surface projections, similar to soil grains resting on
a surface in a windstream, is substantial. Therefore, lift must be
recognized together with drag in determining an equilibrium or crit
ical condition between the soil grains and the moving fluid at the
threshold of movement of the grains" [36, p. 403].
2.5.3 Chao and Sandborn.Chao and Sandborn [39] also performed
experiments on spheres, but, by means of a transducer, they actually
measured the pressure distribution on the upper half of an element.
However, the type of flow they used bore no real relation to that in
nature's streams, and no attempts were made at analyzing velocities and
the like. Einstein used the flow of water in a channel with a measured
logarithmic velocity distribution. On the other hand, Chao and Sandborn
placed lead shot on a flat surface and blew a stream of air down onto
the particles, the air diverting horizontally at the flat surface.
Their conclusions included, "The present experimental results are of
primary interest in demonstrating that a problem exists. . More
extensive research is needed before there can exist a better understand
ing of the mechanism" [39, p. 203].
2.6 Shapes of bodies studied
With the exception of the vast amounts of work done on wings,
struts, and related areas, most effort has beenaimed at circular
cylinders and spheres, especially in work related to the sediment
problem. There have been, however, numerous indications of shape as
a parameter. Rouse defines sphericity [40, p. 777] and indicates it
as a probable factor. Studies have been made on shape influence,
including one by Krumbein [41] which specifically considers nonspherical
particles. Other works also exist for considering the effect of non
spherical particles.
Most design procedures today treat the bed material by conver
sion to a representative bed consisting of uniform spherical particles.
The conversion is generally based on factors such as particle density,
sizes, exposed areas, and volumes. Thus, there are means for relating
natural beds to the uniform spheres to be treated in this work.
2.7 Relation to work of dissertation
The preceding background material was presented in an effort to
point out those factors with special bearing on what the author is try
ing to accomplish in this dissertation. Generally, an attempt will be
made to advance the knowledge of lift through analytical work.
Analysis of the threedimensional cases corresponding to the
work of Einstein and El Samni and Chepil will be made. These multiple
particle systems allow consideration of particle interactions. Poten
tial flow theory will be used to study velocities on the surfaces, these
velocities then being linked to known logarithmic distributions.
Additionally, some twodimensional work and applications will be reported.
Support for the use of irrotational flow is gained from the lift predic
tions made on airfoils, as well as from a knowledge that viscosity is
not necessary to explain lift. Of course, in the airfoil cases as well
as here, the predicted drag is zero, obviously incorrect, but the
description of the magnitude of the lift is the primary goal here.
It should be very strongly noted, however, that whereas these earlier
works used the potential to describe the total flow field, a much
different use will be made herein. Essentially, potential flow theory
will be used as a guide in attempting to develop means of predicting
hydrodynamic lift. Success here might lead eventually to better
analytical approaches to some sediment problems, although it will
also be used in larger scale problems dealing with dunes and ripples.
Summarizing, it can be said that the author is relying on the
success of earlier approaches to related studies in an effort to find
a guide to a better understanding of hydrodynamic lift in sediment
transport.
CHAPTER III
LOGARITHMIC VELOCITY DISTRIBUTION
3.1 Development of logarithmic velocity distribution
Due to the use which will be made of the velocity distribution,
the author wishes to present briefly some background and to study some
specific points. Historically, three modern approaches to velocity
distributions in steady, uniform turbulent flow have arisen. The three
are the following: Prandtl [42), who introduced the concept of the
mixing length (related to the mean free path of particles) with momen
tum conserved; G. I. Taylor [43], who considered vorticity to be con
served along the mixing length; and von Karman [44], who developed a
similarity hypothesis for the problem.
Prandtl's derivation, beginning with the expression for shear
stress in fluid, is frequently cited in texts, such as [45]. It
neglects the viscous forces and considers only the socalled Reynolds
stresses, after 0. Reynolds [46]. The final expression can be written
as below.
V = B + ln (31)
v f k
where B = C + In k
k = equivalent sand roughness
The socalled von Karman constant, K, has been verified
experimentally, as noted by Bakhmeteff [47], to equal 0.40. The
values for B were obtained from experiments by Nikuradse [48,49],
who made tests on smooth pipes and then on pipes roughened by gluing
uniform sand grains of size k to the wall. Results in the transition
range between smooth and rough flow regimes by Colebrook [50] showed
different Bvalues for nonuniform roughness.
Hydrodynamically smooth and rough flow regimes are generally
defined in terms of the wall Reynolds number (vfk/v), Rew. The smooth
f ew
range, Rew less than 3.5 to 5, has a viscous sublayer of depth compar
able to bed roughness size and viscous forces play a dominant role.
In the rough range, Rew greater than 70, the roughness elements have
fully penetrated the sublayer, and viscous effects are negligible.
Most sediment problems in nature occur where flow is in the
rough range. This will be the range considered in the present work.
Substitution of the appropriate values for the rough range yields
S= 8.48 + 2.5 In (32)
vf k
Although this expression was derived from work in circular
pipes, it has been noted on many occasions that it is similarly
applicable in open channels.
3.2 Special problems of present expressions
Working with measurements near rough surfaces, Hwang and
Laursen noted that ". .. the adequacy of this logarithmic equation
can be debated, especially near the boundary, and the indeterminacy
of the zero datum presents difficulties. ." [51, p. 21]. These
essentially are the two problems with use of (32).
3.2.1 Theoretical bed.In a bed in which roughness elements
protrude, some reference datum must be decided upon. This reference
is usually referred to as the theoretical plane bed, and, if sufficient
information is available, its location can be computed. Generally,
satisfactory results are obtained by use of this datum, though problems
may arise for very irregular surfaces.
The basis for computation is to place the bed at such a level
that the water volume contained beneath it equals the volume of rough
ness elements protruding from the bed. Justification for use of this
definition can be gotten by reviewing the work of Einstein and El Samni
[33]. Using hemispheres of diameter k and placed in a hexagonal pat
tern, they found that for satisfactory results it was necessary to
assume the theoretical bed at a distance of 0.2 k below a tangent to
the tops of the spheres. Their arrangement appeared as shown in
Figure 6.
Realizing that the area considered (inside the dotted lines)
has dimensions 3k by v/3 k the following is obtained.
s s
Volume under theoretical bed = 3, k2 (33)
s b
1 3
Volume of particles =2 s k (34)
2 s
Equating these two and solving yields
(35)
yb = 0.302 ks
Plan
1 Flow
Theoretical Bed
y = 0 in equation (37)
0.2 ks by Einstein
and El Samni
Elevation
Figure 6. Arrangement of spheres and theoretical bed in
EinsteinEl Samni work.
Therefore, the distance down from the upper tangent is
0.198 k in excellent agreement with the 0.2 k found by experiment.
It might be noted that the value of 0.2 k seems to be satisfactory
even when the roughness pattern is more irregular. This was shown
experimentally by Einstein and El Samni [33]. Other instances could be
reviewed, but let it suffice to state that the above definition of the
/
theoretical bed is satisfactory.
At this point the conditions for the Chepil [36] case are shown
in Figure 7. The computations for the theoretical bed are also shown,
since this will be needed later. Considering the region inside the
dotted lines, with dimensions 3V1 Ks by 3Ks, the following is obtained:
Volume under theoretical bed = 9, K Yb (36)
s b (3)
1 3
Volume of particles = K (37)
6 s
Equating these and solving yields
Y = 0.0335 K (38)
b s
Therefore, for this pattern the distance down from the upper tangent to
the hemispheres is 0.467 K very near the plane bed itself.
3.2.2 Conditions near the bed.Experience indicates that
expression (32) is quite valid away from the bed. However, inspection
reveals that problems occur near the wall. As y approaches zero, the
value of v/vf, and hence of v, approaches minus infinity. This obvious
discontinuity presents a real problem in understanding what happens at
the bed, where sediment movement begins.
0
0
Flow
r
I '1
0
lan
Flow
__rw

'1
3K
I
Theoretical Bed
ri_...L\
Elevation
Figure 7.
Arrangement of spheres and theoretical bed
in Chepil's work.
P
3.3 Use of proposed adjusted velocity distribution
Christensen [52] is suggesting use of a slightly different
velocity distribution which he has developed. The new expression is
as follows:
v = 8.48 + 2.5 In (Y + 0.0338) (39)
vf k
Calculation will reveal one very desirable feature of (39),
that being prediction of a zero velocity for y equal to zero. The
added constant factor will continue to have a large effect very near
the wall, but as y/k increases, the effect of the additional term will
quickly become negligible, yielding essentially the same expression
as (32). This is correct, since (32) has proved adequate away from
the wall.
3.3.1 Comparison of distributions at wall.In addition to
the effect of zero velocity at the wall as compared with the infinite
value predicted by (32), other features can be noted by looking at
the change of velocity with distance from the wall, or dv/dy shown
below.
dv f
Former: = 2.5
dy y (310)
For y 0 dv/dy 
y O dv/dy 0
dv 1
Proposed: = 2.5 v ] (311)
dy f y + 0.0338k
For y = 0, dv/dy = 74.0 v /k Finite
y dv/dy 0
In both cases, the derivative approaches zero for large
yvalues, as it should. However, at the wall, a definite discontin
uity exists by the former method. The new proposal predicts a finite
change at the bed, a far more reasonable development. The author
feels that the proposed distribution, with its continuous curve, will
enable better studies of action near a rough bed.
3.3.2 Comparison of distribution with increasing y.The
increasing value of (y/k) will eventually negate the effect of the
added term in (39). Table 1 shows the difference in velocity indi
cated'by (32) and (39) for some (y/k) values. The difference in
hydrodynamic force, proportional to the velocity squared, is also
shown. It can be seen that the effect of the added term in (39) is
dissipated very rapidly.
TABLE 1
COMPARISON OF PROPOSED AND FORMER DISTRIBUTIONS
Variation between Force
(32) and (39) y/k Variation
1% 0.99 2%
3% 0.61 6%
5% 0.29 10%
3.4 Determination of k and vf from experimental data
Essentially, in either equation (32) or (39), it is neces
sary to find k and vf to fully describe the profile in a given
situation. If the distribution is first measured to be logarithmic,
k and vf may then be determined by considering the velocity at two
different depths [45].
Call the total flow depth d. Consider then a depth pd (p,
a fraction), where the velocity is v Consider also a second depth
p
qd with velocity v .
k, the equivalent sand roughness
vf, the friction velocity
From these measurements the following can be written, using first
expression (32).
Depth 1:
v
2 = 8.48 + 2.5 In Md (312)
vf k
Depth 2:
v
= 8.48 + 2.5 In q (313)
vf k
Combination of these equations and elimination of vf yield
29.7 pd (314)
(v )/(v v )
( p '
Elimination of k yields
V. v
vf = (315)
2.5 In (P)
q
Frequently the depths chosen are for p = 0.90 and q = 0.15.
Use of the proposed expression in (39) for equations such as
(312) and (313) above, does not lend itself as readily to direct
elimination of variables and solution for k and vf. Therefore, some
consideration must be given to the relative values of the terms.
As noted earlier, for terms where y/k is greater than one, the effect
of the added 0.0338 is less than 1 per cent and can presumably be
neglected. If indeed p = 0.90 and q = 0.15 were used, for most
practical cases y/k would be far greater than one for both depths,
and the expressions for k and vf would be precisely those found in
(314) and (315).
In any case, the mechanism is available for computation of k
and vf through observations. These parameters can then be used in
computing forces, velocities, and the like according to equations
to be subsequently derived.
3.5 Effect of sidewalls
As will be seen later, the experimental flume to be studied is
rather narrow. For this reason, possible effects of the sidewalls on
the velocity will be discussed here.
Sayre and Albertson mentioned this problem, when, in reference
to their analysis, they indicated that one assumption which they made
was . that the channel is of sufficient width, or that the bed
roughness is so great relative to the sidewall effect" [53, p. 124].
Rouse [54, pp. 276277] points out the effect that sidewalls have,
through secondary flows, in varying the isovels and in depressing the
region of maximum velocity below the surface of the flow.
In order to provide some quantitative means of evaluating the
sidewall effect, aside from the qualitative evaluation of experimen
tally obtained velocity profiles and isovels, the following approx
imate analysis is presented. Figure 8 indicates the assumptions util
ized listed below.
Assumptions:
1. Idealized isovel picture, enabling passing of
line through corners.
2. Idealized pressure distribution (varying linearly,
as indicated in Figure 8).
3. T1 proportional to yl in some way; T2 proportional
to Y2 in the same way.
4. Logarithmic velocity distribution in center line
profile.
Roughness:
Idealized
Isovels
Roughness:
Probable True
Tdistribution
Idealized
Shear
Distribution
Sketches for sidewall effect.
Figure 8.
The following expressions can be written relating the shear
stresses.
1.0 1
=  = tan f (316)
2.0 Y2
T1.0 h + T2.0(W h tan a) = yh WSb (317)
This can be written as
T1.0 + tan ') = WSb (318)
Equation (316) results from assuming 71 and 72 to both be
proportional to.their respective yvalues in the same way. The other
expressions merely equate the weight of fluid acting on the wetted
perimeter and the shear forces resisting it.
Obviously the value of tan Q0 is of primary importance'in this
analysis of sidewall effect. The equation enabling evaluation of the
angle comes from equating expressions for the velocity at a point such
as A in Figure 8, computed first with reference to the channel bottom
and then with reference to the side. For this analysis the form (32)
will be used for the velocity.
[8.48 + 2.5 In [ ] A = [8.48 + 2.5 In []} T (319)
k1 2 p
Rewriting this and using the similar proportionality of
T1 and T2 to yl and y2 yields
[In 29.7 ] = [In 29.7 .] L 2 (320)
1 2
100
F(M)
'I = y/k
Evaluation of constant for sidewall analysis.
100
Figure 9.
Introducing y/k = T enables the writing of
,/I F(I,) = /F(T2 ) (321)
where
F(C) = [In (29.7 1)] 7 /T
The curve shown in Figure 9 indicates that the Ffunction can
be approximated by a function of the form
F(J) = alb (322)
This would enable the writing of
1/b
1 27 (y_1/k) y_ k2
S_ 2 =(323)
Sy2/k2) 2 1
Using the fact that = tan C, it can be written that
Y2
l(1/2b)
tan a = () (324)
2
Evaluation from Figure 9 reveals that bcequals about 0.664. Thus,
0.248
k2
tan = 1) (325)
The last equation provides an opportunity for an
approximate evaluation of the sidewall's effect on the isovels and
hence on the center line profile. The likelihood of simulating two
dimensional flow in the center region can then be estimated.
40
Essentially, the work above is an attempt to check beforehand
the probable validity of all the previous discussion on logarithmic
distributions in this chapter. The chapter, in general, was intended
to point out phases of the velocity distribution applicable to the
present work.
CHAPTER IV
TWODIMENSIONAL WORK: EARLIER RESULTS
4.1 Shapes studied
In earlier work the author has studied some twodimensional bed
shapes. In one case a convenient potential function (that formed from
an infinite row of vortices) was chosen, and the grain shape produced
by this was taken. Here the solution was in a closed mathematical form.
In other cases, a series of elliptic cylinders was chosen, and the
solution for the irrotational flow obtained, in terms of the stream
function, by finitedifference methods. Some of the results obtained
and methods used will be noted briefly here as a prelude to the three
dimensional work of this paper.
4.2 General methodsvelocity and pressure results
The potential flow solution for flow around the given shape was
first obtained, using Laplace's equation involving the stream function.
The solution was used to compute velocities along the upper surface of
the shape. Then, applying Bernoulli's equation for a streamline along
the upper surface streamline, the pressure distribution could be found.
The grain shape obtained from the potential of an infinite series of
vortices is shown in Figure 10, and Figures 11 and 12 show the computed
pressure distribution on these repeating grains. It can be seen that
Figure 11 indicates, in the difference between the upper curve and its
Main Stream Flow
cosh Y = 2 + cos X
PO
vT
/ Stagnant
/ b/d = 0.561
\
NS.
Figure 10. Grains placed in rough bed configuration.
c,,
Figure 11.
()bottom
Y bottom
PO
(t)
y top
Pressure distribution on twodimensional grain.
z+
Y
)
(z + =bo m
y bottom y
+7T
(z + )
y top
Figure 12. Piezometric head distribution on twodimensional grain.
P
Y
,,,
dotted counterpart, the buoyant lift, while the remaining pressure
decrease on the upper surface contributes to hydrodynamic lift.
In the experimental work, the actual quantity to be measured
is the difference in piezometric head between upper and lower surfaces.
Since the situation below the grain is hydrostatic, the difference can
be taken between a point on the surface and any point in the region
below.
Results for velocitysquared distributions for the elliptic
cylinders can be found in Figure 13.
4.3 Lift integration
Integrating the vertical pressures over the grain surface will
yield the lift. The lift per unit area, X, can then be expressed by
S=c 2 u2 (41)
L2 t
where CL = lift coefficient
ut = velocity at top of grain = utop
This lift coefficient was found to have a value of 0.500 for the hyper
bolic cosine grain shapes, and the results for the elliptic cylinders
are found in Figure 14. Note the parameter b/a represents the thickness
tolength ratio of the ellipse. The numerical solution of the Laplace
equation requires a completely bounded solution space. Studies of flow
solutions for single elliptic cylinders and for the hyperbolic cosine
grains indicated that the streamlines became essentially straight lines
45
1.0
0.3
22
utop
0.8
5.0 \3.0 \ 2.0 \
0
0 1.0
x/a
Dashed Line: cosh bed shape
Figure 13. Distributions of surface hydrodynamic pressure
decreases.
0.5 1.0
Figure 14. Lift coefficient CL.
Lj
5.0 b/a
at distances of two to three times the particle height. Hence,
a solution space was chosen which was four times the height of the
particle in question.
4.4 Calculation of predicted lift by use of
measured logarithmic velocity profile
Figure 15 illustrates the basis for relating the potential
predictions to the actual turbulent flow. The general procedure can
be outlined in a few steps. These steps are cited under the assump
tion that it is established that the turbulent velocity profile is
indeed logarithmic in nature and that measurements have been made.
1. For the relationship (39)
yYb
S= 8.48 + 2.5 In ( + 0.0338)
f k
evaluate vf and k as outlined in Section 3.4. Notice the use of
(y yb) to indicate that this refers to distance above the theo
retical bed.
2. Then, using (39) at the top of the grain, evaluate ut
from the measured data.
3. Calculate the lift by equation (42).
2
u
Unit area: X = CL p (42)
Thus, by relating the velocity profiles, it is possible to
make a prediction of the lift.
Potential
Velocity
Profile
Logarithmic
Velocity
Profile
Theoretical Bed
Figure 15. Theoretical bed used in relating logarithmic and potential velocity profiles.
The discussion above and Figure 15 indicate the case for the
hyperbolic cosine grains and a logarithmic velocity profile, but any
other.bed shape or velocity profile could have been involved.
4.5 Application to experimental results
At the time of the earlier work, the author referred to some
measurements made by Vanoni and Hwang [55] over a part of an alluvial
bed which had been allowed to form twodimensional ripples under flow.
The bed was artificially stabilized and then pressure measurements were
taken. Measurements were used corresponding to Vanoni's C series, with
a flow depth of 0.350 feet above the mean bed.
The same methods as indicated earlier were used, solving first
for potential flow in the solution space indicated in Figure 16. It was
assumed, for ease of solution, that the streamlines had horizontal
tangents at the two ends of the space. The solution also was developed
using the dividing streamline in both cases as a portion of the bound
ary, thus assuming a region of no flow beneath the streamline. Fig
ure 17 indicates the results of the computations, which seem to approx
imate the measurements well.
This chapter was presented to provide a link between the
earlier work and the work of this dissertation,
Flow
I I
Si SI
I DLviding Streamline ,
......... I d
Actual Bed
Figure 16. Definition sketch for experimental application.
0.8
Theory  
Measured
pp
1 2
pu
P um \\
0.4 
S!/
0.0
0 0.2 0.4 0.6 0.8 1.0
Pr: reference pressure Bed Location
u : mean velocity
m
Figure 17. Comparison of theoretical and measured values.
CHAPTER V
THREEDIMENSIONAL NUMERICAL SOLUTIONS
5.1 General
Previously it was indicated that irrotational flow theory would
be utilized to study velocities on a threedimensional bed surface
toward the end of making an analytical evaluation of hydrodynamic lift
on such a surface. The results of some previous similar twodimensional
studies were presented in Chapter IV. In this dissertation emphasis is
placed on a comparison of the analytical results obtained by the writer
with those obtained by Einstein and El Samni [33,34] and Chepil [36].
Evaluation of lift in these cases by the proposed method will require
solution of the potential flow equation for the two arrangements of
hemispheres as indicated in Figures 6 and 7.
A potential flow solution implies the solution of Laplace's
equation (51) in the given flow space.
2 2 2
2+ 2 + 2 = 0 (51)
2 2 6 2
bx by z
Laplace's equation is of elliptic type and hence its solution
is fully determined by conditions on the boundary enclosing the solu
tion space. On the boundary must be specified either values of the
potential function, cp, or its derivative normal to the boundary.
5.2 Problem formulation
5.2.1 Choice of solution method.There are numerous ways of
attempting a solution for Laplace's equation in a given case.
Robertson [56] gives a listing of several such methods, some of which
were considered for the solutions to be made herein.
The method of separation of variables was eliminated because
no convenient coordinate system is available to represent the boundary
involved in these cases. The method of integral equations,in which the
differential equation (51) is expressed instead in integral form and
a solution made,was eliminated because some past works indicated that
threedimensional studies by this means had resulted in exorbitant
computation time on the computer. Another method which has been
employed is the method of integral transforms where the number of var
iables in the equation is reduced to two, and the resulting two
dimensional problem solved by numerical methods. These twodimensional
results are then inverted back to the threedimensional case to obtain
the desired final result. Tranter [57,58] shows some examples of this
approach. This method, however, has not been frequently used, for it
is dependent upon finding the proper transform kernel, which is not
usually possible. The method frequently used to obtain desired body
shapes in a flow is that of the summation of singularities and their
effects, such as sources and sinks distributed throughout the flow
volume. The distributions and relative strength of these singular
ities would be varied until their accumulative effect produced a surface
across which no flow occurred that was sufficiently similar to the
54
hemispherical surface being treated here. While the latter method
showed promise, there were numerous problems involved, such as a choice
of distribution and types of singularities and possible stability
problems in the numerical process required to vary these factors and
approach a solution. Having considered these other alternatives, the
author finally chose the method of finite differences for the solution
since this method enables obtaining the desired accuracy, while, at
the same time, providing a mathematically stable numerical analysis.
The potential flow case being considered is one in which the
depth of flow is essentially infinite. However, the choice of finite
differences method of solution makes treatment of such an infinite flow
space impossible. Therefore, the geometry of flow patterns of other
cases was reviewed with an eye to choosing a depth of flow which would
enable solution of the problem without adversely affecting the desired
results.
5.2.2 Depth of flow space.Work in two dimensions reported in
Chapter IV indicated that a depth of four times the height of the bed
element was adequate to approach the straight streamlines associated
with free stream flow. To further check this, consider the flow around
a single sphere, where, from [29]
Ua3
p =  cos 0 + Ur cos 0 (52)
2r
This is the same case described in Section 2.4.6. The expression for
the tangential velocity component is shown below, along with its value
when = rT/2, or points immediately above the center of the sphere.
3
1 F =p (+U sin 0)(1 + a ) (53)
r 2r
For 9 = n/2,
3
7 = U(l + a) (54)
2r
Equation (54) expresses the horizontal velocity across the
8 = T/2 axis above the center of the sphere. In a free stream such
velocity would simply be U. Hence, the velocity at any distance r can
be compared with U as an indicator of how closely flow at that depth
approaches the free stream. Using (54), it can be seen that the
velocity shown there is only 1 per cent greater than U when r = 3.7a.
Thus, 4a would provide less than 1 per cent deviation.
If (52), representing one plane in an axisymmetric problem,
is converted to Cartesian coordinates x and z for y = 0, (55) is
obtained.
Ua x
a x3 + Ux (55)
Ss3/2
where S =(x2 + z2)(2)(3)
From this, expressions for the x and zvelocity components can be
found, along with the ratio of the latter to the former. If this is
done, and z is chosen as 4a, with x as +a, it is found that the
vertical velocity component there is less than 1 per cent of the
horizontal component. Thus, at this point also, y = 4a gives a good
approximation of free stream flow, with streamlines horizontal.
For these reasons as well as earlier studies, a flow space of depth
four times the hemisphere height was chosen for the solution space
for (51).
Choosing a depth for the flow space assumes streamlines at
that elevation to be horizontal. Thus, at any greater depth than this
the indicated flow lines above that level will also be horizontal, and
the solution along a lower boundary will remain the same. It might
also be noted that the choice of such a solution space eliminates any
consideration of free surface effects, such as wave resistance. Flow
of greater depth relative to particle size may then be superimposed
on the flow space solved for with the knowledge that the solution of
distribution of velocity and pressure along the lower boundary will
remain essentially the same. The error here is determined by how
closely the stream surface at the chosen elevation approaches a flat
plane. It should also be noted that the depth chosen would form an
even better approximation in the cases of more closely packed hemi
spheres being treated here, since the stream surfaces in these cases
more rapidly approach free stream conditions.
5.2.3 Boundary conditions.Due to the periodic nature of the
hemispheres in the solutions, numerous conditions of symmetry are avail
able which enable reduction of the solution to handling a typical
repeated portion of the total flow space. Figures 18 and 19 represent
the portions chosen for the solutions of this dissertation. It can
be seen that the two drawings are different in that Figure 18 could be
divided by "symmetry" one more time. For purposes of the finite
differences solution, however, this was not done. Two alternative
approaches were available, but neither seemed as practical as simply
including two hemispherical portions in the solution space. One of
these alternatives was assuming an equipotential surface between the
two hemispherical portions and making finite differences solutions in
which this assumed potential surface was to be checked and itself
adjusted and treated as a variable in a solution. However, this would
have involved a complete solution of the problem for each assumed sur
face and would have required, therefore, exorbitant computation time.
The other alternative was to choose a plane between the two hemispheres
and use it along with conditions of "symmetry" to reduce by onehalf
the number of points involved in the differences solution. In the case
of Figure 18, the proximity of the hemispheres introduces complexities
which offset the benefits gained in computation time. For the Chepil
case, however, the latter approach could be used easily 'as indicated
in Figure 19.
The constant xplanes passing through the hemispheres are
equipotential surfaces. All other boundaries in the solution spaces
represent surfaces across which there is no flow, with the exception
of the added symmetry boundary shown in Figure 19. This means that the
derivative of the potential function normal to the given surface is
equal to zero. Thus, the following conditions hold.
On constant yplanes through hemispheres:
a = 0 (56)
y
On upper zsurface:
= 0 (57)
az
xconstant plane:
p = p2 = constant
Figure 18. Solution space for closely packed hemispheres.
xconstant plane:
"foldedsymmetry" axis
\
z =4a 
xconstant plane:
cp= p = constant
y f i^^'
Solution space for Chepil's arrangement.
Figure 19.
On lower zsurface:
zp
7z 0
(58)
On hemispherical surfaces:
= 0 (59)
6h
The selection of the symmetry plane for the Chepil case is
based on a type of "folded" symmetry. This is indicated in Figure 20
and equations (510) and (511).
x +xo
1 0
y=0
Constant 2plane
Figure 20. Foldedsymmetry boundary.
K'P'
y = yl
y0
Yo
'PB
YO
x Xo
1 0
PA 'P = ( B P2) (510)
PB = (2 + 1) A (511)
The problem at hand is therefore to solve Laplace's equation
in the regions indicated,subject to the foregoing boundary conditions.
5.3 General finite differences approach
The method of finite differences consists of replacing a function
continuous within a region by its values at certain points in that
region. The equation to be solved is then expressed in difference form
at each of these points, and the problem becomes one of solving a set
of simultaneous equations corresponding to the points chosen. As noted
earlier, one of the reasons for choosing this method for the present
work is that, due to the nature of the elliptic equations, ". . there
is no problem with stability of difference equations approximating
elliptic partial differential equations" [59].
Some consideration was given to the type of coordinate system
to be used. The spherical coordinate system had some advantages in
expressing the noflow condition at one hemispherical boundary, but was
less efficient at other boundaries and at the second hemisphere. Hence,
the ordinary rectangular Cartesian coordinatesx, y, and zwere chosen.
A general sevenpoint scheme was chosen for the difference solution, as
shown in Figure 21, with equal increments of length h in all three
coordinate directions. For this scheme, the difference equation becomes
that of (512), as found in works by Allen [603, Collatz [61), and
Allen and Dennis [62]. Section 5.4 will give details of how this
is obtained.
E (i 6 cp = 0
i=l
6
(p = Cp
S6 i=l i
L=1
(512)
(513)
6
where E pi represents the 6 points around p .
i=l
Figure 21.
Sevenpoint finite difference scheme.
Similarly, if a net, or lattice, were placed over the entire
solution space, an equation could be written for each point of the net
of the same type as (512). Thus, Laplace's equation could be written
in the form of a matrix equation as
AX = B (514)
where A is the matrix of coefficients of cp, X is the matrix containing
p elements, and B is a matrix determined by the boundary conditions.
Therefore, inversion of the matrix A should enable a solution for X,
the potential field. However, A is usually a very large matrix which
has many zero entries. Fox [63, p. 185] states that, "There is no point
in ever evaluating an inverse Ai for the purpose of solving equations
of the form AX = B. Elimination and back substitution or its compact
equivalents are always faster. . On the subject of sparse matrices
(many zero entries) such as occur here, Fox [63, p. 189] says that,
"Iterative methods are used . for matrices of large order but with
many zero coefficients. .. "
The foregoing comments bring to mind the methods often termed
relaxation methods. These procedures were popularized before the rise
of computers by such people as Southwell [64], Shaw [65], and Allen (60].
As digital computers have become more available, variations in the
relaxation processes have led to expressions more applicable to the
computer. That method wherein the iterative methods are applied to
the points of the net singly and in an orderly repetitive fashion,
through numerous iterations, has become known as the Liebmann method
[66], which employs equation (513) when a rectangular Cartesian system
is used.
In order to accelerate convergence of the iterative method, the
process known as overrelaxation [67] is used, employing an overrelax
ation factor, w, and utilizing an equation of the form of (515).
6
cp0 = {i=1 + (1 u) c (515)
where cp' indicates the value from the preceding iteration. Obviously,
0
as in equation (513), the value of cp from (515) will show little
change as convergence is neared.
If the value of w is equal to one, the method is called the
Liebmann method, while the name extrapolated Liebmann is applied for
cases where w is not equal to unity. The latter case will be employed
here. Much work has been done, at least in twodimensional cases, in
finding an optimum w to give most rapid convergence [59]. Approximate
indications are that the optimum wvalue for the cases herein, lies
below about 1.80. However, no great amount of work will be done to
refine this value, as such work could easily involve very extensive
time.
The solutions of this dissertation will therefore be carried
out using a sevenpoint finite difference scheme in rectangular Cartesian
coordinates, applying the extrapolated Liebmann iterative method. The
ensuing sections will develop the needed relationships for use in
special situations.
5.4 Finite differences equations: interior space
This section will present the equations needed to handle all
the different situations which arise in the solutions. While many
applications of finite differences have been made in two dimensions,
few cases other than problems involving simple cubes, and the like,
have been treated in three dimensions. This condition has caused some
new procedures to be employed in the following work.
For clarity, certain conventions will be followed in presenting
the equations and their descriptions. First, each point will be given
three subscripts representing the three coordinate directions, with I
denoting the xdirection, J denoting the ydirection, and z represented
by K. Hence, the subscript I+1 implies the next point in the positive
xdirection beyond a point at I. Additionally, where one or more of
the variables is a constant for the investigations of a special relation
ship, the figure will be drawn in two dimensions, eliminating the third,
constant dimension.
5.4.1 General lattice point.Equation (513) was presented
as representing those cases where six adjacent points are available,
all at a distance h. However, frequently one or more points lie at
some distance other than h. For this reason, the general equation will
be developed for the case of Figure 22, where all six lengths are dif
ferent from h.
z
I,J,K+
y
chsh
I1,J,K 2h / 4h
S3h 1h PI+1,J,K
S6h
J1, ,JKK
I,J Kl
Figure.22. General lattice point.
Difference expressions can be formed for this case.
Forward difference:
S I+1+JK IJ,K
1h
Backward difference:
= II,J,K I1J,K
'2h
(516)
(517)
From these the following equation for the second derivative can be
obtained.
I+1J K I,J,K
1h
YI,J,K Y1,J2K
C2h
(1 +)h
(0l +0'2)h
SI+I,J,K IJK
I (1+ 2 )
SI,J,K I1,J,K
2 (l +2)
Similar expressions can be found for the other directions.
I ,J+1,K I,,K
S(PIJK+1 IJ),K
SQ'( +Q' )
5 5 6
SIJ,K ,J1,K
 I,J,K IJ,K1
6 (5 6)
(520)
(521)
Adding these three expressions yields an approximation for
the Laplacian.
(C) =
^ A &
u^ P
1 2 2
h
2 x 2
(518)
(519)
1h2 2
hr
2 6 2
12 2
2 2
1 2 2 I+1,J,K I,J,K IJK IlJK
h V p= ai(cI+a ) a (a +2)
2 1 +2 2 12
+ IJ+1K I,J,K I,JK I.,J1,K
%3(a%'') 3 (a%+a4)
+ I,J,K+l YIJ,K %IJ,K CI,J,K1 0 (522)
5 (5 6 6 5+6
Since the desire at a given point is to solve for the potential, or cp
value there, the equation above must be solved for cq(I,J,K).
( I+1,J, K +I1,J,K I,J+1,K
+ + +
CIJ = 1 1 1 (523)
al(al+ 2 2()l+(y 2a ) + 3(a'Y3+)
1 1 1
4(043+40 ) +5( 5+6) + 6(Ca5+016)
Note that if all legs are of length h, which implies that all crvalues
are 1.0, equation (523)reduces to (513), as it should.
Equation (523) is applicable for all interior points of the
solution space. The six points may contain among them points that
lie on a regular lattice point, on a hemispherical surface point where
p is being calculated, on or beyond a noflow boundary, on or beyond
the "foldedsymmetry" boundary of Figure 19, or on one of 'the two
equipotential planes at the xextremities of the solution space.
The next few sections will cover these cases as well as the grading of
the lattice. For clarity, the point cPIJ,,K will be termed the object
point, while the six surrounding points will be called adjacent points.
The case for an adjacent point, lying on a regular lattice point or on
a hemispherical surface point where c is being calculated, will be
omitted, as these simply involve the substitution of the present value
of cp for that adjacent point. The same is true when an adjacent point
lies on one of the equipotential planes.
5.4.2 Object point on planar noflow boundary.The term
planar is intended to exclude the hemispherical boundaries, which will
be treated in Section 5.5. Therefore, the planar surfaces indicated
are the constant ybounding planes and the constant zbounding planes,
as illustrated in Figures 18 and 19. If the object point lies on one
or two of such planes, one or more of the adjacent points will lie
outside the solution space, and it must be replaced by an equivalent
expression involving points on the interior of the space. For this
purpose, use is made of the fact that the derivative of cp normal to the
planes is zero. Then the value at such an exterior point can be
expressed as equaling the corresponding point inside the space. Fig
ure 23 illustrates what is meant. Here, as later, CB represents the
cpvalue at a surface point.
Constant zplane Constant zplane
I ,J+1,K
y=o 0'I,J,K y=0 I,J,K
I I
CP= I,J+1,K = B
(a) (b)
Figure 23. Examples for object point on planar boundary.
The preceding requirements can be summarized as follows, where
the equations imply replacing the left member by the right member for
use in equation (523).
For z = 0:
I,J,K = ,J,K+ (524)
For z on upper bound:
I,J,K+1 I=I,J,K1 (525)
For y = 0:
I,J1,K I,J+1,K (526)
For y = a (or y = 3a):
I,J+1,K = YI,J1,K (527)
Note that in the last two equations the terms may represent regular
lattice points, as in Figure 23a, or points on the hemispheres, as
in Figure 23b.
5.4.3 Object point on "foldedsymmetry" boundary.This
problem arises only in the case with hemispheres spaced wider apart,
as in Figure 7. One of the points, in the sevenpoint difference
scheme lies beyond the boundary, but the symmetry of the flow pattern
makes it possible to express this value at this point by a value inside
the given space. Equations (510) and (511) show the expressions
needed. The value thus obtained for this external point is then
inserted into its proper place in the sevenpoint formula.
5.4.4 Adjacent point on hemispherical surface.In some
instances, the lattice point being treated may have as adjacent points
one, two, or as many as three points which lie on a hemispherical bound
ary. At first, the author hoped to expand some of these points in
terms of other values in the field, while still applying the normal
derivative boundary condition (59) to most surface points. It was
found, however, that this approach created stability problems in some
portions of the iterative solutions. Therefore, the application of
the difference equations for (59) was extended to obtain a value for
all adjacent points falling on a hemisphere, with a small number of
exceptions. Discussion of the noflow condition and its use will occur
in ensuing sections.
The single exception occurs when both adjacent ypoints lie on
hemispherical surfaces. In order to obtain values for these two points,
a process which will be called arc'interpolation will be used. This
entails interpolating a value for the point from two surface points
where a value exists. The two points selected lie in the same zplane,
and the interpolation is made by the arc length between the three
points. These lengths along the arc are proportional to the angles
shown in Figure 24, and the ratio of lengths can be replaced by the
ratio of the angles, written in (528).
Constant zplane
ch
Il,J,K 0 I"
S3h
Point 0: cpI,
iPB
PI1,JI,K.
F2
Figure 24. Arc interpolation for surface value.
PIlJK,J1,K I1,J,K
1 2
(528)
A similar expression can be written for the point on the other
hemisphere, and the two values are then inserted into the general
sevenpoint equation.
It would have been possible to use the difference expressions
at these ypoints. It was felt, however, that in the regions where
both points lie on a hemisphere, the accuracy of the arc interpolation
method is of the same order as the expressions which would have to be
used for (59) in such cases. Therefore, arc interpolation was chosen
as a more convenient manner of obtaining values for the points which
were still consistent with the remainder of the potential field.
The evaluation of cp at the six adjacent points has been discussed,
with the exception of surface points to which the boundary condition
equation is applied, which will be covered in Section 5.5 and following.
5.4.5 Graded lattice.The greatest accuracy is desired near
the boundaries to enable velocity calculations on those surfaces. Those
points which are farther away from the hemispheres can be handled,
therefore, with a larger mesh spacing, h. This is more true the
nearer the approach to a free stream condition. An advantage is gained
in computatio, as using a larger lattice cube at some distance can
greatly reduce the number of points at which a solution is sought.
In the present cases, it was decided to use a larger spacing in the
upper onehalf of the solution space. Allen and Dennis [68] used a
scheme similar to that which is shown in Figure 25 and which will be
used here.
First, consideration will be given to those points where the
distances involved are all h or 2h (ac = 2 = 1). The problem can be
considered on three zplanes. First is the plane CDTH. Here, letters
will be used for points rather than the usual convention, as it may
allow greater clarity. For points such as C, D,. H, or J, six adjacent
points exist at a distance 2h; for points like Q, R, S, T, and U, the
six needed points lie at a distance h. Similarly, for all points in
plane EFKL, there are six points 2h away.
Q
S
U
J^
C
.4
R
 .
p
I 1
E
M
1* I
K
Figure 25. Grading the lattice.
D
h
y
1
'
I
I

There are two types of points in the zplane intermediate to
the above two; points such as V and those such as G. For points like V,
the following expression can be quickly written:
2e = 'A + 29V (529)
ox2 h2
The other four needed points lie in plane CDEF and can be utilized,
since the Laplacian is invariant with respect to a rotation of axes,
and thus
.y . + 2D E *F V (530)
by2 az2 (/2 h)2
Summing the two preceding equations leads to the final difference
expression.
2(cpA+pG) + C + D + E+ PF 8yPV = 0 (53i)
For points such as G, a triple McLaurin series expansion [68]
enables a result as below, with neglect of terms 0(h ).
C + D + E + cF + H + TJ + (PK + L. 8(G = 0 (532)
Due to the irrational numbers (square root of three) arising
from the hexagonal patterns of the hemispherical elements,.the succeed
ing planes of constant x treated in the grading process are not always
equidistant. This causes littleproblem for those points where
6 adjacent points are used, as these are treated precisely as indicated
by equation (523), the general equation. However, in the two cases
in the intermediate plane, equations (531) and (532) must be altered.
In the case of points such as V, rewrite (529) for the case where G
is at a distance h and A at a point a2h away.
2h22 c2(l)A +G l  4 (533)
6K2 L2(1+0 ) 2 2
Adding this to (530) and equating to zero produces the desired result.
C D + } 4 1 + p = 0 (534)
Points such as G require more attention. The triple McLaurin
expansion used to obtain (532) is such that, due to symmetry, all
oddorder terms (first derivative, third derivative, and so on) cancel
between the terms. However, an unequal xspacing causes those terms
involving a 6x to remain in the equation, though those oddorder terms
in the other directions still disappear. Here, the 0(h3) terms will
be neglected so that of the following form can be written for cp(x,y,z)
at some point 6x, 6y, and 6z from co.
cp(x,y,z) = Cp + 6x L + 6y . + 6z 6C
2 2 2
(1x)12 J + (6y)2 + 0(z)2 }
x by z
2 2
+ 6x 6y 2! + 8x 6z 62 + 6y 6z 2 (
b bay axaz ayaz
(535)
where all derivatives are taken at point 0. Expansions similar to this
are discussed by Olmsted [69].
The crossderivative terms cancel due to symmetry in directions
other than x, as do the first derivatives except for the xderivative,
which can be approximated by
V1 P (536)
(1+02 )h x G
and the second xderivative, which can be approximated by
1 2 2PV 0p q0G
h 2 =+ ( + (537)
x2 G 0(1+2) 1+2 2
Then (535) yields, when applied successively to all eight of the
adjacent points used by object point G, the following:
~ 2
i 8pG 4h2 (V22 G 2(2 l)h2 62cp
bx
4(a21)
I+ 12 V cP = 0 (538)
where cpi indicates the sum of the eight adjacent points. Setting
the Laplacian equal to zero yields an expression the same as (532),
except for including the last terms in the equation above.
One further point, such as N, is slightly different from V,
though the ideas are similar.
h2 G N1 2N (539)
hy
where (Nl) is the point lying a distance h from N and not shown
in Figure 25.
In plane DEKJ:
2 2 YD + + y + 4N
S, D E K J (540)
ax2 z 2 (,r2 h)2
The resulting difference equation is
2(CG+PN1) + D + CE + CPK + J 8N = 0 (541)
Treatment of points such as N becomes somewhat different when the
xintervals are not equal. A McLaurin expansion will also be applied
in this instance, and,as for point G, the first and second derivatives
will respect to x remain in the difference equation, as shown below.
.2( PG +N) + + (PE + K + J 8N
2
2h(O1)  h2("21) 2 = 0 (542)
2 6x 2 2
x
Evaluation of these two derivatives requires values at a point where
no value is being computed, as seen in Figure 25. Thus, such a
value is obtained by a linear interpolation between two known points,
such as D and E.
The foregoing completes the discussion on graded nets and some
special problems they create. Application of the grading process will
take place away from the hemispherical elements, as a finer lattice
is desired there.
5.5 Finite difference equations: hemispherical boundary
5.5.1 General.The only boundary condition which has not yet
been discussed is the normal derivative condition on the hemispherical
surfaces and expressed in equation (59). The surfaces of the two
hemispherical portions being treated can be expressed by
2 2 2
x + y + z = 1 (543)
(xv3)2 + (y1)2 + z = 1 (544)
Both of these are for spheres of unit radius. Since the gradient of
a scalar function indicates the normal to a surface represented by
that function, the desired derivative is found by taking the gradient
of the above surfaces, yielding
x + y + z + = 0 (545)
ax 6y az
(x) J + (y1) a + z 0 (546)
ax 6y dz
In order to assure conformation to the conditions of the
problem, (545) and (546) must be incorporated into the differences
solution. The surface points encountered will be thus treated, the
only exceptions being those few points mentioned in 5.4.4. The aim
will be to express the first derivatives at the surfaces in difference
form and use.them in equations (545) and (546). The ensuing sections
will describe the formation of these terms.
5.5.2 Xdirection derivative.In this section the general
means for expressing the boundary derivatives will be indicated.
Consider the situation of Figure 26.
PBI CP
II ,,J,K ^I+I,J,K
Constant zplane 1
\I+2,J,K
y
Figure 26. Xderivative condition at boundary.
Indicated in Figure 26 is afictitious cp value located outside the
solution volume. This point will be utilized in the derivative and
eliminated by expressing it in terms of points in the field. The
difference expression for the derivative at the boundary can be written
as (547).
pI Spher
Sphere 1: =4
Sphere 2:
2 h
ax 2a 2h
(547)
A Taylor's expansion can be written for the fictitious point, expand
ing it about the point (I+1). Neglecting terms of order h3 yields
the following.
Sphere 1:
SI+2 B
S= 9I 2 {rh
I+C1 1 (1+al1)h J
1
2h2 1
l 1 I+2 B/c +11
2: 2 L 1 (548)
(1 + 1
Then, from (547) and (548), the derivative is approximated by
12 I
23 1
h B + 2y,+1 + 2i I1+2 (549)
ax 1+011I+1 1+P+ 1+2
1 1
Similarly, on Sphere 2, the fictitious point can be found by expand
ing about the point (11), yielding an equation like (548). Subse
quent substitution into (547) and rearrangement yield
Sphere 2:
3 2 21
=h 7O 0 +B 3 2 21c (550)
+2 B +Q I2 I
Not in all cases, however, do two lattice points exist in a
direction away from the boundary, thus causing a need for some other
means of finding the derivative. Two cases exist: first, the next
field point might lie on the next hemisphere, or it might lie on a
lattice point with the second field point lying on a sphere. In both
instances, the choice here is to use a simple linear expression involv
ing a oneway forward or backward difference rather than a centered
difference. These difference expressions simply involve taking the
difference in value between the two points and dividing it by the
distance between the two.
T+1 B
Sphere 1: h = +
alh
(551)
Sphere 2: = B Yh
5.5.3 Ydirection derivative.For determining the derivative
in the ydirection, the same ideas exist as for the xderivatives.
However, since the points needed for an expansion in the ydirection
generally do not lie at a lattice intersection, special means must be
used to develop these points. For an indication of the problem, see
Figure 27. For simplicity, call the needed points PJ and PJ1. Their
value will be discussed. First, a series expansion can be written for
the fictitious points shown.
Sphere 1:
SPJ1 cp 2 cp + PJ1 2PJ
p = PJ 2h ii 2h + 2 2 (552)
h
Rearranging yields the derivative
Sphere 1:
h 2h =C B  + 2PJ  (553)
Similar expansion of the fictitious point at the second sphere results
in (554).
zconstant
(a)
zconstant
(b)
Figure 27. Yderivatives.
A+1
A1
Sphere 2:
_'_PJ 3PB PJ1
h h 2 =  2PJ + T (554)
The values indicated by PJ and PJ1 must be established by
interpolation, since they do not, in general, lie at a lattice point.
The interpolation will utilize points in the lattice with expansion
to be made about the nearer of two surrounding points in the xdirection,
shown in Figure.27a as PI.,J+1,K and cPI+1,J+1,K( C and PA).
First, consider the case where Ph is less than 1 h. Here, PJ
o
can be expressed as follows, for Sphere 1.
Sphere 1:
B PA+1 'C
PJ = CPA B( +1. C
A A '+{ +
0
CA+1 PC 21 1
+ B2 A (555)
aI + X + B
o
where A: (I+1,J+1,K)
A+l: (I+2,J+1,K)
C: (I,J+1,K)
A similar expression can be written by the other point needed, PJ1,
involving points one increment in the ydirection. In addition, the
same process yields an expression for points located from the second
hemisphere, indicated in Figure 27b.
Sphere 2:
+ + A+
$ +0+1
0
+ B 0 (556)
a + + a
o
where A: (I1,J1,K)
Ai: (I2,J1,K)
C: (I,J1,K)
Again a very similar expansion provides a value for PJ1.
To gain accuracy the expansions used are to be developed about
the point nearest the desired point. In the case where a is less
o
than B, a slightly different approach becomes more convenient due to
the proximity of the spherical surface. For that reason, the point at
a distance a will be expressed by expanding about PJ (or PJ1), and
subsequently an equation obtained for the latter point. This is
illustrated below.
Sphere 1:
9A PJ PJ PC
P = Pa o + 7 (_
0 + +V
o
o % C
PJ"l 0 Cpc + { .CPo (558)
0 a +
Again, in the same manner, an equation can be obtained for the second
sphere.
(+1 =1 1 + (559)
o o o
Both equations (558) and (559) are the same expressions to
be used for PJ1 at either sphere, with cPC and cpA replaced by points
one increment further away from the spheres.
All the immediately preceding equations depend for their
expansions upon the existence of certain lattice points. In some
instances, points such as A+l (or Al) lie on a spherical boundary.
The above equations may be used, however, with a proper value of .
In the case where the points such as A fall on a boundary surface,
a linear interpolation between points A and C will be used. This is
a firstorder approximation, assuming a oneway derivative between
points A and C equal to the difference in their values divided by the
distance between them. The same approach will be used for either PJ
or PJ1 from either sphere where the conditions warrant.
The foregoing equations for cp/6y, (553) and (554), are
based on those cases where two points such as PJ and PJ1 actually exist.
In some cases, they will not both be available, and these cases must be
treated in a different manner. Figure 28 demonstrates the two cases
to be met, shown with reference to obtaining the derivative on the
first sphere.
Constant zplanes
y
\I ,B 2
I CDB 2=PJ
J+1
CPB1
(a) (b)
Figure 28.
Special points for yderivative.
The case in Figure 28a will be treated by first evaluating PJ
by a linear interpolation between the two points (I,J + 1,K) and
(I+1,J+1,K). Then, the value of PJ is used in a forward difference
expression to form
h a PJ OB1
h = 
(560)
In the situation of Figure 28b, the same expression is used,
with PJ being determined by the process of arc interpolation, described
in an earlier section.
One further instance where a special case arises, involves
points near, but not at, y = 0 or y = 1, the two planar noflow bound
aries. Here, the value of PJ1 may correspond to a point on another
hemisphere outside the solution space, such as shown in Figure 23b.
As in that case, PJ1 will equal CB. The equations involving these
SJ+2
+
L
