Title Page
 Table of Contents
 List of Tables
 List of Figures
 List of Symbols
 Related background on lift force...
 Logarithmic velocity distribut...
 Two-dimensional work: earlier...
 Three-dimensional numerical...
 Results and comparisons
 Conclusions and future work
 Biographical sketch

Title: Hydrodynamic lift in sediment transport.
Full Citation
Permanent Link: http://ufdc.ufl.edu/UF00090227/00001
 Material Information
Title: Hydrodynamic lift in sediment transport.
Series Title: Hydrodynamic lift in sediment transport.
Physical Description: Book
Creator: Benedict, Barry A.
 Record Information
Bibliographic ID: UF00090227
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: alephbibnum - 000953375
oclc - 16919111

Table of Contents
    Title Page
        Page i
        Page ii
    Table of Contents
        Page iii
        Page iv
        Page v
    List of Tables
        Page vi
    List of Figures
        Page vii
        Page viii
        Page ix
    List of Symbols
        Page x
        Page xi
        Page xii
        Page xiii
        Page xiv
        Page 1
        Page 2
        Page 3
    Related background on lift force studies
        Page 4
        Page 5
        Page 6
        Page 7
        Page 8
        Page 9
        Page 10
        Page 11
        Page 12
        Page 13
        Page 14
        Page 15
        Page 16
        Page 17
        Page 18
        Page 19
        Page 20
        Page 21
        Page 22
        Page 23
        Page 24
    Logarithmic velocity distribution
        Page 25
        Page 26
        Page 27
        Page 28
        Page 29
        Page 30
        Page 31
        Page 32
        Page 33
        Page 34
        Page 35
        Page 36
        Page 37
        Page 38
        Page 39
        Page 40
    Two-dimensional work: earlier results
        Page 41
        Page 42
        Page 43
        Page 44
        Page 45
        Page 46
        Page 47
        Page 48
        Page 49
        Page 50
        Page 51
    Three-dimensional numerical solutions
        Page 52
        Page 53
        Page 54
        Page 55
        Page 56
        Page 57
        Page 58
        Page 59
        Page 60
        Page 61
        Page 62
        Page 63
        Page 64
        Page 65
        Page 66
        Page 67
        Page 68
        Page 69
        Page 70
        Page 71
        Page 72
        Page 73
        Page 74
        Page 75
        Page 76
        Page 77
        Page 78
        Page 79
        Page 80
        Page 81
        Page 82
        Page 83
        Page 84
        Page 85
        Page 86
        Page 87
        Page 88
        Page 89
        Page 90
        Page 91
        Page 92
        Page 93
        Page 94
        Page 95
        Page 96
        Page 97
        Page 98
        Page 99
        Page 100
    Results and comparisons
        Page 101
        Page 102
        Page 103
        Page 104
        Page 105
        Page 106
        Page 107
        Page 108
        Page 109
        Page 110
        Page 111
        Page 112
        Page 113
        Page 114
        Page 115
        Page 116
        Page 117
        Page 118
        Page 119
        Page 120
        Page 121
        Page 122
        Page 123
        Page 124
        Page 125
        Page 126
        Page 127
        Page 128
        Page 129
        Page 130
        Page 131
        Page 132
        Page 133
        Page 134
        Page 135
        Page 136
        Page 137
        Page 138
        Page 139
        Page 140
        Page 141
        Page 142
        Page 143
        Page 144
        Page 145
    Conclusions and future work
        Page 146
        Page 147
        Page 148
        Page 149
        Page 150
        Page 151
        Page 152
        Page 153
        Page 154
        Page 155
        Page 156
        Page 157
        Page 158
        Page 159
        Page 160
        Page 161
        Page 162
        Page 163
        Page 164
        Page 165
        Page 166
        Page 167
        Page 168
        Page 169
        Page 170
        Page 171
        Page 172
        Page 173
        Page 174
        Page 175
        Page 176
        Page 177
        Page 178
        Page 179
        Page 180
        Page 181
        Page 182
        Page 183
        Page 184
        Page 185
        Page 186
        Page 187
        Page 188
        Page 189
        Page 190
        Page 191
        Page 192
        Page 193
        Page 194
        Page 195
        Page 196
        Page 197
        Page 198
        Page 199
        Page 200
        Page 201
        Page 202
        Page 203
        Page 204
        Page 205
        Page 206
    Biographical sketch
        Page 207
        Page 208
