|
UFL/COEL-95/012
EVALUATION AND DEVELOPMENT OF PARTICLE
SALTATION MODELS
by
Michael R. Krecic
Thesis
1995
EVALUATION AND DEVELOPMENT OF PARTICLE SALTATION MODELS
By
MICHAEL R. KRECIC
A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF
FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE
DEGREE OF
MASTER OF SCIENCE
UNIVERSITY OF FLORIDA
1995
ACKNOWLEDGMENTS
I would like to especially thank my advisor and supervisory committee chairman, Dr.
Daniel M. Hanes, for his support during my master's research here at the University of
Florida. I would also like to thank the members of my committee, Dr. Robert J. Thieke
and Dr. Ashish J. Mehta, who have helped me through many research and class problems.
Also my thanks go to Dr. Renwei Mei for his help in understanding some of the force
terms.
Thanks are also extended to the staff of the Coastal and Oceanographic Engineering
Department and the Coastal and Oceanographic Engineering Laboratory. Special
acknowledgment is given to Becky, Lucy, Helen, John, and Vernon.
A special thank you is given to my fellow young people, Santi and Leah, Mike,
Marshall, Al, Kenny and Kim, Annette, Craig, Justin, Paul, Mark and Allison, Susan, Chris
and Monica, Mark, Yigong, Bill, Rob and Cindy, Darwin, Eric, Wayne, and Jie.
Finally I would like to thank my parents, Russ and Trudy, for supporting me as I
meander through life. They have always encouraged me to do anything I wanted. That's
probably a good thing since I have yet to decide what I want to do with my life. Also I
would like to thank my two brothers, Matt and Dan, for not ribbing me too hard for
extending my professional student status.
TABLE OF CONTENTS
Page
ACKNOW LEDGM ENTS ........................................................................................... ii
LIST OF TABLES ..................................................................................................... v
LIST OF FIGURES....................................................................................................... vi
LIST OF SYM BOLS..................................................................................................... ix
CHAPTER 1: INTRODUCTION .................................................. ........................... 1
CHAPTER 2: PREVIOUS W ORK ...............................................................................
M echanism s of Grain Saltation .................................................... .......................... 4
Incipient M otion..........................................................................................................4
Particle Saltation .................................................................................................. 6
Force Analysis ......................................................................................................... 8
Gravity Force ....................................................................................................... 8
Added M ass Force ............................................................................................... 9
Drag Force........................................................................................................... ... 10
Shear Lift Force ................................................................................................. 10
Basset History Force .......................................................... .............................. 12
M agnus Lift Force...................................................... .......................................... 13
Applicability of Forces For Different Reynolds Number Ranges.............................. 14
Predictive M odels...................................................................................................... 15
M urphy and Hooshiari (1982) .................................................. ........................ 15
van Rijn (1984) .................................................................................................. 17
W iberg and Sm ith (1985) ..................................................... ........................... 19
Nino and Garcia (1992).......................................................................................... 23
Lee and Hsu (1994).............................................................................................. 24
CHAPTER 3: SALTATION MODEL FORMULATION ............................................ 27
Single Grain M odel................................................................................................ 27
Equation of M option ............................................................................................ 27
Boundary Conditions.................................................... .................................... 31
M ethod of Solution ............................................................................................ 32
Sensitivity Analysis................................................................................................. 32
Random Trajectory M odel......................................................................................... 38
Determination of Initial Angles............................................................................... 40
B oundary C onditions.............................................................................................. 42
Sensitivity Analysis........... ........................................................................ ......... 43
CHAPTER 4: SALTATION MODEL RESULTS ........................................................45
Single G rain M odel................................................................................................ 45
W iberg and Sm ith (1985) Shear Lift....................................................................... 45
Saffm an (1965) Shear Lift ...................................................................................... 49
Adding M agnus Effect..................................................... ................................. 53
R andom Trajectory M odel......................................................................................... 55
Com prison to Previous D ata Set.......................................................................... 55
Comparison of Grain Ejection Methods.............................................. .............. 56
CHAPTER 5: CONCLUSIONS ................................................................................ 68
APPENDIX A: FORMULATION OF FORCE COMPONENTS.............................. 71
APPENDIX B: BASSET HISTORY FORCE EVALUATION......................................73
APPENDIX C: MATLAB PROGRAMS USED FOR MODEL DEVELOPMENT........ 75
R EFER EN CES............................................................................................. .............. 90
BIOGRAPHICAL SKETCH ....................................................................................... 93
LIST OF TABLES
Table Page
3.1 Distribution of Initial Angles ......................................................... ...................41
4.1 Average Saltation Characteristics For Fernandez Luque And van Beek (1976).........56
4.2 Average Saltation Characteristics From Random Trajectory Model........................ 64
4.3 Saltation Characteristics With Grain Ejected With Shear Flow ............................... 67
5.1 Summary of Forces In M odels ................................................ ......................... 69
LIST OF FIGURES
Figure Page
2.1 Critical Shields' Parameter As A Function of Archimedes Number (Sleath, 1984) .....5
2.2 Shear Lift Results From Velocity Gradient Across Grain.......................................... 10
2.3 Magnus Force Results From Velocity Gradient Across Grain And Viscous Effects.. 13
2.4 Comparison of Typical Trajectory and Murphy, et al. (1982) Model..................... 17
2.5 Observed Trajectory and Model Comparison (van Rijn, 1984).............................. 19
2.6 Sign of Form Lift for Different 'P1 Angles of a Body Relative to Flow.................... 21
2.7 Comparison of Predicted and Observed Grain Saltation Trajectory..........................22
2.8 Comparison Between Predicted and Observed Trajectory (Nino et al., 1992) ...........24
2.9 Trajectories With Same Angle and Different Velocity (Lee and Hsu, 1994) .............25
3.1 Force Definition Sketch of Grain Saltation......................................... .............. 29
3.2 D definition of B ed Slope.......................................................................................... 29
3.3 The Saltating Grain Initial Position.................................................................. ..... 32
3.4 Particle Trajectories With A Particle Size Of 0.18 Centimeters................................. 34
3.5 Particle Trajectories With A Particle Size Of 3.1 Centimeters................................. 35
3.6 Particle Trajectories Sensitivity to Particle Rotation ............................................ 36
3.7 Particle Trajectory Sensitivity to Initial Trajectory Angle........................................ 37
3.8 Particle Trajectory Sensitivity to Initial Velocity....................................................... 37
3.9 Comparison Of Variable Bed Roughness Values .............................................. 39
3.10 Simulated Trajectories of Particles With Different Densities ............................... 39
3.11 Typical Random Trajectories With Initial Velocity 5ucr........................................42
3.12 Particle Trajectories Of Same Initial Velocity And Angles Symmetric About The
Vertical Plane ......................................................................................................................43
4.1 Model Comparison With Fernandez Luque and van Beek (1976) ...........................46
4.2 Reynolds Number and Drag Coefficient As Calculated By The Model...................46
4.3 Sim ulated Particle V elocities.................................................................................. 47
4.4 Sim ulated Particle Accelerations ...........................................................................47
4.5 Sim ulated Particle Forces....................................................................................... 48
4.6 Simulated Particle Basset History Force...................................................................48
4.7 Model Comparison With Fernandez Luque and van Beek (1976) ........................... 50
4.8 Reynolds Number And Drag Coefficient As Calculated By Model............................ 51
4.9 Sim ulated Particle Velocities.................................................................................... 51
4.10 Simulated Particle Accelerations ................................................................. 52
4.11 Sim ulated Particle Forces..................................................................... .............. 52
4.12 Simulated Particle Basset History Force............................................................... 53
4.13 Simulated Particle Trajectory With Magnus Effect ................................................. 54
4.14 Simulated Particle Combined Lift Force Effect..................................................... 55
4.15 Simulated Trajectories With Grain Ejected With And Against Shear Flow.............. 57
4.16 Simulated Trajectories With Grains Ejected With Shear Flow............................. 57
4.17 Simulated Trajectories With Initial Velocity, uc................................................. 58
4.18 Simulated Trajectories With Initial Velocity, 2ucr.............................. .............. 59
4.19 Simulated Trajectories With Initial Velocity, 3ucr.............................. .............. 59
4.20 Simulated Trajectories With Initial Velocity, 4ucr............................................... 60
4.21 Simulated Trajectories With Initial Velocity, 5ur................................ ........... .. 60
4.22 Simulated Trajectories With Initial Velocity, 6ucr................................................. 61
4.23 Simulated Trajectories With Initial Velocity, 7ucr............................................61
4.24 Simulated Trajectories With Initial Velocity, 8ucr.............................. .............. 62
4.25 Simulated Trajectories With Initial Velocity, 9ucr............................................ 62
4.26 Simulated Trajectories With Initial Velocity, 10ucr.................................................. 63
4.27 Simulated Trajectories With Initial Velocity, 100ucr......................................... 63
A.1 Particle Relative Velocity ................................................................................. 71
LIST OF SYMBOLS
A cross-sectional area of particle
Ar Archimedes number
b proportionality constant
ci proportionality constant
c2 proportionality constant
CD drag coefficient
CM added mass coefficient
CL lift coefficient
d grain diameter
Fa added mass force
Fb Basset history force
Fd drag force
Fg submerged weight
FiMagnus Magnus lift force
Flshear shear lift force
g acceleration due to gravity
ge effective gravity
H saltation trajectory maximum height
I moment of inertia
ks equivalent bed roughness
L saltation trajectory maximum length
Lb saltation length for random model
AL added saltation length
Re Reynolds number
Re* grain Reynolds number
s specific gravity of particle
S surface area of top half of sphere
Ts time particle travels to complete saltation
u saltating particle longitudinal velocity
ucr critical velocity for incipient motion
Uf fluid velocity
Ufrop fluid velocity at top of particle
UfBot fluid velocity at bottom of particle
UATop relative fluid velocity at top of particle
UABot relative fluid velocity at bottom of particle
uo initial longitudinal particle velocity
u* fluid shear velocity
V volume of saltating particle
Vr relative saltating particle velocity
w saltating particle vertical velocity
wo initial particle vertical velocity
x longitudinal coordinate
z vertical coordinate
zo zero fluid velocity level
P bed angle
0 horizontal plane angle
0c nondimensional Shields' parameter
v fluid kinematic viscosity
tL fluid dynamic viscosity
Q particle angular velocity
no initial particle angular velocity
TI flow angle of attack
T2 related to airfoil centerline curvature
p density of fluid
ps density of saltating particle
Tb bed shear stress
Tcr critical bed shear stress
0 azimuth angle
01 arbitrary surface angle
(2 arbitrary surface angle
AO difference between 01 and 02
Abstract of Thesis Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Master of Science
EVALUATION AND DEVELOPMENT OF PARTICLE
SALTATION MODELS
By
Michael R. Krecic
August 1995
Chairman: Daniel M. Hanes
Major Department: Coastal and Oceanographic Engineering
The sediment transport mode of saltation is very important in developing bedload
transport theories. If saltating particle trajectories are simulated, a bedload concentration
profile may be obtained. Theoretically, this can be done by determining the amount of
time particles spend in any arbitrary layer above the theoretical bed. Then the
concentration profile may yield a rate of sediment transport. Presently there are
theoretical models used to predict the trajectories of particles.
The theoretical models have been supported by experimental observations.
Laboratory observations of the saltation process in adjustable flumes have been reported.
Saltation characteristics such as trajectory height, length, and longitudinal particle velocity
are measured through the use of high speed photography and standard video imaging
techniques. The data presently available range from saltating sand grains to gravel.
These data are used to evaluate the predicative capabilities of five different theoretical
models. Comparisons are performed by plotting the measured particle trajectories and the
predicted trajectories. Comparisons are also made by investigating the proper application
of forces in the predictive models. The models are able to predict trajectories fairly well
except they usually include too many adjustable parameters.
When examining the predictive models it is found that the shear lift force has a greater
effect on large particles than on small particles whereas the opposite is true for the Basset
history force. The Magnus effect is found to have a large impact on particle saltation.
The initial value of angular velocity is not well known and needs further evaluation. The
particles that are ejected against the flow and with the flow as compared to particles
ejected with the flow only have greater average saltation lengths, heights, and velocities by
approximately ten percent.
A new predictive model is proposed to relate experimental observations and theory.
This model borrows from those previous models as well as adding some new features.
The particle's angular velocity is continuously updated and particles may be ejected into
the shear flow. Sensitivity analyses are performed to determine the relative importance of
the parameters involved. The model is also compared to experimental observations.
It is also deemed necessary to measure actual particle angular velocities to more
appropriately represent the particle's trajectory. In addition, an experiment involving
particles ejecting into a shear flow needs to be conducted. If this future experiment
validates particle ejection into the flow, the shear lift force may need to be reformulated.
CHAPTER 1
INTRODUCTION
Sediment transport may occur in several different modes. These modes depend on
the bed shear stress and the characteristics of the sediment. They include rolling, saltation,
suspension, and sheet flow. With increasing shear stress the transport progresses as stated
above. Rolling is where the grains maintain contact with the bed and are driven by grain-
to-grain forces. Suspension occurs when the bed shear velocity exceeds the particle fall
velocity. This causes the grains to remain supported in the fluid column by turbulence and
other fluid forces. Sheet flow may occur if the shear stress is very high and results in a
layer of grains, usually 5-10 grain diameters thick, moving above the bed in unison. All of
these modes of transport may occur simultaneously if the shear stress is sufficiently large.
The transport mode saltation is reviewed herein. The term saltation was first used by
Gilbert (1914) and comes from the Latin word "saltare" meaning to leap or dance.
Saltation occurs when the shear stress is great enough to allow a grain to overcome
gravity but not so great as to allow the grain to go into suspension. That is, saltation is
analogous to a ballistic trajectory.
Many have argued that sediment transport may be divided into two types: bedload
and suspended load. These two types of sediment transport are distinguished by the
dominant force mechanisms. Grain-to-grain forces are the dominant mechanisms in
bedload transport and fluid forces dominate suspended load transport; therefore, rolling
would be considered bedload transport and obviously suspension would be considered
suspended load transport. Saltation, on the other hand, does not fit neatly into either of
these categories. A grain is driven into saltation by another grain impacting it or by the
grain rolling and striking a neighbor grain in front of it, forcing the striking grain to move
into the fluid column. This scenario prescribes a grain-to-grain interaction; hence,
saltation may be considered bedload transport. Once the grain leaves the bed, however,
fluid drag, inertia, and other fluid forces act on the grain. Therefore, saltation could be
considered as suspended load transport. It is unclear how to categorize saltation in the
context of a bedload/suspended load dichotomy.
Perhaps more important than what category saltation fits, is the fact that it plays a
significant role in sediment transport. According to Bridge and Dominic (1984), saltation
is the dominant mode of bedload sediment transport. Rolling and sliding of grains does
not occur as often. This mode usually appears between individual saltation events (Bridge
and Dominic, 1984).
Saltation has been seen while conducting physical models. A flume, with an
adjustable bed slope, is used for these experiments. The bed roughness is set by
permanently placing particles on the base of the flume. Then bedload grains are allowed
to move freely in a shear flow. If the shear stress is great enough, saltation occurs.
Photography and video have been used to record grain flow characteristics. The initial
velocity and angle of a saltating grain is estimated through the use of photography (Abbott
and Francis, 1977) and video (Nino and Garcia, 1992). These quantities are only
estimates because it is difficult to see a grain leave the bed. The photographs and video
track a grain travel through its trajectory. Thus, a set of saltation trajectories are obtained
along with some initial condition estimates. Attempts are then made to match the paths
with numerical models.
The usual technique for developing a predictive model is to define all of the forces
acting on a grain, using numerical simulation, which includes the 'law of the wall' fluid
velocity profile and assuming ejection of grains from the bed with an initial velocity and
angle. The forces that act on a grain include gravity, drag, inertia, lift, Magnus lift, and
Basset history force. Some combination of these forces is used by most researchers
depending on the particle size and flow intensity. An equation of motion is developed and
solved numerically to yield a particle trajectory. From many trajectories, a concentration
profile may be obtained. Finally, a bedload transport model is developed from the
concentration and velocity profiles (Wiberg and Smith, 1989).
The focus of this study is to evaluate the existing models on their ability to predict the
trajectories of saltating grains. The next step then is to develop a new saltation model.
This model incorporates the best aspects of the earlier models and adds some new
features. The single grain model yields a particle trajectory as well as presenting figures
that demonstrate the relative importance of forces acting on the particle. In addition, a
random trajectory model is made as an extension of the single grain model where
numerous saltation paths are found. The grains are ejected against and with the shear
flow. This is the first time a particle saltation model will include grains ejecting into the
flow.
CHAPTER 2
PREVIOUS WORK
Mechanisms of Grain Saltation
Incipient Motion
Incipient motion occurs when the shear stress is great enough for the sediment to be
carried as bedload. That is, the grain begins to roll or slide along the bed surface. If the
shear stress increases further, the grain will gain more momentum and eject from the bed.
The grain will then travel in a ballistic manner and eventually strike the bed. It will either
rebound from the bed again or eject another grain through a momentum transfer. This
process continues until the shear stress lessens.
A nondimensional Shields' parameter is used to determine incipient motion. It is
defined as
0 = ( b equation 2.1
p (s 1)g d
where Tb is the critical shear stress, p is the density of the fluid, s is the specific gravity of
the grain, g is gravity, and d is the grain diameter. Shields (1936) plotted Oc values versus
grain Reynolds number (Re*),
u,d
Re, = equation 2.2
where u is the shear velocity of the fluid and v is the cinematic viscosity of the fluid.
where u. is the shear velocity of the fluid and v is the kinematic viscosity of the fluid.
For small grain Reynolds numbers, the Shields parameter varies linearly with the
Reynolds number. As the Reynolds number increases, the critical value goes through
some transitions with a minimum of approximately 0.03 at Re. equal to 10. For large
values of Re*, the critical value increases to 0.058 (Fredsoe and Deigaard, 1992).
Figure 2.1 shows the Shields' curve plotted against Archimedes number (Ar),
Ar = "2 d equation 2.3
where ps is the density of the particle.
100
_10
.-\
10-2
100 101 102 103
Archimedes Number
Figure 2.1 Critical Shields' Parameter As A Function of Archimedes Number (Sleath,
1984)
The Archimedes number can be a more convenient formulation to use if the shear
velocity is the desired quantity. That is, the original Shields' curve plotted versus the
grain Reynolds number requires iteration to obtain the shear velocity since it appears on
both axes. With Archimedes number replacing the grain Reynolds number, the shear
velocity is obtained directly.
Particle Saltation
When critical shear stress has been surpassed, sand grains will continue to build
momentum. As previously stated, the grain will roll along the bed and strike another
grain. This may cause the rolling grain to leave the bed and enter the shear flow. Once in
the fluid column, the grain will follow a ballistic trajectory unless the flow shear stress is
great enough to induce suspension. This situation can be attributed to drag forces. The
flow can also lead to lift forces due to the shearing across the grain surface that cause the
grain to leave the bed as well. The hydrodynamic forces, shear lift and Magnus lift can
account for this phenomenon (Nino, Garcia, and Ayala, 1992).
Einstein (1950) examined saltation in deriving his stochastic bedload formula. He
assumed that the particle's saltation length was proportional to grain size. It was later
found that the saltation length is inversely proportional to grain size.
Bagnold (1973) investigated saltation in water. He defined saltation as particles
transported by consecutive 'hops' over a bed driven by a fluid. He stated that the particles
'hop' by contacting other grains and the bed. The initial momentum comes from a grain
contacting the bed and rebounding higher into the flow.
Fernandez Luque and van Beek (1976) performed experiments in a flume with
different bed slopes. They were able to measure the mean critical bed shear stress, rate of
bedload transport, average particle velocity, and the average length of individual
saltations. This was accomplished through the use of high speed photography. Two
different sand types, gravel and magnetite were studied. When they compared their model
results, they concluded that a lift force was needed to explain the observed saltation
characteristics.
Abbott and Francis (1977) also conducted a series of experiments of particle saltation
in water. They too were able to obtain mean trajectory lengths and heights in addition to
the mean particle velocity with high speed photography. The particles that were studied
ranged in size from five to ten millimeters. Natural and artificial particles were both
utilized. The conclusion drawn again was that there were grounds for the existence of a
force that opposes gravity when a particle is in a vertical gradient of horizontal velocity.
Chepil (1958) examined the ratio of drag and lift forces on a fixed particle under the
influence of wind. At different wind intensities, the drag and lift forces on the particle
were calculated to determine a ratio. In addition to the lift force Chepil (1958) calculated,
a spinning lift force was also investigated. White and Schulz (1977), while inspecting
saltation in air, found that particles traveled higher and farther than if they were in a
vacuum. First they tried to explain the additional height and length with drag only. This
force proved insufficient. They concluded that the extra height and length was a result of
a Magnus or spinning lift force based on the work done previously by Rubinow and Keller
(1961). The latter derived the Magnus lift force and moment on a spinning sphere in a still
fluid.
Saltation may also occur when many layers of grains are in motion on the bed, in what
is commonly called sheet flow. Sheet flow conditions occur when the shear stress is
substantially greater than that required for grains to go into suspension. When sheet flow
is occurring, all modes of transport, rolling, saltation, and suspension may be taking place
as well. Sheet flow can cause grains to saltate through an event analogous to popping
popcorn being ejected from an air popper. The popped kernels are 'ejected' from the air
popper at different angles and speeds. As they strike the popcorn already in the receiving
bowl, they in turn eject popcorn at random angles and speeds. A similar effect can be
applied to sheet flow with sand grains striking each other. That is, the grains would be
ejected from the bed randomly. This can lead to grains saltating.
There is some debate as to what angles the saltating grains eject from the bed. Owen
(1964) examined grain saltation in air. When looking at the initial conditions for a grain
ejecting from the bed, he assumed that particles had nearly a vertical initial trajectory.
Since the height of a particle is related to the square of the vertical velocity, only those
particles that leave the bed nearly vertically with a certain initial speed will be able to
saltate (Owen, 1964). Others like White and Schulz (1977) believed that grains eject at
angles ranging from 30 to 69.5 degrees. This was also observed in air.
Other researchers have tried to derive a set of equations to describe the motion of
a particle in saltation from bed ejection to bed impact based on the fluid forces. Those
who developed models include Murphy and Hooshiari (1982), van Rijn (1984), Wiberg
and Smith (1985), Nino and Garcia (1992), and Lee and Hsu (1994). All of the models
explored by these researchers are from the Langrangian perspective. These models are
described in the section titled "Predictive Models."
Force Analysis
The search for forces has been motivated by the fact that computed trajectories based
on measured initial velocities and angle did not agree well with measured trajectories.
Models based on drag and gravity proved insufficient, so others were sought. Some of
these other forces are agreed to be significant while others can be ignored in certain
situations. Below is list of forces that researchers have used in saltation modeling.
Gravity Force
A force that all researchers can agree is the influence of gravity on the path of a
particle. Gravity always opposes the vertical ejection of a grain from the bed and usually
helps the horizontal ejection of a grain from the bed. The effect of gravity is usually
written as a submerged weight,
FG = (P p)gV equation 2.4
where V is the volume of the grain (van Rijn, 1984). Obviously, this force increases with
increasing grain size because of its dependence on the volume of a grain.
Added Mass Force
This force arises from the relative accelerations of the grain and the fluid. A
submerged body induces accelerations on a fluid if the body is moving with an acceleration
relative to the surrounding fluid. It can be thought of as a body having an 'added mass' of
fluid attached to its own mass when a submerged body induces an acceleration to some of
the surrounding fluid (Patel, 1989).
In general, the added mass force may be written
V fdu (duv dua\
FA du = p -du equation 2.5
Sdt dt dt)
where uf is the fluid velocity, u is the grain velocity, and CM is the added mass coefficient.
The added mass coefficient is defined as a ratio of the additional mass of fluid that is
accelerated with the particle to the mass of the displaced fluid by the particle (Dean and
Dalrymple, 1984). For a sphere, CM equals 0.5.
For a steady, shear flow, the added mass force may be reduced to
FA = dCV equation 2.6
dt
where the force acts opposite the grain's acceleration. The inertial effects are excluded.
Drag Force
A drag force is a fluid force exerted on a body. More specifically, it is a net force in
the direction of the fluid relative to the sphere due to pressure and shear forces on the
body. The drag force on a particle may be written as
F, = CDA pV,2 equation 2.7
2
where CD is the coefficient of drag, A is the cross sectional area of the grain normal to the
force, and Vr is the relative velocity (Fredsoe and Deigaard, 1992). The drag coefficient is
a strong function of Reynolds number and shape; thus, it is generally not constant. There
have been many empirical formulas such as those found by Graf (1984) and Morsi and
Alexander (1972). As to which formula is the best depends on the range of Reynolds
numbers that are encountered.
Shear Lift Force
This force is a direct result of a shear flow. The particle or sand grain develops a
pressure gradient across it and as a result causes a lift force. This phenomenon can be
attributed to the Bernoulli effect where the lift force acts up if a positive gradient is
present and acts down if a negative gradient is present. That is, the lift force acts up if
Ufrop is greater than Ufot and down if UfBot is greater than ufTop.
fr (top) ur (top)
uf (bot)
ur (bot)
Figure 2.2 Shear Lift Results From Velocity Gradient Across Grain
The shear lift force may be written as
F, (shear) = CLA p(u2To u ) equation 2.8
where CL is the lift coefficient and A is the cross-sectional area of the grain normal to the
force. The UATop and UABot are the relative velocities evaluated at the top and bottom of the
grain, respectively (Wiberg and Smith, 1985).
The lift coefficient has been related to the drag coefficient by Chepil (1958). He
conducted experiments involving evenly-spaced hemispheres and allowed wind to flow
over them. He then proceeded to calculate drag, lift, and the ratio of lift to drag for
different wind speeds. It was determined that ratio was approximately equal to 0.85;
therefore, CL equals 0.85 times CD (Chepil, 1958).
Based on the above definitions, in a wall-bounded shear flow, the lift force will act
to lift the grain off of the bed when the local fluid velocity is greater than the grain
velocity. If the grain velocity is greater, the lift force will cause to the grain to accelerate
toward the bed as shown in figure 2.2.
Another definition of the shear force was developed by Saffman (1965). This
expression is for a sphere moving in a viscous flow. It is
S du 2
F,(shear)= CLpt2d2V 'z equation 2.9
Sdz )
where CL is a lift coefficient but different than Chepil's (1958). Saffman determines his lift
coefficient as a function of the grain Reynolds number, Re*; CL equals 1.6 when Re. is less
than 5 and equals 20 when Re* is greater than 70. CL varies linearly between those Re*
ranges.
The two shear lift forces are very close even though the lift coefficients are different
(Lee and Hsu, 1994). The Saffman expression does not have a negative lift associated
with it since the velocity gradient is always positive. That is, this lift force definition
always has positive lift in the plane normal to bed even if the grain is at a higher velocity
than the surrounding fluid (Lee and Hsu, 1994).
Basset History Force
A. B. Basset (1888) first acknowledged that a particle's history had a role in the
present particle path; hence, the force bears his name. Mei (1995), as told to the author,
described the force as
... derived from the diffusion of vorticity generated at the surface of the
particle at a rate proportional to the particle's relative acceleration. Since
the diffusion rate is finite, this means that the force is dependent on the
history of the particle motion. (Mei, Personal Communication, 1995)
Consider a sphere in a steady fluid and then instantaneously increasing the flow to
some higher value. It takes some 'time' for the boundary layer on the sphere to adjust to
the new flow intensity. This 'time' is accounted for through the Basset history force.
The Basset force is defined as
Sd 2 du duf
F f- dt dt dt equation 2.10
0 S
where pL is the dynamic viscosity, Ts is the total saltation time, and T is a dummy variable
(Mei, 1994). This force has both a vertical and horizontal component. For steady, shear
flow, the fluid acceleration is zero. So the term can be simplified to be a function of the
grain acceleration.
Magnus Lift Force
In addition to the shear lift force, the Magnus force results from the velocity gradient
across the grain, too. This force, named for Heinrich Magnus who first discovered the
phenomenon in 1853, is a pressure force due to the circulation around a spinning sphere
(Murphy and Hooshiari, 1982). As a result of viscous effects, a rotation is supplied to the
particle. This force accounts for the different types of pitches in baseball such as the
curveball and the slider (Munson, Young, and Okiishi, 1990). The shear flow induces a
rotation which causes the grain to saltate higher and further than without the inclusion of
this term. If the velocity gradient is positive, the grain will rotate clockwise (cw); hence,
an upward lift. If the velocity gradient is negative, the grain will rotate counterclockwise
(ccw); hence, a downward lift. See figure 2.3 below.
ur (top) ur (top)
grain ) FL grain FL
uf (bot) uf (bot)
Figure 2.3 Magnus Force Results From Velocity Gradient Across Grain And Viscous
Effects
For a moving sphere in a shear flow, the force is expressed as
F, (Magnus)= 8 d3pV Q-- equation 2.11
8 2 dz
with Q as the angular velocity of the grain with units in rad/s (White and Schulz, 1977).
This force was developed from the work of Rubinow and Keller (1961). Rubinow and
Keller (1961) looked at rotating a moving sphere in a still viscous fluid with low Reynolds
numbers only. Based on their analysis, the Magnus force was independent of viscosity.
From this they derived the form of the Magnus force and the moment acting on the
sphere. The force has been described previously. The moment has the form
Moment = I d= -_ 3 1 dUf equation 2.12
dt 2 dz
where I is the grain's moment of inertia. This moment equation is solved simultaneously
with the equations of motion to constantly adjust the grain rotation to that induced by only
the fluid. Note that the dynamic viscosity appears in equation 2.12. It acts to dampen the
effects of the initial particle rotation to that of the fluid. It must be said that the Magnus
force and moment equations are only valid for small Reynolds numbers. They are used in
this shear flow analysis only for the purposes of determining the effect the Magnus force
has on improving the correlation between experiments and theory (White and Schulz,
1977).
Applicability of Forces For Different Reynolds Number Ranges
The drag, added mass, shear lift, Magnus lift, and Basset history forces have been
formulated as previously published in the literature. The drag force, as defined here, is
valid for all Reynolds number flows. The form of the shear lift force was originally verified
for turbulent flow. Also, theoretically a solid sphere in an inviscid fluid has an added mass
coefficient value of 0.5; so, the added mass force defined herein is applicable to large
Reynolds number flows, too. From Rubinow and Keller (1961), the Magnus force was
derived for low Reynolds number flows only. The Basset history force is also not well
known for low Reynolds number flow (Nino and Garcia, 1992). For analysis purposes,
the forces are extended to the entire range of Reynolds numbers. This follows the
previous research procedures as done by Wiberg and Smith (1985) and Nino and Garcia
(1992) which have yielded satisfactory results.
Predictive Models
Murphy and Hooshiari (1982)
As mentioned earlier, Murphy and Hooshiari (1982) developed a model of particle
saltation. They obtained their own data set based on marble saltation with a diameter of
1.57 cm and a density of 2.49 g/cm3. The basis of their model was to neglect those forces
proportional to the relative velocity so that the result is a particle moving in a vacuum
under the influence of an equivalent gravity, ge,
2g(p, p)
ge = ( p) equation 2.13
2p, +p
The simple dynamic model, as it is called, yielded saltation characteristics such as height,
length, and duration of the trajectory. They are for height
2
H= equation 2.14
2g,
for length
L 2uo equation 2.15
ge
and for duration
2w
T 2wo equation 2.16
ge
where Uo is the initial horizontal velocity of the sphere and wo is the initial vertical velocity
of the sphere. They also included an additional length to the length of a trajectory by
including the effect of added mass due to fluid shear. This length is given by
equation 2.17
where
b=240p(p, p)
(2p, +p)2
equation 2.18
equation 2.19
with u* as the fluid shear velocity. The extra length to the trajectory comes from the
logarithmic fluid velocity profile given by
Uf = 2.5u.ln ..j
equation 2.20
where z is the vertical position of the center of the sphere and zo is the vertical position
where the fluid velocity is zero. The initial conditions, initial vertical and horizontal
velocity, were both based on the first observable, short, vertical displacement of the
marble from the bed.
The breakdown in the simple dynamic model is the fact the drag and lift both increase
as the initial vertical velocity of the sphere increases. According to this model, the
saltation height increases as the square of the initial vertical velocity. This means that the
height is not accurately predicted. Murphy and Hooshiari (1982) acknowledge this aspect.
They, however, maintain that the length of a trajectory is predicted accurately. They
needed to add the drag force back into the particle equation to more accurately match
wWubb' 1I l,
ALI -tah-
C1 gd 2i
C +W2
their obtained data set. After doing this, they believe the lengths of the trajectories are
matched fairly well.
2.5 I
2
E
E1.5
0o
I
1
0.5
Length (cm)
Figure 2.4 Comparison of Typical Trajectory and Murphy, et al. (1982) Model
From figure 2.4, the model appears to describe the path of a marble fairly well;
however, it may be too much of a simplification to exclude the lift forces entirely since
those forces may determine the height of a particle trajectory.
van Rijn (1984)
A model developed by van Rijn (1984) included added mass, gravity, drag, and a
shear lift force. The lift force used was the one developed by Saffman (1965) as described
previously in the section "Force Analysis." Recall that the expression was valid only for a
sphere moving in a viscous flow. The shear lift force has two components: one in the
plane normal to the bed and one in the plane horizontal to the bed. The Magnus force was
0 O
oo
00
0-
0 Experiment
m Model
neglected because Saffman (1965) showed that for a viscous flow, the Magnus force was
an order of magnitude smaller than the shear lift force (van Rijn, 1984).
The fluid velocity was given as the standard logarithmic profile with the velocity equal
to zero at the bed height, zo
z0 = 0.11 + 0.03k, equation 2.21
where ks is the equivalent bed roughness. The equivalent bed roughness is adjusted to
allow the model to match observed trajectories (van Rijn, 1984).
The initial conditions were obtained through others' observations. The initial
velocities of the grain used were gathered by Abbott and Francis (1977). They found that
the average initial velocity was approximately 2u*. White and Schulz (1977), looking at
saltation in air, found that the velocity varied in the range u* to 2u*. The take-off angles
from the bed were observed to vary from 300 to 700 (White and Schulz, 1977).
The theoretical bed is located at 0.25 times the grain diameter, d, below the layer of
grains that constitute the plane bed. The saltating grain is located 0.6 times the grain
diameter above the theoretical bed. Now the equation of motion can be solved
numerically.
A typical model trajectory is shown in figure 2.5 with an observed trajectory provided
by Fernandez Luque and van Beek (1976). The particle diameter, d, is 1.8 mm, s is 2.65,
u* is 0.04 m/s, the kinematic viscosity is 0.01 cm2/s, and the initial velocity is taken as 2u*.
The adjustable parameters in the model are the lift coefficient, bed roughness, and initial
velocity and angle. For the figure, the bed roughness, ks, is a grain diameter and the lift
coefficient, CL, is equal to 20. The initial angle is 450. The Fernandez Luque and van Beek
trajectory is probably wavy because of the upward fluid forces due to turbulence (van
Rijn, 1984).
I
E2.5
a,
I
c= 2
0
U)
C,
1.5
0
z
0 5 10 15 20 25 30
Nondimensional Length, L/d
Figure 2.5 Observed Trajectory and Model Comparison (van Rijn, 1984)
Van Rijn (1984) was able to conclude from his model that small particles are
dominated by drag and therefore have long, flat trajectories. Large particles are
dominated by the shear lift force, particularly near the bed, and therefore have short, high
trajectories.
Wiberg and Smith (1985)
Wiberg and Smith (1985) derived an equation of motion that included gravity, added
mass, drag, shear lift force, Magnus lift force, and a form lift force. All but form lift have
been described previously under "Force Analysis." There is only a vertical component of
shear lift according to Wiberg and Smith (1985). The shear lift force coefficient that is
used is held constant at 0.2 as a representative value of the lift to drag ratio. In addition,
the Magnus force is defined here as
F (Magnus)= 27t pvQ equation 2.22
which is the theoretical form of an infinite, rotating cylinder in a two-dimensional flow
(Batchelor, 1967). They use this as an approximation of the force on a sphere. Wiberg
and Smith (1985) then proceeded to state that the Magnus force will cause the grain to
saltate too long and high when compared to observed results. Therefore a negative lift at
the top of the grain trajectory is needed to keep the hop length from increasing to such a
high value (Wiberg and Smith, 1985).
In response to the concerns expressed above, they defined a form lift force which acts
on non-spherical bodies. This lift force is the same as the one that acts on airplane wings.
Based on the general expression for form lift on an elliptic cylinder, the force is defined as
FL(form)= 4cpvc, sin(N1y + N2) equation 2.23
where c2 is proportional to the length of the centerline of the airfoil, YJ is the flow's angle
of attack relative to the leading edge, and T2 is a function of the curvature of the airfoil's
centerline. Here they take Y2 to be equal to zero and adjust the circulation so that the
trailing stagnation point is at the rear of the cylinder's long axis. As is shown in figure 2.6,
this force can be positive or negative based on the orientation of the particle (Batchelor,
1967).
FL>O FL
Figure 2.6 Sign of Form Lift for Different T' Angles of a Body Relative to Flow
Of course to use the form lift force, an elliptic grain was assumed. Wiberg and Smith
(1985) oriented the grain in a manner such that the form lift is a maximum at the beginning
of the trajectory and allows the grain to rotate so that the form lift is negative at the top of
the trajectory. For the downward portion of the trajectory, the grain orientation does not
change significantly (Wiberg and Smith, 1985).
The Magnus and lift forces are calculated together with one value of circulation that
accounts for the relative magnitude and signs of the two individual circulations. Once the
combined forces per unit length are calculated, the value is multiplied by the diameter of
the grain and by a coefficient. This coefficient is determined by comparing the model
result to the height of known trajectories, namely, Francis (1973). Figure 2.7 displays the
model result versus a Francis (1973) photographed trajectory with d = 0.7 cm, s = 1.3,
and u. = 6.6 cm/s. The velocity profile is the standard "law of the wall" profile and the
grain is initially positioned at one-half a grain diameter above the theoretical bed.
:E
rn
.0
a)
a:
0 2 4 6 8 10 12
Length (cm)
Figure 2.7 Comparison of Predicted and Observed Grain Saltation Trajectory
The initial vertical and horizontal grain velocities are determined by decomposing the
forces on a grain into radial and circumferal components in polar coordinates. An
expression for the centrifugal force is obtained as a result of integrating the circumferal
force with respect to the angle the pivoting axis of the grain makes with the vertical. The
centrifugal force expression and the one for the radial force are balanced when the
pivoting angle is some critical value. The initial velocities are obtained from the resulting
equation. When this is done, the equation of motion is solved with a fourth order Runge
Kutta method (Wiberg and Smith, 1985).
From their model, Wiberg and Smith (1985) concluded that the saltation
characteristics, height and length, both increase with increasing particle density. This was
not supported by Abbott and Francis (1977). Without more data, Wiberg and Smith
(1985) state that it is very difficult to determine if the problem with the model is with the
force combination or the with the way the grain is oriented for the form lift force (Wiberg
and Smith, 1985).
Nino and Garcia (1992)
A gravel saltation model was developed by Nino and Garcia (1992) based on the data
they accumulated from Nino, Garcia, Ayala (1992). They include the forces of added
mass, gravity, drag, shear lift, Magnus lift, and Basset history. These forces appear as
previously described in section "Force Analysis". The shear lift force form is the same as
for Wiberg and Smith (1985) with CL equal to 0.2. The Magnus lift force is of the same
form as used by White and Schulz (1977) with the fluid velocity gradient as a part of the
force. The Basset history force is included for the first time of the models examined. This
force is not integrated by "brute force" means in the equation of motion, but it is
approximated using the analysis of Brush et al. (1964).
The fluid velocity profile is the turbulent velocity profile over a rough boundary given
as
uf=2.5 ln(30z) equation 2.24
where the roughness height is taken as equal to one grain diameter. The initial velocity of
the gravel was determined by examining the previously obtained video of gravel saltation
provided by Nino et al., (1992). The particle is initially positioned one-half grain diameter
above the theoretical bed. The angular velocity of the gravel was determined by choosing
the best value that allowed the model to match an observed gravel trajectory. A typical
trajectory is plotted in figure 2.8.
1.2
.2)1.1
c")
I 1
0 1 2 3 4 5 6 7
Nondimensional length, L/d
Figure 2.8 Comparison Between Predicted and Observed Trajectory (Nino et al., 1992)
Nino et al. (1992) were able to conclude the following: 1)The shear lift force does
indeed add to the agreement between theoretical and observed trajectories, 2)The Basset
history force had little effect on the predicted trajectory for gravel and 3)Particle rotation
does bring the predicted trajectory closer to the observed, however, only for some cases.
(Nino et al., 1992). The problem might be a result of the constant rotational value
assigned to the gravel instead of having the particle rotation adjust to the flow.
Lee and Hsu (1994)
Lee and Hsu (1994) closely followed the model developed by van Rijn (1984). The
main exception is that they included the Magnus force and compared the two shear lift
forces explored by Wiberg and Smith (1985) and Saffman (1965). The initial conditions
were very similar except that the initial position of the grain is at 0.5 times grain diameter
above the theoretical bed.
The two shear lift forces were calculated for a run and the values were comparable
even though the lift coefficients are very different. The Magnus force was explored by
trying different values to best match Fernandez Luque and van Beek (1976). Lee and Hsu
(1994) came to the conclusion that the trajectory has three distinct regions: rising, central,
and falling branches. For each region, there is an angular velocity associated with it. The
greatest rotation is in the rising branch and decreases through the other two branches.
- 2u* at 45 deg
2.5 ,'
-- 5u* at 45 deg
/ \
S2h S
o 2
E 1.5
0
z
0 2 4 6 8 10 12 14 16 18 20
Nondimensional Length, L/d
Figure 2.9 Trajectories With Same Angle and Different Velocity (Lee and Hsu, 1994)
When excluding the Magnus force altogether, an interesting result is obtained. As a
result of the shear lift force having two components (van Rijn, 1984), two particles fired at
the same angle with different velocities will exhibit a strange behavior. The particle with
the higher initial velocity will not go as high nor as far as the particle with the smaller
velocity. This is shown in figure 2.9. One particle has an initial velocity of two times the
shear velocity while the other has five times the shear velocity. Both were fired at the
26
same angle. The grain diameter is 1.36 mm with a shear velocity of 3.6 cm/s, CL equal to
20, and a ks of 2d. From this model it appears that maybe the shear lift force should act
only in the plane normal to the bed as done in Wiberg and Smith (1985) to achieve a more
physically realistic result.
CHAPTER 3
SALTATION MODEL FORMULATION
Single Grain Model
Equation of Motion
The objective of the single grain model is to yield a particle saltation trajectory. The
model developed by the writer contains many of the same forces as were used in the
previously discussed models. The forces that are considered important are added mass
(Fa), gravity (Fg), drag (Fd), shear lift (Flshear), Basset history (Fb), and Magnus or spin lift
force (F1Magnus). Recall that Saffman (1965) and Wiberg and Smith (1985) derived two
forms of the shear lift force. Both can be used in this model by simply exchanging one for
the other. This allows for an easy comparison of the two forms. The Basset history force
is used but in a different form than what is presented in Nino and Garcia (1992). They
borrowed from the simplifications of Brush et al. (1964) whereas in the writer's model the
integro-differential equation is handled directly. The form of the Magnus force is the same
as the form developed by White and Schulz (1976). Their form includes the constant
updating of the particle's angular velocity to that of the rotation induced by the fluid. All
the forces listed have a normal and a parallel component to the bed except for the shear lift
force. Since there exists a shear only in the plane parallel to the bed, the shear force has
only the component normal to the flow.
The forces are more readily seen in figure 3.1. The equation of motion for a saltating
grain may be divided into longitudinal and vertical components. The forces acting on a
saltating grain are described by the following equations:
tpV = -CDA
d 2 d
PSDAPu-(u)+w2(-u,)
-p VC du d
P du 6dt dt
dt 0r o -,-
S 2 p+-1 du, w
+-d22
8 2 dz (2uuy +W
+p(s )gV sin
and
dw 1 /\2
pVCdt = -CA p (u-u) +w(w)
8 2 dz (uf +w2
+2-CL~~TO ~0~- ps -1)gc u j3 ,
+ ACL(Urop -Uo )- p(s 1)gVcos P
where
where
equation 3.1
equation 3.2
2u = (U _-fop +W2
UATop ".- W2
equation 3.3
u2 Bot e-)2 +W2
UABot U UfBot-
equation 3.4
The variables, urrop and UfBot, are the fluid velocities evaluated at the top and bottom of the
grain, respectively. The quantities, u and w, are the particle's speed parallel and normal to
the bed, respectively. All these assume no w-component of fluid velocity.
Equation 3.1 is the bed-parallel component and equation 3.2 is the bed-normal
component of the equation of motion. The terms for equation 3.1 are, in order, mass
times acceleration, drag, added mass, Basset history, Magnus lift, and gravity. The terms
for equation 3.2 are mass times acceleration, drag, added mass, Basset history, Magnus
lift, shear lift, and gravity.
Fl(shear)
e Fl(Magnus)
Vr
/ Fd& '
wO ; Fb
uO 0.5d
Figure 3.1 Force Definition Sketch of Grain Saltation
Bed
Grains ---
Horizontal Plane
Figure 3.2 Definition of Bed Slope
The shear lift force currently in the model is the one derived by Wiberg and Smith
(1985). To compare the trajectories that result from the model with the shear lift of
Wiberg and Smith (1985) and those that result with the Saffman (1965) shear lift
expression, equation 3.2 needs to be modified. The shear lift term is replaced by
F = -C ,p0.5d2 U2+W2 C duf u-u,
r dz ) ( uf)2 +w2
equation 3.5
where this CL is related to the grain Reynolds number, Re* (Saffman, 1965). When using
the Wiberg and Smith (1985) shear lift force, the lift coefficient, CL, is taken to be 85% of
the drag coefficient, CD.
The added mass coefficient, CM, is taken to be 0.5 for a sphere. The coefficient of
drag, CD is calculated based on the empirical equations of Morsi and Alexander (1972).
Their empirical curves may be divided into the following regimes:
24
C =
Re
22.73 0.0903
C = + +3.69
D Re Re2
29.1667 3.8889
C +1.222
Re Re2
46.5 116.67
C = _+0.6167
D Re Re2
98.33 2778
C Re Re2 + 0.3644
S Re Re2
148.62 4.75e04
C Re 0.357
Re Re2
for Re<0.1
for 0.1
for 1
for 10
for 100
for 1000
equation 3.6
equation 3.7
equation 3.8
equation 3.9
equation 3.10
equation 3.11
-490.546 57.87e04
490.546 57.87e4 0.46 for 5000
S Re Re2
-1662.5 5.4167e06
CD 16. + +e6 0.5191 for 10000
Re Re2
The piecewise coefficient of drag formulation is chosen because it is very complete
and does not require much effort to use. It can, however, be replaced with a universal
equation valid for the entire range of Reynolds numbers.
The fluid velocity profile used in the model is the standard logarithmic "law of the
wall" profile with the velocity equal to zero at
zo = 0.11 + 0.03k, equation 3.14
U,
This profile is the same as van Rijn (1984) used. Note that the fluid velocity only has a
horizontal component. This fluid velocity can be easily modified in the model should it be
desired.
Boundary Conditions
The bed level is assumed to be 0.5 times grain diameter below the top of the grains as
is shown in figure 3.3. The initial position of the saltating grain is at one-half a grain
diameter above the theoretical bed. Although the grains aligned in the bed may saltate at
sufficiently high shear stresses, the model herein only yields the path of a grain at the given
position above the theoretical bed.
In order to solve the equations of motion given previously, the grain's initial vertical
and horizontal velocity components are needed. The values of these components come
from White and Schulz (1977). They found that the velocities varied from u* to 2u*. They
also found that lift-off angles varied from 30 to 69.5 degrees. In addition, when the
z
S .. Saltating
/ IGrain
0.2d zo
Figure 3.3 The Saltating Grain Initial Position
Magnus force is considered, an initial particle angular velocity is needed. This value is
adjusted to yield a best match to a known trajectory.
Method of Solution
The equations of motion along with the moment equation are defined as first order
ordinary differential equations. A fourth order Runge Kutta approach is used to yield a
vertical and horizontal grain velocity and an angular velocity. The vertical and horizontal
velocity components are numerically integrated with a simple Simpson's Rule approach to
obtain a particle trajectory.
The model includes some adjustable parameters. The bed roughness can be taken to
be any multiple of the grain diameter. The aforementioned initial velocities and angular
velocity also may be adjusted. Also, the initial grain ejection angle from the bed can range
from 0 to 180 degrees with the former meaning with the shear flow and the latter meaning
directly into the shear flow. This is the first saltation model to allow for the inclusion of
ejecting grains from the bed into the fluid flow.
Sensitivity Analysis
This section looks at the relative effect force combinations and initial conditions
have in determining the trajectory of a saltating grain. The parameters that are held
constant unless otherwise specified through the force sensitivity test are grain diameter,
specific gravity, initial particle velocity, initial angle, bed roughness, and the coefficient of
lift. They are 0.18 cm, 2.65, 2u*, 45 degrees, 2d, and 0.2, respectively. These are the
values that were either used or observed by Fernandez Luque and van Beek (1976). In
addition, it needs to be stated that the bed slope was taken to be zero in this case. The
critical velocity for incipient motion is found by calculating Archimedes number and using
a look-up table to find Shields' parameter. From this, the critical velocity is determined
with the relation,
Ucr -= equation 3.15
where Tcr is the critical shear stress.
The most basic form of the model contains the drag, added mass, and gravity forces.
The shear lift force developed by Wiberg and Smith (1985) and the history force from
Basset (1888) are successively added to the model. Figure 3.4 shows the results. For this
case, u. = 4 cm/s with uo = wo = 2u*.
The shear lift force increases the length of the trajectory by about one grain diameter
while the Basset history force appears to have a slightly lesser effect on the trajectory.
The results with the shear lift force included are consistent with our intuition. Only a small
velocity gradient develops across the grain because of its size. As a result, the shear lift
has a relatively small effect. The Basset history force is not as easy to predict. The height
is increased by less than one-tenth of a grain diameter while the length increased by
approximately one-half grain diameter when it is included. The increase in length is due to
the fact that the grain travels higher into the fluid column and thus attains a greater
velocity. This increase in length may occur because the vortices have a large impact on
small particles.
0.5 I I X I \ I
0 0.5 1 1.5 2 2.5 3 3.5 4
Nondimensional Length, L/d
Figure 3.4 Particle Trajectories With A Particle Size Of 0.18 Centimeters
The next test was to run the model with a larger particle size. A diameter of 3.1 cm
was chosen, roughly in the gravel regime. The other parameters maintained their same
values with the exception of u* = 22.84 cm/s and P = 0.07 to match the experimental
conditions of Nino et al. (1992). Figure 3.5 shows that the shear lift force had more of an
effect than it did with the smaller grain. This result is expected because the velocity
gradient is large on a larger particle. The Basset force, however, did not posses the same
significance as it did in the other case. This force only slightly increased the particle
trajectory. The same result was obtained in Nino and Garcia (1992) with regard to the
Basset history force. Both forces do demonstrate that they act to increase saltation height
and length as a result of being included in the model.
1.4
No Shear Lift, No Basset
.. .. ..
0 2,. "4.6 \
8- I ."
/ -
08 \
Figure 3.5 Particle Trajectories With A Particle Size Of 3.1 Centimeters
The next force to consider adding to the formulation is the Magnus lift force. To
review, this force results from two effects. First, a grain may be transferred an angular
velocity from the shear flow. Since the fluid velocity is greater at the top of a grain than at
the bottom of a grain, it may acquire a rotation. Second, a grain on the bed may obtain an
angular velocity from another grain striking it. The magnitude of the rotation is dependent
on the speed and the placement of the blow that the grain striking the bed delivers.
Figure 3.6 displays the trajectories that result by varying the initial angular velocity.
The diameter of the particle is 0.18 cm with ks=2d, Cl=0.2, u*=4 cm/s, uo=2u*, and
wo=2u*. It is obvious and expected that the trajectories should both increase in length and
height with increasing initial angular velocity. The effect of the Magnus force on a
particle's path is quite significant. This is a discouraging result since little is known about
what angular velocities are appropriate for saltation.
/
I I
S
0.8 -
0.6
0 1 2 3 4 5 6 7 8 9 10
Nondimensional Length, L/dinitial angles are 30, 45, and 60
Figure 3.6 Particle Trajectories Sensitivity to Particle Rotation
As stated previously, the initial conditions of the model greatly effect the resulting
salvation trajectory. The initial angle and velocity of the grain is the next area of interest.
Again the grain size used is 0.18 cm. Figure 3.7 displays the results of varying the take-
off angles and maintaining a constant initial velocity. The initial angles are 30, 45, and 60
degrees. The initial longitudinal velocity and initial vertical velocity are both 2u*. Figure
3.8 shows what happens if the velocities are varied and the angles are held constant at 45
degrees. As the lift-off angles increased, the grains saltated farther. Also, the higher the
initial velocity, the longer the trajectory. These results are obviously from the fact that the
grain attains a higher velocity from the fluid with a larger angle and higher initial velocity.
I"1
Ca
in
( 0.8
E
0
o
Z -,1
Nondimensional Length, L/d
Figure 3.7 Particle Trajectory Sensitivity to Initial Trajectory Angle
*a
I0.9
r
.0)
I<
cZ
0.8
.C
E
' a .
0.7
z
0.5 "
0 0.5 1 1.5 2
Nondimensional Length, L/d
Figure 3.8 Particle Trajectory Sensitivity to Initial Velocity
Another parameter that may be adjusted is the bed roughness. This value is usually
taken to be one or two times the grain diameter. The bed roughness affects the fluid
velocity intensity. Different bed roughness values were run with the model to determine
their effect on the particle's trajectory. The initial velocity used was 2u* for each
component with a particle size of 3.1 cm. The coefficient of lift, CL, equals 0.2, u. =
22.84 cm/s, and s = 2.65. These were chosen to match one of the cases of Nino and
Garcia (1992). Gravity, drag, shear lift, and Basset were the acting forces. The
trajectories are shown in figure 3.9. It shows that the higher the bed roughness, the
shorter the saltation length.
The particle density may also be adjusted; however, since the focus is on sand grains,
the density is held constant at 2.65 g/cm3. This model can also be used for particles other
than sand grains. As a way of showing the effect of particles of different densities under
the same fluid conditions, the model was run with s = 2.65 and s =1.3. As seen in figure
3.10, the less dense the particle, the farther it will saltate.
To summarize, the forces acting on the particle vary significantly with grain size. The
shear lift force and the Basset history force vary the greatest of the forces examined. The
Magnus effect increases saltation length and height greatly as initial angular particle
velocity increases. Finally, the larger the initial velocity and angle, the greater the saltation
length and height.
Random Trajectory Model
The random trajectory model is just an extension of the previous model. Instead of
ejecting only one grain, the random model ejects 42 grains over a specified range of initial
angles. Then for ten different initial velocities, a set of 42 trajectories would exist. The
S1-
0.7 -
0.6
0.5
0 1 2 3 4 5 6 7 8
Nondimensional Length, L/d
Figure 3.9 Comparison Of Variable Bed Roughness Values
2.5
Ss=1.3
S- s =2.65
2-
F!
"-
0
01
c 1.5
z
0.5
0 5 10 15 20 25 30 35 40
Nondimensional Length, L/d
Figure 3.10 Simulated Trajectories of Particles With Different Densities
resulting trajectories could then be used to determine a concentration profile and finally a
bedload model.
Determination of Initial Angles
The angles were determined through the maintenance of equal increments of surface
area over the top half of a sphere. The objective was to have 21 increments of 0 and 21
increments of 0, using standard spherical coordinates. The angle, 0, is the angle between
the x and y axis and the angle, ), is the azimuth angle when a three-dimensional axes is
considered. Since this is a two-dimensional problem, however, only the 21 increments of
4 are considered. This angle ranges from 0 to 7t/2. The following equation is used for
surface area:
AS = sinqAA equation 3.16
where AS is the quasi-surface area because the angle, 0, is absent. This value should
remain constant while 4 changes. Or this equation can be written
sinodo = -d(coso) equation 3.17
which is also a valid expression. Since 0 ranges from 0 to 7c/2, coso ranges from 0 to 1.
Therefore for 21 increments of 0, coso is divided into 21 segments. Table 3.1 displays the
actual initial lift-off angles.
The angle phi can be found by integrating over the range 01 to 02. The following
shows how to obtain this angle.
2Sds equation 3.18
S= sin Odo = cos 1 -cos02 equation 3.18
Since S is a constant, it may be stated that
cos(, -cos2, =sin),(02 e,)
or solving for 0*
equation 3.20
ScosKd, -cos)2,
(. = arcsin -2 ----1
02 -1
where the angle 0* has no physical significance. It is only the angle that satisfies the
surface area equation.
Table 3.1 Distribution of Initial Angle
cos4
0
1/21
2/21
3/21
4/21
5/21
6/21
7/21
8/21
9/21
10/21
11/21
12/21
13/21
14/21
15/21
16/21
17/21
18/21
19/21
20/21
1
(deg)
90
87.27
84.53
81.78
79.02
76.23
73.40
70.53
67.61
64.62
61.56
58.41
55.15
51.75
48.19
44.41
40.37
35.95
31.00
25.21
17.75
0
(deg)
88.42
85.82
83.11
80.37
77.60
74.79
71.95
69.05
66.10
63.08
59.97
56.76
53.44
49.96
46.29
42.38
38.15
33.46
28.09
21.46
8.84
sin,*Ao
0.0476
0.0477
0.0476
0.0475
0.0476
0.0477
0.0476
0.0476
0.0476
0.0477
0.0476
0.0476
0.0476
0.0477
0.0476
0.0477
0.0475
0.0477
0.0476
0.0476
0.0476
----
I I
equation 3.19
;s
The surface area calculations are also shown in table 3.1. The values are nearly the
same and therefore validate this approach. All of these angles are found prior to running
the model.
Boundary Conditions
Once the angles have been determined, the initial velocities are needed. For this
study, the initial velocities were taken to be one to ten times the critical velocity for
incipient motion. That is, for case one, the initial velocity was ucr, for case two, 2ucr, and
so on. These numbers are purely arbitrary Statistically speaking, the grains are given an
equal chance to leave the bed at all of these angles and velocities, i.e., uniformly
distributed.
Besides ejecting with the flow, grains also may eject from the bed initially going
against the flow. This may occur during sheet flow described earlier in the section
"Particle Saltation." Of the models published, this is the first time that a model has
3.5-
S2.5
.
1.5
0.5
-2 0 2 4 6 8 10 12 14
Nondimensional Length, Lb/d
Figure 3.11 Typical Random Trajectories With Initial Velocity 5ucr
Figure 3.11 Typical Random Trajectories With Initial Velocity 5Ucr
1 2.5
0)
Co
0)
E
-D
c 1.5
z
0.5'
-5 0 5 10 15
Nondimensional Length, L/d
Figure 3.12 Particle Trajectories Of Same Initial Velocity And Angles Symmetric About
The Vertical Plane
included grain ejection against the fluid. The angles previously found are translated to the
other plane by making the x-component of the initial velocity negative. This yields the
total of 42 different angles at ten different velocities.
Sensitivity Analysis
A typical result is plotted in figure 3.11. The relevant quantities are d = 0.1 cm, ks =
2d, CL = 0.85CD, 9o = 20 rev/s, and an initial velocity magnitude of 5ucr. Note that the
grains that are ejected into the shear flow at the same angle and velocity as ones that are
ejected with the fluid have greater saltation lengths and heights. This can be more readily
seen in figure 3.12. This figure shows grains ejected from the bed at the same angle with
respect to the bed. The difference exists because of the way the shear lift is defined by
Wiberg and Smith (1985). When the grain is ejected against the flow, it has a negative x-
44
component. This causes the magnitude of the shear force to be much greater than the
force acting on the grain ejected with the flow.
CHAPTER 4
SALTATION MODEL RESULTS
Single Grain Model
Wiberg and Smith (1985) Shear Lift
The Magnus Effect is neglected for the first comparison with a given data set to see
how well the shear lift force predicts the saltation trajectory. The trajectory data used for
the model comparison is provided by Fernandez Luque and van Beek (1976). The shear
lift force that was first tried in the model was Wiberg's and Smith's (1985) lift. From
Fernandez Luque and van Beek (1976), the grain diameter was 0.18 cm with a fluid shear
velocity of 4 cm/s. The initial conditions for the model were uo = 8 cm/s, wo = 8 cm/s, 3 =
0, CL = 0.85CD, and ks = 2d. The trajectory is shown in figure 4.1. It is easily seen that
the shear lift force alone does not adequately describe the motion of the particle as
observed by Fernandez Luque and van Beek (1976).
Figures 4.2 through 4.6 show some of the relevant velocities and forces as predicted
by the model. The fluid velocity is always in the rough, turbulent range as shown by the
Reynolds number. The drag coefficient was between 0.67 and 0.86. Those quantities are
graphically shown in figure 4.2 where time is nondimensionalized with the time it takes for
the grain to complete its trajectory, Ts. Figure 4.3 shows how the grain, fluid, and relative
velocities varied over the trajectory length of the grain saltation. The plot reveals that the
grain velocity is continuously increasing to equal that of the velocity of the fluid. So as a
result, the relative velocity is decreasing. This means that the drag and lift forces are both
decreasing over the grain trajectory. This is shown in figure 4.5. The negative lift results
b
4.5
4
-c
.
13.5
E
I 3-
"o
0
o 2-
z
1.5
1
0 5 10 15 20 25
Nondimensional Length, L/d
Figure 4.1 Model Comparison With Fernandez Luque and van Beek (1976)
3001 1 II1
E 250
z
-o
c200
cc
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Nondimensional Time
Figure 4.2 Reynolds Number and Drag Coefficient As Calculated By The Model
.. Fernandez Luque &
van Beek (1976) Range
- Model, d=0.18cm,ks=2d,
CI=0.85Cd, u*=4cm/s
~I
0 1 2 3 4
Nondimensional Length, L/d
Figure 4.3 Simulated Particle Velocities
C
c)
0
0
. -2
E
Z-4
-6 k.
2 3 4
Nondimensional Length, L/d
Figure 4.4 Simulated Particle Accelerations
-- uf/u*
- ug/u*
..... v/U*
SHorizontal Acceleration
Vertical Acceleration
. . . . . . . . . . . .. ,. . ..
48
0.8- Lift Force, FI/Fg
0.7 Drag Force, Fd/Fg
0D
LL 0.6
LL
a 0.5
0.4
E 0.3 -
CO
.2 0.2
( -
E 0.1
0
z 0
-0.1
-0.2
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Nondimensional Length, L/d
Figure 4.5 Simulated Particle Forces
0.25
0.2
1.
0
L. 0.15-
C:
I- 0.1 -
E
o
z
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Nondimensional Length, L/d
Figure 4.6 Simulated Particle Basset History Force
from the fact that the difference between the fluid velocity at the top of the grain and the
grain velocity is less than the difference between the fluid velocity at the bottom of the
grain and the grain velocity. Also note that there is a "jump" in the lift force at the very
beginning of the curve as presented in the figure. This results from the way the lift force is
defined once again. At first the bottom of the particle has no fluid velocity. It then
immediately acquires a bottom velocity as the particle leaves the bed and causes the spike
in the curve.
The horizontal acceleration of the particle decreases as the particle velocity
approaches the fluid velocity as shown in figure 4.4. The vertical acceleration remains
fairly constant except near the bed. The accelerations are nondimensionalized by
multiplying the acceleration by TJu*. The particle decelerates significantly in the rising part
of the trajectory and continues to do so until the particle approaches the bed. The particle
accelerates again near the bed because of the aforementioned lift force. This shear lift
force is significant when compared to the drag force. The downward aspect of the lift
force is clearly displayed in figure 4.5. It is evident from this figure that the lift force is not
providing enough upward thrust for the particle to saltate as observed by Fernandez
Luque and van Beek (1976). The Basset force is smaller than both the lift force and the
drag force for much of the saltation as demonstrated in figure 4.6. Note that the Basset
force continually increased as the particle saltated.
Saffman (1965) Shear Lift
The lift force described by Wiberg and Smith (1985) was replaced by the Saffman
(1965) lift force to determine which allows the model to match the given trajectories more
closely. The same initial conditions were used as with the other shear lift force case. The
grain diameter was once again 0.18 cm, the grain roughness was two times the grain
diameter with the initial vertical and horizontal velocities both two times the shear
velocity. The shear velocity was 4 cm/s. The lift coefficient was constant at 20.
Figures 4.7 through 4.12 demonstrate the grain's properties with the inclusion of the
Saffman (1965) lift force. From figure 4.7, it is clear that this lift force has more of an
impact on the particle's trajectory than its predecessor. Its predicted trajectory exceeds
that of the range observed by Fernandez Luque and van Beek (1976). The height is
exceeded by a grain diameter while the length is by about five grain diameters. The
Reynolds number and therefore the drag coefficient both change dramatically as the
particle travels.
The grain's velocity approaches and eventually surpasses the fluid velocity near the
end of the trajectory. This occurs because the shear lift force increases near the bed and
causes an actual upward lift near the end of saltation unlike the lift defined by Wiberg and
Smith (1985). The Saffman (1965) lift keeps the grain higher in the fluid column and
allows it to maintain a higher velocity.
5.5- .. Fernandez Luque &
van Beek (1976) Range
5 Model, d=0.18cm,ks=2d,
CI=20, u*=4cm/s
4.5
2-
1.- 4
I(
S3.5
3
CO 3-
E
2.5 ." "
2 /
0 5 10 15 20 25 30 35
Nondimensional Length, L/d
Figure 4.7 Model Comparison With Fernandez Luque and van Beek (1976)
400
. 350
E
z
'300
0-
( 250
r:
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
0.8
0.75
0.7
0.65
0.6
0.55
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Nondimensional Time
Figure 4.8 Reynolds Number And Drag Coefficient As Calculated By Model
14
uf/u*
12 ug/u*
.*. v/u*
10-
.8
0 4
6 \
4" / "" ,
0 4 / "'....
0 5 10 15 20
Nondimensional Length, L/d
Figure 4.9 Simulated Particle Velocities
15
10
5
0
t--
E
< -5
0
-15
Z
-20-
-25-
-30
0 5 10 15 20
Nondimensional Length, L/d
Figure 4.10 Simulated Particle Accelerations
41 1 --1-1
L 3
-a
LL
0)
S2.5
LL
I, 2
o(
C
0
a 1.5
E
0 1
C
25 30 35
0 5 10 15 20 25 30 35
Nondimensional Length, L/d
Figure 4.11 Simulated Particle Forces
Horizontal Acceleration
Vertical Acceleration
. ... ........ ................................................ .
I I
0 5 10 15 20 25 30 35
Nodimensional Length, L/d
Figure 4.12 Simulated Particle Basset History Force
Once again the horizontal component of the grain's acceleration changes for the
majority of the saltation as the grain is decelerating. The vertical acceleration varies
significantly near the bed. It remains constant for the central portion of the trajectory and
is decelerating and accelerating during the rising and falling portion of the trajectory,
respectively. The lift force is significantly larger than the drag force according to figure
4.11. The shear lift force is greatest near the bed since this is where the velocity gradient
across the grain is the largest. Also note that the curve for the lift force does not contain a
spike as it did in the previous case. This is because the Saffman (1965) lift is defined by
the velocity gradient directly instead of using relative difference. The Basset history force
has the same features as with the lift described by Wiberg and Smith (1985). The force
increases with increasing trajectory length.
Adding Magnus Effect
Since the Saffman (1965) shear lift force seems to overpredict the range from
Fernandez Luque and van Beek (1976), the Magnus force is only applied to the case
54
where the Wiberg and Smith (1985) lift force is concerned. Recall that the latter case did
not come close to matching the observed trajectory range. The Magnus force was added
to the equation of motion to help with the correlation. In addition, the moment equation
was also added to be solved simultaneously with the equation of motion. The initial
conditions were maintained with one exception. An initial angular velocity of 33 rev/s was
added. The model comparison with the data range is displayed in figure 4.13. The model
with the newly added Magnus force predicts the grain trajectory fairly well and certainly
much better than without it. The Magnus force significantly increases the overall lift effect
on the particle especially in the rising part according to figure 4.14. For this particular
case, the Magnus force increased the saltation height by 2.3 grain diameters and the length
by 21 grain diameters. With the Wiberg and Smith (1985) form of the lift force, the
Magnus force needs to be included. This is the same conclusion agreed upon by Wiberg
and Smith (1985).
4.5
0 Fernandez Luque &
4 o van Beek (1976)
- Model, 33 rev/s
Model, Shear lift
3.5 .- o
0 Only
3
0- / \o
0
I /
/ \
S2.5 /
E 2
Sl L ,
a) / /
/ \
/ \
/ \
o\
0 5 10 15 20 25 30
Nondimensional Length, L/d
Figure 4.13 Simulated Particle Trajectory With Magnus Effect
1
LL
0.5
o 0 .5 ........................ .. ..........
o 0
z
-0.5
0 5 10 15 20 25 30
Nondimensional Length, L/d
Figure 4.14 Simulated Particle Combined Lift Force Effect
Random Trajectory Model
This model incorporated the use of the Wiberg and Smith (1985) shear lift force and
the Magnus lift force. Recall that this model ejected grains from the bed over a range of
42 angles: 21 against the shear flow and 21 with the shear flow. The initial velocities were
taken to be one to ten times the critical velocity for incipient motion.
Comparison to Previous Data Set
The random model was first run to try and predict the observed saltation range of
Fernandez Luque and van Beek (1976). From their experiment, the shear velocity was set
at 4 cm/s. The magnitude of the initial velocity was set equal to 11.3137 cm/s so that a
particle ejected at 45 degrees with respect to the bed will have a longitudinal and vertical
component of 2u* or 8 cm/s. This was done because White and Schulz (1977) observed
lift-off velocities of 2u. and this initial velocity was used by van Rijn (1984) when he
compared his model to this data. The grain size and bed roughness were 0.18 cm and
0.36 cm, respectively. An initial angular velocity of 20 rev/s was used in the model
simulation. This value was chosen because the trajectory had been matched satisfactorily
with 33 rev/s at 45 degrees. So 20 rev/s was chosen to compensate for using angles larger
than 45 degrees.
The resulting trajectories can be seen in figure 4.15. The saltation paths do fit into
the established experimental region. The trajectories of only the grains ejecting with the
flow were also examined. This allowed for the quantification of the effect of ejecting
grains against the flow. The results are displayed in table 4.1.
Table 4.1 Average Saltation Characteristics For Fernandez Luque And van Beek (1976)
How Grain Ejected Avg. Length (L/d) Avg. Height (H/d) Avg. Velocity (U/u*)
With And Against 20.94 3.75 6.74
Flow
With Flow 17.85 2.96 6.32
The average saltation length and height for Fernandez Luque and van Beek (1976)
were 24.5 and 2.8, respectively. The model run with grains ejected with the flow only was
closer to matching the height while the other combined effect was closer to the length.
The discrepancy may lie in the inaccuracy of the initial angular velocity.
Comparison of Grain Ejection Methods
The random trajectory model was also run by trying several different initial velocities.
As a simplification, the initial velocities were taken to be multiples of the critical velocity
for incipient motion. The objective of this was to generate a set of trajectories with grains
S3.5
t-
I 3
07
C
2.5
E
C 2
o
z
1.5
1
0.5
-5 0 5 10 15 20 25
Nondimensional Length, Lb/d
Figure 4.15 Simulated Trajectories With Grain Ejected With And Against Shear Flow
0 5 10 15 20 25
Nondimensional Length, Lb/d
Figure 4.16 Simulated Trajectories With Grains Ejected With Shear Flow
ejected both with and against the flow and with grains ejected only with the flow to
compare their average saltation characteristics. That is, to determine the difference
between the two cases' average predicted length, height, and velocity over a range of
initial velocities. It is necessary to do this because a model has not been run before with
grains ejected into the fluid shear.
Figure 4.17 displays the results of the first 42 trajectories with an initial velocity of ucr.
Figures 4.18 through 4.26 show the trajectories that result from initial velocity values of
2ucr through O1ucr. In addition an initial velocity of 100ucr is considered to show an
extreme case. The grain diameter for all cases is 0.1 cm since this is the upper limit for
sediment to be considered as sand. The lift coefficient, CL, is equal to 0.85CD. This
means that the Wiberg and Smith (1985) lift force was incorporated into the model. The
equivalent bed roughness is equal to two times the grain diameter. The initial angular
velocity was chosen to be 20 rev/s.
1.1
"-o
"Z 0.9
Co
S0.8
E
0
zn7
-0.5 0 0.5 1 1.5 2 2.5
Nondimensional Length, Lb/d
Figure 4.17 Simulated Trajectories With Initial Velocity, ucr
I 1.4
0
-Il
CO
t-
E 1
0
Z
-1 0 1 2 3 4
Nondimensional Length, Lb/d
Figure 4.18 Simulated Trajectories With Initial Velocity, 2ucr
-1 0 1 2 3 4 5
Nondimensional Length, Lb/d
Figure 4.19 Simulated Trajectories With Initial Velocity, 3ucr
0.51 \1 I A '
-2 0 2 4 6
Nondimensional Length, Lb/d
Figure 4.20 Simulated Trajectories With Initial Velocity, 4ucr
r- I
0.5'
-2 0 2 4 6 8
Nondimensional Length, Lb/d
Figure 4.21 Simulated Trajectories With Initial Velocity, 5ucr
10 12 14
10
A- ------------------ ----- ----- -----
*1 I i I I I
3.5
3
o
T 2.5
"(
S 2
E
0
0-
OI
0.5
-2 0 2 4 6 8 10
Nondimensional Length, Lb/d
Figure 4.22 Simulated Trajectories With Initial Velocity, 6ucr
4.51 I
-o
3
-r,
2.5
C,
.E 2
*"
0
z
-o
c=
z
S 1\\\\\ 4 1I
12 14 16
-5 0 5 10
Nondimensional Length, Lb/d
Figure 4.23 Simulated Trajectories With Initial Velocity, 7ucr
20
4.5
4
S3.5
C
0
o
CO
S2.5
r
0
"Zi .
0.5' I I % '
-5 0 5 10 15
Nondimensional Length, Lb/d
Figure 4.24 Simulated Trajectories With Initial Velocity, 8ucr
5.5
5
4.5
-o
4-
3.5
z 2
O"5 3
E
5 2.5
O
2
1.5
1
0.5
-5 0 5 10 15
Nondimensional Length, Lb/d
Figure 4.25 Simulated Trajectories With Initial Velocity, 9ucr
20 25
25
I
63
7
6
I
0r
a)
14
O
2
1-
-5 0 5 10 15 20 25 30
Nondimensional Length, Lb/d
Figure 4.26 Simulated Trajectories With Initial Velocity, 10ucr
100
90
80
p 70
2-
20--
-5 0 5 10 15 20 25 3060
Nondimensional Length, Lb/d
Figure 4.26 Simulated Trajectories With Initial Velocity, lOucr
100C
C 50
E 40
r.
0
z 30
20
10
-200 0 200 400 600 800 1000
Nondimensional Length, Ud
Figure 4.27 Simulated Trajectories With Initial Velocity, 100ucr
The trajectories of the grains ejected against the flow consistently go higher and
farther than those grains that initially travel with the flow. This is a result of the higher
shear lift force magnitude for those grains shot into the flow. The average saltation
height, length, and particle velocity for each initial velocity are summarized in table 4.2.
An obvious result that may be seen from the table is that both the average saltation height
and length increase with increasing initial velocity.
The average particle velocity is the average longitudinal velocity of the saltating grain.
The velocity has been nondimensionalized by the fluid shear velocity which was 2.317
cm/s for all of the cases. This was equal to the critical velocity for incipient motion for
grain size of 0.1 cm. Like the other averaged quantities, the average particle velocity also
increased with increasing initial velocity. This is an expected result because the grain is
Table 4.2 Average Saltation Characteristics From Random Trajectory Model
Initial Vel. (n*Ucr) Avg. Length (Lb/d) Avg.Height (Hb/d) Avg.Velocity (u/u*)
1 3.18 0.98 3.34
2 4.70 1.35 4.03
3 6.24 1.73 4.61
4 7.85 2.12 5.14
5 9.55 2.53 5.63
6 11.29 2.96 6.09
7 13.22 3.39 6.52
8 15.14 3.83 6.93
9 17.14 4.28 7.32
10 19.20 4.73 7.69
100 526.67 61.59 30.95
allowed to travel higher into the fluid column and therefore acquire more of the fluid's
velocity. As a result of this momentum transfer, the grain travels farther.
These values that are chosen may or may not be physically realistic in grain saltation.
As stated previously, none of the models nor experiments have dealt with the ejecting of
grains into the flow. This makes it difficult to determine if this random trajectory model is
correct in ejecting grains into the fluid flow. If it is accurate to eject grains against the
flow, it may not be at the same velocities as those that are ejected with the fluid. All of
these points need to be looked at to either validate the above approach or modify it.
The saltation characteristic averages for the grains ejected both with and against the
flow are compared with the values obtained by ejecting the grains with the shear flow
only. The initial conditions were the same as those for the previous case. The grains were
Table 4.3 Saltation Characteristics With Grain Ejected With Shear Flow
Initial Vel. (n*Ucr) Avg. Length (L/d) Avg. Height (H/d) Avg. Velocity (u/u*)
1 2.86 0.87 3.20
2 4.35 1.16 3.92
3 5.86 1.46 4.55
4 7.44 1.77 5.15
5 9.17 2.10 5.72
6 10.92 2.44 6.23
7 12.88 2.79 6.82
8 14.87 3.15 7.34
9 16.97 3.52 7.86
10 19.14 3.90 8.36
100 615.50 58.19 49.72
0 1 2 3 4 5 6 7 8 9 10
Initial Velocity, n*ucr
Figure 4.28 Average Length And Height Of Particle Saltation From Random Model
20 1
O Average Length O
18- + Average Height
Average Velocity 0
16
o 0
a)
> 14 -
6 o
.0
12-
.~ 0
T-
a10
i) o
a)
0
a) K
.-o o
E 6 0
C X
w -4.
5 6
Initial Velocity, n*ucr
Figure 4.29 Average Saltation Length And Height From Particles Ejected With Flow
Table 4.3 Saltation Characteristics With Grain Ejected With Shear Flow
Initial Vel. (n*Ucr) Avg. Length (L/d) Avg. Height (H/d) Avg. Velocity (u/u*)
1 2.86 0.87 3.20
2 4.35 1.16 3.92
3 5.86 1.46 4.55
4 7.44 1.77 5.15
5 9.17 2.10 5.72
6 10.92 2.44 6.23
7 12.88 2.79 6.82
8 14.87 3.15 7.34
9 16.97 3.52 7.86
10 19.14 3.90 8.36
100 615.50 58.19 49.72
ejected at velocities ranging from Ucr to 10ucr with the angles varying from 0 to 90 degrees
with respect to the bed. Like before a case for 100 ucr is also examined.
The average characteristics are summarized in table 4.3. It can be seen from the table
that the quantities are smaller by approximately ten percent compared to the values in
table 4.2. This would be expected because of the shear lift force effects on particles
ejected against the flow as described earlier. This is true for all but the 100Ucr case where
the average length and particle velocity are greater than in table 4.2. Still, the difference is
not that great to warrant abandoning this new approach. The process of ejecting grains
into the flow appears to have less of an effect than changing the initial angular velocity by
5 rev/s in the Magnus force as demonstrated in Chapter Three.
CHAPTER 5
CONCLUSIONS
The objective of this study was to evaluate existing saltation models and develop a
new model that would hopefully bridge the gap between theory and observation. The new
model incorporates the best aspects of the previous models. In addition a model was
developed to eject grains from the bed with varying angles and initial velocities under the
assumption of grains having an equal likelihood of leaving the bed at all of those angles.
Finally, the model was compared to experimental observations.
Many of the models examined were quite similar. Most contained the same basic
forces and combined them in different manners. This can be seen in table 5.1. All of the
models include gravity, added mass, and drag forces. Meanwhile, Wiberg and Smith
(1985) developed the only model with a form lift in addition to the shear and Magnus lift
forces. Also, Nino and Garcia (1992) had the only model of those reviewed that included
the Basset history force. The model that was developed in this study included all of the
forces listed in table 5.1 except for the form lift force. It was very flexible in that many
parameters could be changed to adapt to many different situations.
A more complete summary is given below:
1. Five particle saltation models were examined in this study. They were: the
Murphy and Hooshiari (1982), the van Rijn (1984), the Wiberg and Smith (1985), the
Nino and Garcia (1992), and the Lee and Hsu (1994) theoretical models. The models
were compared with data sets either obtained by the researchers themselves or with
trajectories provided by Fernandez Luque and van Beek (1976) or Abbott and Francis
(1977).
Table 5.1 Summary of Forces In Models
Forces Murphy van Wiberg Nino and Lee and Current
and Rijn and Garcia Hsu Model
Hooshiari (1984) Smith (1992) (1994) (1995)
(1982) (1985)
gravity X X X X X X
added mass X X X X X X
drag X X X X X X
shear lift X X X X X
Magnus lift X X X X
form lift X
Basset X X
history____
2. The models did show some ability to predict the observed particle trajectories;
however, each had some discrepancies. The Murphy and Hooshiari (1982) model was
unable to accurately predict saltation height because they ignored forces related to the
relative velocity except for the drag force. The van Rijn (1984) model had both
longitudinal and vertical components of the shear lift even though the shearing was only
in the plane normal to the bed. The Wiberg and Smith (1985) model included a form lift
where the particle needs to rotate a certain number of times for the desired effect. In
addition, their model had the particle's rotation remain constant throughout the saltation.
The Nino and Garcia (1992) model also had the rotation remain constant. The model
developed by Lee and Hsu (1994) had the same problem as van Rijn (1984) and assigned a
different angular velocity for three distinct regions of the particle saltation trajectory.
3. The model developed in this study included gravity, added mass, drag, shear lift,
Magnus lift, and Basset history forces. The shear lift force is the same as defined by
Wiberg and Smith (1985) but can be exchanged for the Saffman (1965) lift. When the
Magnus lift force is applied, the moment equation from White and Schulz (1977) is used
to continually update the particle's rotation.
4. The forces have been assumed to apply for the entire range of Reynolds numbers
in this study.
5. The major problem is obtaining a realistic initial angular velocity of a saltating
particle, which is currently considered an adjustable parameter. The Magnus force
expression is very sensitive to the value of the initial angular velocity
6. The random model assumes that the lift-off angles are uniformly distributed. Also,
it allows the grains to be ejected into the shear flow.
7. The average saltation characteristics for grains ejected against and with the flow
differed from those ejected with the flow by approximately ten percent.
8. An experiment needs to be run with grains being ejected into the oncoming flow,
since this phenomena has not been investigated before. How to do this without disturbing
the flow is still not as yet resolved.
APPENDIX A
FORMULATION OF FORCE COMPONENTS
U -Uf
Figure A. 1 Particle Relative Velocity
The formulation of longitudinal and vertical force components is given here. This
section describes, for example, how the components of the drag force were developed
for the grain saltation model.
Drag Force Component Derivation:
The relative velocity is
Vr = (u-u +w2
equation A.1
From figure A. 1
cos =
u(u-uf) +w2
equation A.2
W
sinO =
(U- f)2 +W2
equation A.3
The drag force is defined as
Fd =CDA pV,2
2
equation A.4
The longitudinal component is
Fd = -CDA plVIVcos
2
1 p(u 2 w(u- r +w2 u -Mu
=-CDA'p (-U) +fw2 (U) +2 W-2 U
2 (Juu-u) +w2
=-CDA 2 (ju-u)2 +w2 (u-u )
equation A.5
The vertical component is
1
F = -CDA plVIVsinO
2
1 p2u- u, 2 22 2 W
=-CDAP (u-uf) +w2 [(u-uf+w2 W
2 V(u-uf) +w
equation A.6
=-CDA 2 pl uufY)2+ W (W)
APPENDIX B
BASSET HISTORY FORCE EVALUATION
The Basset history force integral expression is given here as it appears in the model, a. k.
a., the "brute force" approach.
The Basset history force is defined as
diu
Fb = const dt dt
0 !t-
equation A.7
Let al, a2, a3, etc. represent the particle's acceleration over two time steps with al as the
initial acceleration. For the first time step, At, the Basset history integral may be written
( da = a+ a2 =
00 2 0
= -(a + a2X At- At
= (a, + a2 PA
equation A.8
For the next time step, the Basset history integral may be written
a1 + a2
JAt dt+
If2At-T
a2 +a3DJ
2At 2At-'
.4 --t -,
= a2 2A t_ t -(a 2+ 2At -t 2
l2 )2f2 2At -r 2 2 At
i a,+ a,,,, :~)liZ~~At
=[(a, + a2 X )+(a2 + a)]At
Therefore in summation notation, the term can be expressed as
(du\
dt
Sa + aj+1 (2At) + j + n-j)+ a +a (2At
j=1 2 2
equation A.9
equation A. 10
APPENDIX C
MATLAB PROGRAMS USED FOR MODEL DEVELOPMENT
The following programs were used for the analysis of the saltation models, and for the
development of a new predictive model using the previous available literature.
"arch.m"
%%Determines taucr from Archimedes # (from van Rijn)
%%for traject.m
global s rho g nu d taucr
medes=((s-l)*rho*g/(nu^2))^(1/3)*d;
if medes<4
taunon=0.24/medes;
elseif medes<=10
taunon=0.14*medes^(-0.64);
elseif medes<=20
taunon=0.04*medes^(-0.1);
elseif medes<=150
taunon=0.013*medes^(0.29);
elseif medes>150
taunon=0.058;
end
taucr=(s-l)*gam*d*taunon;
"coeffd.m"
%%Calculation of drag coefficient for traject.m
global Re Cd
if Re<0.1
Cd=24/Re;
elseif Re
Cd=22.73/Re + 0.0903/(Re^2) + 3.69;
elseif Re<10
Cd=29.1667/Re 3.8889/(Re"2) + 1.222;
elseif Re<100
Cd=46.5/Re -116.67/(Re^2) + 0.6167;
elseif Re<1000
Cd=98.33/Re 2778/(Re^2) + 0.3644;
elseif Re<5000
Cd=148.62/Re 4.75e04/(Re^2) + 0.357;
elseif Re<10000
Cd=-490.546/Re + 57.87e04/(Re^2) + 0.46;
elseif Re<50000
Cd=-1662.5/Re + 5.4167e06/(ReA2) + 0.5191;
end
"fluidvel.m"
%%Calculation of velocity profile for traject.m
global d uf z z0 flow fltop flbot Ks
%flow=2.5*uf*log(29.7*((z+0.5*d)/Ks) +1);
flow=2.5*uf*log(z/z0);
%fltop=2.5*uf*log(29.7*((z+d)/Ks) +1);
fltop=2.5*uf*log((z+0.5*d)/z0);
%flbot=2.5*uf*log(29.7*(z/Ks) +1);
if z>0.5*d
flbot=2.5*uf*log((z-0.5*d)/z0);
else
flbot=0;
end
if flbot<0
flbot=0;
end
"motion.m"
%%Details the equation of motion of a saltating grain
%%The forces included are gravity, drag, added mass, shear lift,
%%Basset history, and Magnus lift force.
%%This is the m-file called in traject.m and solved with a
%%fourth or fifth order Runge Kutta.
function veldot = motion(t,vel);
global dudt dwdt s g tf d flow fltop flbot Cl Cd bassetx bassetz dragx
dragz..
admass btu btw mu V Ar nu Cm Re coeff dudz gravx gravz beta vrel
liftx...
liftz magx magz rho utop ubot history drag shear spin saff
vrel=sqrt((vel(1)-flow)^2+vel(2)^2);
Re=vrel*d/nu;
%%Call coeffd.m to determine drag coefficient based on Morsi and
Alexander
coeffd
%Cd=(24/Re)+7.3/(l+sqrt(Re))+0.25;
if saff=='N'
Cl=0.85*Cd; %%or remark this out if C1=0.2
end
gravx=rho*(s-l)*g*V*sin(beta);
gravz=-rho*(s-l)*g*V*cos(beta);
if history=='N'
bassetx=0;
bassetz=0;
end
if drag=='Y'
dragx=-Cd*Ar*0.5*rho*vrel*(vel(1)-flow);
dragz=-Cd*Ar*0.5*rho*vrel*vel(2);
else
dragx=0;
dragz=0;
end
if shear=='Y'
if saff=='Y'
liftx=0;
liftz=-Cl*rho*sqrt(nu)*(d^2)*vrel*sqrt(dudz)*((vel(1)-flow)/vrel);
else
utop = sqrt((vel(l)-fltop)^2 + vel(2)^2);
ubot = sqrt((vel(l)-flbot)^2 + vel(2)^2);
liftx=0;
liftz=(rho/2)*Ar*Cl*(utop^2 ubot^2);
end
end
if shear=='N'
liftx=0;
liftz=0;
end
admass=l+Cm/s;
if spin=='Y'
magx=(pi/8)*(d^3)*rho*vrel*(vel(3)-0.5*dudz)*(vel(2)/vrel);
magz=-(pi/8)*(d^3)*rho*vrel*(vel(3)-0.5*dudz)*((vel(1)-flow)/vrel);
else
magx=0;
magz=0;
end
if spin=='Y'
dudt=(bassetx + dragx + gravx + liftx +magx)/(rho*s*V*admass);
dwdt=(bassetz +gravz + dragz + liftz +magz)/(rho*s*V*admass);
%%bassetz acts in same direction as grave
domegadt=-coeff*(vel(3)-0.5*dudz);
veldot=[dudt;dwdt;domegadt];
else
dudt=(bassetx + dragx + gravx + liftx +magx)/(rho*s*V*admass);
dwdt=(bassetz +gravz + dragz + liftz +magz)/(rho*s*V*admass);
%%bassetz acts in same direction as grav
veldot=[dudt;dwdt];
end
"random.m"
clear
global dudt dwdt tf btu btw s d g flow fltop flbot Cl Cd bassetx bassetz
gravx..
gravz beta dragx dragz admass mu V Ar nu Cm Re coeff dudz vrel liftx
liftz...
z z0 rho magx magz Ks utop ubot uf history drag shear spin taucr t
ang=[17.75 7.46 5.79
2.92...
2.87 2.83 2.79 2.76
ang=(pi/180)*ang;
4.95 4.42 4.04 3.78 3.56 3.4 3.26 3.15 3.06 2.99
2.75 2.74 2.73];
history='Y';
drag='Y';
shear='Y';
spin='Y';
beta=0;
deltat=1/500;
s=2.65;
d=0.136;
g=981;
nu=0.009796;
Cm=0.5;
rho=0.9985;
mu=nu*rho;
Ar=(pi/4)*d^2;
V=(pi/6)*d^3;
const=(6*pi*mu*
results=[];
yous=[];
((d/2)^2))/sqrt(pi*nu);
%Calculation of velocity profile
gam=rho*g; %%(g/(cm^2-s^2)
%taucr=0.058*gam*(s-l)*d; %%(g/cm-s^2) --From Shield's diagram (0.058)
%%Call arch.m to determine taucr from Archimedes #
arch
ucr=sqrt(taucr/rho); %(cm/s)
uf=ucr; %%Make friction velocity a multiple
%%of the critical velocity
%%Scaling factors
vscale=sqrt((10/3)*(s-l)*g*d);
tscale=(s+Cm)*sqrt((10/3)*d/(s-l)/g);
lscale=(s+Cm)*(10/3)*d;
ufnon=uf/vscale; %%Nondimensional shear velocity
Ks=2*d;
z0=0.11*(nu/uf) + 0.033*Ks;
%%Initial grain rotation
inert=(2/5)*(rho*s)*V*(d/2)^2;
coeff=pi*mu*(d^3)/inert;
for n=1:10 %%This loop is for 10 cases
initnon=n*ufnon;
init=initnon*vscale; %%Dimensional initial grain velocity
loop=l; %%This loop is for firing grains into and with flow
for loop=l:2
for m=l:21 %%This loop is for 21 angle increments (0-90)
t0=0;
tf=deltat; %%Step size
if loop==l
if m==l
u0=sqrt((init^2)/(1+(tan(ang(m)))^2));
w0=abs(u0)*tan(ang(m));
else
direct=sum(ang(l:m));
u0=sqrt((init^2)/ (+(tan(direct))^2));
w0=abs(u0)*tan(direct);
end
elseif loop==2
if m==l
u0=-sqrt((init^2)/(1+(tan(ang(m)))^2));
w0=abs(u0)*tan(ang(m));
else
direct=sum(ang(l:m));
u0=-sqrt((init^2)/(1+(tan(direct))^2));
w0=abs(u0)*tan(direct);
end
velO = [uO wO 2*pi*20]'; % initial conditions
uvel(l)=vel0(1);
z=0.5*d;
%%Call fluidvel.m to determine fluid velocity
fluidvel
fluid(l)=flow;
A=[];
i=2;
jones=0;
%%Determination of relative velocity
vrel=sqrt((vel0(1)-flow)^2+vel0(2)^2);
%%Determination of Reynolds number
Re=vrel*d/nu;
%%Call coeffd.m to determine drag coefficient
%%based on Morsi and Alexander(1972)
coeffd
%%Set lift coefficient
Cl=0.85*Cd;
%%Initial gravity force/buoyancy force
gravx=rho*(s-1)*g*V*sin(beta);
gravz=-rho*(s-1)*g*V*cos(beta);
uO = ? wO = ?
%%Initial drag force
dragx=-Cd*Ar*0.5*rho*vrel*(velO(1)-flow);
dragz=-Cd*Ar*0.5*rho*vrel*vel0(2);
%%Initial shear lift force
utop=sqrt((vel0(1)-fltop)^2 + velO(2)^2);
ubot=sqrt((vel0(1)-flbot)^2 + vel0(2)^2);
liftx=0;
liftz=(rho/2)*Ar*C1*(utop^2 ubot^2);
%%Initial Magnus force
magx=(pi/8)*(d^3)*rho*vrel*(velO(3)-0.5*dudz)*(vel0(2)/vrel);
magz=-(pi/8)*(d^3)*rho*vrel*(vel0(3)-0.5*dudz)*((vel0(1)-...
flow)/vrel);
%%Initial added mass
admass=l+(Cm/s);
%%Initial Basset force is zero
bassetx=0;
bassetz=0;
%%Initial grain accelerations
dudts(1)=(bassetx+dragx+gravx+liftx +magx)/(rho*s*V*admass);
dwdts(l)=(bassetz +gravz + dragz + liftz...
+magz)/(rho*s*V*admass);
%%Initial grain position
dep(l)=0.5*d;
x=0;
len(l)=0;
time(l)=0;
%Calculation of velocities until grain hits bed
while jones==0
[t,vel] = ode45('motion',t0,tf,vel0);
B=[t vel];
z = 0.5*deltat*(vel0(2)+vel(length(vel),2))+z;
x = 0.5*deltat*(vel0(1)+vel(length(vel),l))+x;
len(i)=x;
dep(i)=z;
time(i)=(i-l)*deltat;
if z<0.5*d
jones=l;
end
%%Keep accelerations tabbed
dudts(i)=dudt;
dwdts(i)=dwdt;
%%Update of Basset history force
btu=0; %% Need these conditions
btw=0; %% Need these conditions
if i==2
btu=(dudts(i-l)+dudts(i))*sqrt(deltat);
btw=(dwdts(i-l)+dwdts(i))*sqrt(deltat);
else
for k=l:(i-2)
btu=(dudts(k)+dudts(k+l))*sqrt(deltat)*...
(sqrt(i-k)-sqrt(i-l-k))+btu;
btw=(dwdts(k)+dwdts(k+l))*sqrt(deltat)*...
(sqrt(i-k)-sqrt(i-l-k))+btw;
end
btu=(dudts(i-l)+dudts(i))*sqrt(deltat)+btu;
btw=(dwdts(i-l)+dwdts(i))*sqrt(deltat)+btw;
end
%%Call fluidvel.m to determine fluid velocity
fluidvel
fluid(i)=flow;
%%Update initial conditions
vel0=[vel(length(vel),1) vel(length(vel),2)...
vel(length(vel),3)]';
t0=tf;tf=tf+deltat;
A=[A;B];
uvel(i)=vel0(1);
%%Update Basset history force
bassetx=-const*btu;
bassetz=-const*btw;
i=i+l;
end
%%Plot of particle trajectory
figure(n)
len=len./d; scalel;
dep=len./d; scalel;
plot(len,dep)
hold on
time=time./max(t); %tscale;
vg=zeros(length(len),l)+n;
hey=[vg uvel' fluid'];
yous=[yous;hey];
interim=[vg time' len' dep'];
results=[results;interim];
nines3=[9 9 9];
nines4=[9 9 9 9];
results=[results;nines4];
yous=[yous;nines3];
clear len dep vel B
end
m=l;
end
hold off
xlabel('Nondimensional Length, L/d')
ylabel('Nondimensional Height, H/d')
end
save results
save yous
"traject.m"
%%This m-file determines the trajectory and forces
%%acting on a single particle in saltation
clear
global dudt dwdt tf btu btw s d g flow fltop flbot Cl Cd bassetx bassetz
gravx...
gravz beta dragx dragz admass mu V Ar nu Cm Re coeff dudz vrel liftx
liftz...
z z0 rho magx magz Ks utop ubot uf history drag shear spin taucr saff
t0=0;
deltat=1/500;
tf=deltat; %%Step size
%%Grain size
d=input('Input the grain size, (cm) ');
%%Initial velocities
init=input('Input the initial velocity, (no units) ');
ang=input('Input the initial take-off angle, (0-180 deg) ');
jones=0;
s=2.65;
g=981; %cm/s^2
nu =0.01; %%0.01; %cm^2/s
Cm=0.5;
rho=l; %1; %g/cm^3
mu=nu*rho; %g/(cm-s)
Ar = (pi/4)*d^2;
V = (pi/6)*d^3;
const=(6*pi*mu*((d/2)^2))/sqrt(pi*nu);
%Calculation of shear velocity
gam=rho*g; %%(g/(cm^2-s^2)
%%Call arch.m to determine taucr
arch
%taucr=0.058*gam*(s-l)*d; %11.1 %%(g/cm-s^2) --From Shield's diagram
(0.058)
ucr=sqrt(taucr/rho); %(cm/s)
uf=ucr; %Make friction velocity a multiple
%of the critical velocity
%%Nondimensional scales (Jenkins)
vscale=sqrt((10/3)*(s-l)*g*d);
tscale=(s+Cm)*sqrt((10/3)*d/(s-l)/g);
lscale=(s+Cm)*(10/3)*d;
%%Transform nondimensional velocity to cgs units
ucrp=ucr/vscale;
init=init*ucrp;
init=init*vscale;
%%Determine velocity components
if ang<=90
ang=pi/180*ang;
xdot=sqrt((init^2)/(l+(tan(ang))^2));
else
ang=180-ang;
ang=pi/180*ang;
xdot=-sqrt((init^2)/(l+(taan(ang))^2))
end
zdot=abs(xdot)*tan(ang);
%%Bed slope
beta=input('Input the bed slope, (rise/run) ');
%%Forces to include in saltation model
history=input('Include Basset history force? Y/N, [Y] ','s');
if isempty(history)
history='Y';
drag=input('Include drag force? Y/N, [Y] ','s');
if isempty(drag)
drag='Y';
end
shear=input('Include shear lift force? Y/N, [Y] ','s');
if isempty(shear)
shear='Y';
end
if shear=='Y'
saff=input('Saffman shear lift force? Y/N, [Y] ','s');
if isempty(saff)
saff='Y';
end
end
spin=input('Include Magnus lift force? Y/N, [Y] ','s');
if isempty(spin)
spin='Y';
end
if spin=='Y'
rev=input('Input initial spin rate of grain, (rev/s) ');
vel0=[xdot zdot 2*pi*rev]';
else
rev=0;
vel0=[xdot zdot];
end
rough=input('Input bed roughness as multiple of grain diameter ');
Ks=rough*d;
z=0.5*d;
z0=0.11*(nu/uf) + 0.033*Ks;
%%Call fluidvel.m to determine fluid velocity
fluidvel
%%Initial fluid velcity w/respect to center of grain
fluid(1)=flow;
%%Initial relative velocity at top of grain
top(1)=fltop;
%%Initial relative velocity at bottom of grain
bot(1l)=flbot;
%%Calculation of velocity gradient across grain
dudzs(1)=dudz;
%%Determination of relative velocity
vrel=sqrt((velO(1)-flow)^2+vel0(2)^2);
Vr(l)=vrel;
%%Determination of Reynolds number
Re=vrel*d/nu;
%Cd=(24/Re)+7.3/(l+sqrt(Re))+0.25;
%%Call coeffd.m to determine drag coefficient
%based on Morsi and Alexander (1972)
coeffd
Rno(1)=Re;
Cdrag(l)=Cd;
%%Set lift coefficient
if saff=='Y'
Restar=uf*d/nu;
if Restar<=5
C1=1.6;
elseif Restar>=70
C1=20;
else
Cl=interpl([5 70],[1.6 20],Restar);
end
end
if saff=='N'
%C1=0.2;
C1=0.85*Cd;
end
Clift(1)=C1;
%%Initial grain rotation
inert=(2/5)*(rho*s)*V*(d/2)^2;
coeff=pi*mu*(d^3)/inert;
if spin=='Y'
larsson(l)=vel0(3);
end
A=[];
i=2;
%%Initial Basset force
bassetx=0;
basx(1)=bassetx;
bassetz=0;
basz(1)=bassetz;
%%Initial gravity force/buoyancy force
gravx=rho*(s-l)*g*V*sin(beta);
gravz=-rho*(s-1)*g*V*cos(beta);
%%Initial drag force
if drag=='Y'
dragx=-Cd*Ar*0.5*rho*vrel*(velO(1)-flow);
dragz=-Cd*Ar*0.5*rho*vrel*velO(2);
end
if drag=='N'
dragx=0;
dragz=0;
end
%%Initial shear lift force
if shear=='Y'
if saff=='Y'
liftx=0;
liftz=-Cl*rho*sqrt(nu)*(d^2)*vrel*sqrt(dudz)*((vel0(1)-
flow)/vrel);
else
utop=sqrt((vel0(l)-fltop)^2 + velO(2)^2);
ubot=sqrt((vel0(1)-flbot)^2 + velO(2)^2);
utops(1)=utop;
ubots(1)=ubot;
liftx=0;
liftz=(rho/2)*Ar*Cl*(utop^2 ubot^2);
end
end
if shear=='N'
liftx=0;
liftz=0;
end
%%Initial Magnus force
if spin=='Y'
magx=(pi/8)*(d^3)*rho*vrel*(velO(3)-0.5*dudz)*(velO(2)/vrel);
magz=-(pi/8)*(d^3)*rho*vrel*(velO(3)-0.5*dudz)*((velO(1)-flow)/vrel);
end
if spin=='N'
magx=0;
magz=0;
end
magsx(1)=magx;
magsz(1)=magz;
%%Initial added mass
admass=l+(Cm/s);
%%Initial grain accelerations
dudts(1)=(bassetx+dragx+gravx+liftx +magx)/(rho*s*V*admass);
dwdts(1)=(bassetz +gravz + dragz + liftz +magz)/(rho*s*V*admass);
%%Initial grain position
dep(1)=z; % z=0 from velocity profile calculation
x=0;
len(l)=x;
time(1)=0;
%Calculation of velocities until grain hits bed
while jones==0
[t,vel] = ode45('motion',t,0,tf,vel);
B=[t vel];
%%Calculation of grain position
z = 0.5*deltat*(vel0(2)+vel(length(vel),2))+z;
x = 0.5*deltat*(velO(1)+vel(length(vel),l))+x;
len(i)=x;
dep(i)=z;
time(i)=(i-l)*deltat;
%%Check if grain reached bed
if z<0.5*d
jones=l;
end
%%Keep accelerations tabbed
dudts(i)=dudt;
dwdts(i)=dwdt;
%%Update of Basset history force
btu=0; %% Need these conditions
btw=0; %% Need these conditions
if i==2
btu=(dudts(i-l)+dudts(i))*sqrt(deltat);
btw=(dwdts(i-l)+dwdts(i))*sqrt(deltat);
else
for k=l:(i-2)
btu=(dudts(k)+dudts(k+l))*sqrt(deltat)*...
(sqrt(i-k)-sqrt(i-l-k))+btu;
btw=(dwdts(k)+dwdts(k+l))*sqrt(deltat)*...
(sqrt(i-k)-sqrt(i-l-k))+btw;
end
btu=(dudts(i-l)+dudts(i))*sqrt(deltat)+btu;
btw=(dwdts(i-l)+dwdts(i))*sqrt(deltat)+btw;
end
%%Call fluidvel.m to determine fluid velocity
fluidvel
%%Fluid velocity at middle of sphere
fluid(i)=flow;
%%Relative velocity at top of grain
top(i)=fltop;
%%Relative velocity at bottom of grain
bot(i)=flbot;
dudzs(i)=dudz;
%%Update initial conditions
if spin=='Y'
vel0=[vel(length(vel),1) vel(length(vel),2)...
vel(length(vel),3)] ';
else
vel0=[vel(length(vel),1) vel(length(vel),2)]';
end
t0=tf;tf=tf+deltat;
A=[A;B];
%%Keep record of contributing forces
bassetx=-const*btu;
bassetz=-const*btw;
basx(i)=bassetx;
basz(i)=bassetz;
Vr(i)=vrel;
Rno(i)=Re;
Cdrag(i)=Cd;
Clift(i)=C1;
if shear=='Y'
utops(i)=utop;
ubots(i)=ubot;
end
if spin=='Y'
larsson(i)=vel0(3);
end
magsx(i)=magx;
magsz(i)=magz;
i=i+l;
end
%%Plot of particle trajectory
p=l;
for q=l:18:length(A)
ts(p)=A(q,l);
us(p)=A(q,2)
ws(p)=A(q,3);
p=p+1;
end
ts(p)=A(length(A),1);
ws(p)=A(length(A),3);
u=us/uf;
figure(1)
plot(len./d,dep./d)
xlabel('Nondimensional Length, L/d')
ylabel('Nondimensional Height, H/d')
figure(2)
subplot(2,1,1)
plot(time./max(ts),Rno)
ylabel('Reynolds Number')
subplot(2,1,2)
|
RSS
HIDE