Full Text







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


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.


ACKNOWLEDGMENTS . . . . . . . . . . .... ii

LIST OF TABLES . . . . . . . . . . . . vi

LIST OF FIGURES . .. ......... . . . . . .. vii

LIST OF SYMBOLS . . . . ... .. . . . . .. x

ABSTRACT ...... .. ........... .. ........ xiii

I INTRODUCTION. .. . . . . . . . . ..1


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


3.1 Development of logarithmic velocity
distribution . . . . .. ... .. 25
3.2 Special problems of present expressions . 26


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


4.1 Shapes studied . . . . . . .. 41
4.2 General methods--velocity 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


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 no-flow boundary 68
5.4.3 Object point on "folded-symmetry"
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 X-direction derivative . . . .. 79



5.5.3 Y-direction derivative . . . .. 81
5.5.4 Z-direction derivative . . . ... 87
5.5.5 Adjacent z-points and y-points 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 Einstein-El 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


APPENDIX . . . . . . . . . .. . . 149



REFERENCES . . . . . . . . . . . . . 201

BIOGRAPHICAL SKETCH . . . ... . . .. . . . 207


Table Page

1 Comparison of Proposed and Former Distributions ... .. 32

2 Chepil's Experimental Data . . . . . . ... 112

3 CL for Chepil's Work .................. 114

4 Roughness-Grain Size Ratios . . . . . . ... 115

5 Comparison of Theoretical with Chepil's Work ...... 123


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
Einstein-El 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 two-dimensional grain ..... 43

12 Piezometric head distribution on two-dimensional 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




20 Folded-symmetry boundary . . . . .

21 Seven-point finite difference scheme . .

22 General lattice point . . . . . .

23 Examples for object point on planar boundary

24 Arc interpolation for surface value . . .

Grading the lattice . . . . .

X-derivative condition at boundary .

Y-derivatives . . . . . .

Special points for y-derivative .

Z-point 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 . . . . . . . .









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

. .

* .



. .


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 xy-plane ... .. ..... .. . 136

53 Closely packed hemispheres: flow pattern on surface
viewed toward yz-plane . . .. .. . . . .... 137

54 Velocity suppression near sidewall. ...... .. 140


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 xy-plane

dA1 elemental surface area on hemisphere

g acceleration due to gravity

h full lattice increment

I point index, x-direction

J point index, y-direction

K point index, z-direction

k equivalent sand roughness

K diameter of hemispheres
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 x-direction

PI1 field point in x-direction

PJ field point in y-direction

PJI field point in y-direction

PK field point in z-direction

PK1 field point in z-direction

q total velocity at a point

qb velocity along base of the grain
r spherical coordinate

r roughness-grain 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
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 z-direction

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


Abstract of Dissertation Presented to the Graduate Council
in Partial Fulfillment of the Requirements for the
Degree of Doctor of Philosophy



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.


The analysis begins with a solution of the potential flow over

the two sets of hemispheres, solving here the three-dimensional 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.




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 two-dimensional

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 two-dimensional 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.



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. 64-65].

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 one-half the drag for experiments

with a Reynolds number based on the pipe diameter and the mean veloc-

ity in the range 360-1115 (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 well-known 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.



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 source-sink 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 Milne-Thomson [27]. He attributes the devia-

tions from theory primarily to generation of vorticity along the body.


0.8 -- Measured
-_- Potential Flow
P-Po i I
p v /2 1
o 0.2


0.4 x
0 0.2 0.4 0.6 0.8 1.0

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 two-dimensional study,

as shown in Figure 4.

Uniform Flow

a Long Circular

Figure 4. Jeffreys' cylinder.

Jeffreys developed his work based on the complex potential for

this case,

W = raU coth (--) (2-1)

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 (2-2)
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 (2-3)
1.43 p

Jeffreys noted that J. S. Owens [28] had earlier measured

the velocity required to move pebbles, finding

U2 = 1.65 ga (2-4)

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


Differences from theory are primarily due to two causes.

First, the motion measured was rolling rather than totally lifting,

which (2-5) 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 U-values 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

S= cos e + Ur cos e (2-6)

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 (2-7)

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 (2-8)

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 (2-8) 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 (2-8)

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 (2-9)
v 2 4

where dA = 2a2 sin28 dG


S= a p 2- (2-10)

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
3 u
3 t
F = --- (2-11)
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 (2-12)
1 2

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 (2-10) 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.

- 0
1 2

0 60 120


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 (2-8) and subtracting it from the force of (2-11) to give a

resultant vertical lift. Integration here occurs in a direction normal

to that in (2-9). 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 (2-13)

The lift per unit area is therefore

F -F u
S vvv vb t1 t
S= ,a = p --- (2-P.)
u 2 2

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


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 one-half 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

p = CL P (2-15)

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 natura-l 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 (2-16)

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 (2-16)

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 been-aimed 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 three-dimensional 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 two-dimensional 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




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 so-called Reynolds

stresses, after 0. Reynolds [46]. The final expression can be written

as below.

V = B + ln (3-1)
v f k

where B = C + In k

k = equivalent sand roughness

The so-called 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 B-values 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 (3-2)
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 (3-2).

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 (3-3)
s b

1 3
Volume of particles =2 s k (3-4)
2 s

Equating these two and solving yields


yb = 0.302 ks


1 Flow

Theoretical Bed
y = 0 in equation (3-7)

0.2 ks by Einstein
and El Samni


Figure 6. Arrangement of spheres and theoretical bed in
Einstein-El 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 (3-6)
s b (3-)

1 3
Volume of particles = K (3-7)
6 s

Equating these and solving yields

Y = 0.0335 K (3-8)
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 (3-2) 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.




I '1






Theoretical Bed


Figure 7.

Arrangement of spheres and theoretical bed
in Chepil's work.


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) (3-9)
vf k

Calculation will reveal one very desirable feature of (3-9),

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 (3-2). This is correct, since (3-2) 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 (3-2), other features can be noted by looking at

the change of velocity with distance from the wall, or dv/dy shown


dv f
Former: = 2.5
dy y (3-10)

For y 0 dv/dy -

y O dv/dy 0

dv 1
Proposed: = 2.5 v ] (3-11)
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

y-values, 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 (3-9). Table 1 shows the difference in velocity indi-

cated'by (3-2) and (3-9) 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 (3-9) is

dissipated very rapidly.



Variation between Force
(3-2) and (3-9) 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 (3-2) or (3-9), 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
qd with velocity v .

k, the equivalent sand roughness

vf, the friction velocity

From these measurements the following can be written, using first

expression (3-2).

Depth 1:

-2- = 8.48 + 2.5 In Md (3-12)
vf k

Depth 2:

= 8.48 + 2.5 In q- (3-13)
vf k

Combination of these equations and elimination of vf yield

29.7 pd (3-14)
(v )/(v v )
( p '

Elimination of k yields

V. v
vf = (3-15)
2.5 In (P)

Frequently the depths chosen are for p = 0.90 and q = 0.15.

Use of the proposed expression in (3-9) for equations such as

(3-12) and (3-13) 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

(3-14) and (3-15).

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. 276-277] 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.


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




Probable True


Sketches for sidewall effect.

Figure 8.

The following expressions can be written relating the shear


1.0 1
= -- = tan f (3-16)
2.0 Y2

T1.0 h + T2.0(W h tan a) = yh WSb (3-17)

This can be written as

T1.0 + tan ') = WSb (3-18)

Equation (3-16) results from assuming 71 and 72 to both be

proportional to.their respective y-values 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 (3-2)

will be used for the velocity.

[8.48 + 2.5 In [ ] -A = [8.48 + 2.5 In [-]} T (3-19)
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 (3-20)
1 2



'I = y/k

Evaluation of constant for sidewall analysis.


Figure 9.

Introducing y/k = T enables the writing of

,/I F(I,) = /F(T-2 ) (3-21)


F(C) = [In (29.7 1)] 7 /T

The curve shown in Figure 9 indicates that the F-function can

be approximated by a function of the form

F(J) = alb (3-22)

This would enable the writing of

1 27 (y_1/k) y_ k2
S_ 2 =(3-23)
Sy2/k2) 2 1

Using the fact that = tan C, it can be written that

tan a = (-) (3-24)

Evaluation from Figure 9 reveals that bcequals about 0.664. Thus,

tan = 1) (3-25)

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.


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.



4.1 Shapes studied

In earlier work the author has studied some two-dimensional 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 finite-difference 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 methods--velocity 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



/ Stagnant

/ b/d = 0.561


Figure 10. Grains placed in rough bed configuration.


Figure 11.

Y bottom


y top

Pressure distribution on two-dimensional grain.


(z + =bo m
y bottom y


(z + -)
y top

Figure 12. Piezometric head distribution on two-dimensional grain.



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


Results for velocity-squared 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 (4-1)
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-

to-length 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






5.0 \3.0 \ 2.0 \

0 1.0
Dashed Line: cosh bed shape

Figure 13. Distributions of surface hydrodynamic pressure

0.5 1.0

Figure 14. Lift coefficient CL.

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 (3-9)

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 (3-9) at the top of the grain, evaluate ut

from the measured data.

3. Calculate the lift by equation (4-2).

Unit area: X = CL p (4-2)

Thus, by relating the velocity profiles, it is possible to

make a prediction of the lift.



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 two-dimensional 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,



I DLviding Streamline ,
......... I d

Actual Bed

Figure 16. Definition sketch for experimental application.


Theory - -
1 2
P um \\

0.4 --------


0 0.2 0.4 0.6 0.8 1.0
Pr: reference pressure Bed Location
u : mean velocity

Figure 17. Comparison of theoretical and measured values.



5.1 General

Previously it was indicated that irrotational flow theory would

be utilized to study velocities on a three-dimensional bed surface

toward the end of making an analytical evaluation of hydrodynamic lift

on such a surface. The results of some previous similar two-dimensional

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 (5-1) in the given flow space.

2 2 2
2+ -2 + -2 = 0 (5-1)
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 (5-1) is expressed instead in integral form and

a solution made,was eliminated because some past works indicated that

three-dimensional 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 two-dimensional

results are then inverted back to the three-dimensional 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


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


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]

p = -- cos 0 + Ur cos 0 (5-2)

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.

1 F =p (+U sin 0)(1 + a ) (5-3)
r 2r

For 9 = n/2,

7 = U(l + a-) (5-4)

Equation (5-4) 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 (5-4), 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 (5-2), representing one plane in an axisymmetric problem,

is converted to Cartesian coordinates x and z for y = 0, (5-5) is


Ua x
a x3 + Ux (5-5)

where S =(x2 + z2)(2)(3)

From this, expressions for the x- and z-velocity 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 (5-1).

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 one-half

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 x-planes 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 y-planes through hemispheres:

a = 0 (5-6)

On upper z-surface:

= 0 (5-7)

x-constant plane:
p = p2 = constant

Figure 18. Solution space for closely packed hemispheres.

x-constant plane:
"folded-symmetry" axis


z =4a -

x-constant plane:
cp= p = constant

y f i^^'

Solution space for Chepil's arrangement.

Figure 19.

On lower z-surface:

7z- 0


On hemispherical surfaces:

= 0 (5-9)

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 (5-10) and (5-11).

x +xo
1 0


Constant 2-plane

Figure 20. Folded-symmetry boundary.


y = yl



x -Xo
1 0

PA 'P = ( B P2) (5-10)

PB = (2 + 1) A (5-11)

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 no-flow condition at one hemispherical boundary, but was

less efficient at other boundaries and at the second hemisphere. Hence,

the ordinary rectangular Cartesian coordinates--x, y, and z--were chosen.

A general seven-point 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 (5-12), 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

(p = Cp
S6 i=l i



where E pi represents the 6 points around p .

Figure 21.

Seven-point 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 (5-12). Thus, Laplace's equation could be written

in the form of a matrix equation as

AX = B (5-14)

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 A-i 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 (5-13) 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 (5-15).


cp0 = {i=1 + (1 u) c (5-15)

where cp' indicates the value from the preceding iteration. Obviously,
as in equation (5-13), the value of cp from (5-15) 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 two-dimensional cases, in

finding an optimum w to give most rapid convergence [59]. Approximate

indications are that the optimum w-value 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


The solutions of this dissertation will therefore be carried

out using a seven-point 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 x-direction, J denoting the y-direction, and z represented

by K. Hence, the subscript I+1 implies the next point in the positive

x-direction 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 (5-13) 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.



I-1,J,K 2h / 4h

S3h 1h PI+1,J,K
J-1, ,JKK-
I,J K-l

Figure.22. General lattice point.

Difference expressions can be formed for this case.

Forward difference:

S I+1+JK- IJ,K

Backward difference:

= II,J,K I-1J,K



From these the following equation for the second derivative can be


I+1J K I,J,K

YI,J,K Y-1,J2K

(1 +)h
-(0l +0'2)h

I (1+ 2 )

SI,J,K I-1,J,K
2 (l +2)

Similar expressions can be found for the other directions.

I ,J+1,K I,,K

SQ'( +Q' )
5 5 6

SIJ,K ,J-1,K

- I,J,K IJ,K-1
6 (5 6)



Adding these three expressions yields an approximation for

the Laplacian.

(C) =
-^ A &
u^ -P

1 2 2
2 x 2



1h2 2
2 6 2

12 2
2 2

1 2 2 I+1,J,K I,J,K IJK I-lJK
-h V p= ai(cI+a ) a- (a +2)
2 1 +2 2 12

+ IJ+1K I,J,K I,JK I.,J-1,K
%3(a%'') 3 (a%+a4)

+ I,J,K+l YIJ,K %IJ,K CI,J,K-1 0 (5-22)
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 +I-1,J,K I,J+1,K

+ + +

CIJ =- 1 1 1 (5-23)

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 cr-values

are 1.0, equation (5-23)reduces to (5-13), as it should.

Equation (5-23) 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 no-flow boundary, on or beyond

the "folded-symmetry" boundary of Figure 19, or on one of 'the two

equipotential planes at the x-extremities 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 no-flow 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 y-bounding planes and the constant z-bounding 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

cp-value at a surface point.

Constant z-plane Constant z-plane

I ,J+1,K

y=o 0'I,J,K y=0 I,J,K

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 (5-23).

For z = 0:

I,J,K- = ,J,K+ (5-24)

For z on upper bound:

I,J,K+1 I=I,J,K-1 (5-25)

For y = 0:

I,J-1,K I,J+1,K (5-26)

For y = a (or y = 3a):

I,J+1,K = YI,J-1,K (5-27)

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 "folded-symmetry" boundary.--This

problem arises only in the case with hemispheres spaced wider apart,

as in Figure 7. One of the points, in the seven-point 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 (5-10) and (5-11) show the expressions

needed. The value thus obtained for this external point is then

inserted into its proper place in the seven-point 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 (5-9) 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 (5-9) was extended to obtain a value for

all adjacent points falling on a hemisphere, with a small number of

exceptions. Discussion of the no-flow condition and its use will occur

in ensuing sections.

The single exception occurs when both adjacent y-points 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 z-plane,

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 (5-28).

Constant z-plane

I-l,J,K 0 I"
Point 0: cpI,



Figure 24. Arc interpolation for surface value.

PIlJK,J-1,K I-1,J,K
1 2


A similar expression can be written for the point on the other

hemisphere, and the two values are then inserted into the general

seven-point equation.

It would have been possible to use the difference expressions

at these y-points. 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 (5-9) 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 one-half 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 z-planes. 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.







- -.


I- 1



1* I


Figure 25. Grading the lattice.









There are two types of points in the z-plane 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 (5-29)
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 (5-30)
by2 az2 (/2 h)2

Summing the two preceding equations leads to the final difference


2(cpA+pG) + C + D + E+ PF 8yPV = 0 (5-3i)

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 (5-32)

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 (5-23), the general equation. However, in the two cases

in the intermediate plane, equations (5-31) and (5-32) must be altered.

In the case of points such as V, rewrite (5-29) for the case where G

is at a distance h and A at a point a2h away.

2h22 c2(l)A +G l -- 4 (5-33)
6K2 L2(1+0 ) 2 2

Adding this to (5-30) and equating to zero produces the desired result.

-C D + } 4 1 + p = 0 (5-34)

Points such as G require more attention. The triple McLaurin

expansion used to obtain (5-32) is such that, due to symmetry, all

odd-order terms (first derivative, third derivative, and so on) cancel

between the terms. However, an unequal x-spacing causes those terms

involving a 6x to remain in the equation, though those odd-order 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


where all derivatives are taken at point 0. Expansions similar to this

are discussed by Olmsted [69].

The cross-derivative terms cancel due to symmetry in directions

other than x, as do the first derivatives except for the x-derivative,

which can be approximated by

V1 P (5-36)
(1+02 )h x G

and the second x-derivative, which can be approximated by

1 2 2PV 0p q0G
h 2 =+ -( + (5-37)
x2 G 0(1+2) 1+2 2

Then (5-35) 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
I+- 12 V cP = 0 (5-38)

where cpi indicates the sum of the eight adjacent points. Setting

the Laplacian equal to zero yields an expression the same as (5-32),

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 N-1 2N (5-39)

where (N-l) 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 (5-40)
ax2 z 2 (,r2- h)2

The resulting difference equation is

2(CG+PN-1) + D + CE + CPK + J 8N = 0 (5-41)

Treatment of points such as N becomes somewhat different when the

x-intervals 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

2h(O-1) -- h2("2-1) -2 = 0 (5-42)
2 6x 2 2

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 (5-9). The surfaces of the two

hemispherical portions being treated can be expressed by

2 2 2
x + y + z = 1 (5-43)

(x-v3)2 + (y-1)2 + z = 1 (5-44)

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 (5-45)
ax 6y az

(x-) -J + (y-1) a + z 0 (5-46)
ax 6y dz

In order to assure conformation to the conditions of the

problem, (5-45) and (5-46) 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 (5-45) and (5-46). The ensuing sections

will describe the formation of these terms.

5.5.2 X-direction derivative.--In this section the general

means for expressing the boundary derivatives will be indicated.

Consider the situation of Figure 26.

-II ,,J,K ^I+I,J,K

Constant z-plane 1


Figure 26. X-derivative condition at boundary.

Indicated in Figure 26 is a-fictitious 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 (5-47).
pI Spher
Sphere 1: =4

Sphere 2:

2 h
ax 2a 2h


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
2h2 1
l 1 I+2 B/c +11
2: 2 L 1 (5-48)
-(1 + 1

Then, from (5-47) and (5-48), the derivative is approximated by

1-2 I
23 1
h B- + 2y,+1 + -2i I1+2 (5-49)
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 (1-1), yielding an equation like (5-48). Subse-

quent substitution into (5-47) and rearrangement yield

Sphere 2:

3 2 2-1
=h -7O 0 +B 3 2 21c (5-50)
+2 B +Q I-2 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 one-way 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 = +

Sphere 2: = -B Yh

5.5.3 Y-direction derivative.--For determining the derivative

in the y-direction, the same ideas exist as for the x-derivatives.

However, since the points needed for an expansion in the y-direction

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 (5-52)

Rearranging yields the derivative

Sphere 1:

h 2h- =C B - + 2PJ - (5-53)

Similar expansion of the fictitious point at the second sphere results

in (5-54).





Figure 27. Y-derivatives.



Sphere 2:

_'_PJ 3PB PJ1
h h 2- -= - 2PJ + T (5-54)

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 x-direction,

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
can be expressed as follows, for Sphere 1.

Sphere 1:

B PA+1 'C
PJ = CPA B( +1. -C
A A- '+{ +

CA+1 PC 21 1

+ B2 A (5-55)
aI + X + B

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 y-direction. In addition, the

same process yields an expression for points located from the second

hemisphere, indicated in Figure 27b.

Sphere 2:

+ + A-+
$ +0+1

+ B 0 (5-56)
a + + a

where A: (I-1,J-1,K)

A-i: (I-2,J-1,K)

C: (I,J-1,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
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:

P = Pa -o + 7 (_
0 + +V

o % C
PJ"l- 0 -Cpc + { .CPo (5-58)
0 a +

Again, in the same manner, an equation can be obtained for the second


(+1 =1 1 + (5-59)
o o o

Both equations (5-58) and (5-59) 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 A-l) 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 first-order approximation, assuming a one-way 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, (5-53) and (5-54), 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 z-planes


\I ,B 2



(a) (b)

Figure 28.

Special points for y-derivative.

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 = --


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 no-flow 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



University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs