Citation
Evaluation and development of particle saltation models

Material Information

Title:
Evaluation and development of particle saltation models
Series Title:
UFLCOEL-95012
Creator:
Krecic, Michael R., 1970-
University of Florida -- Coastal and Oceanographic Engineering Dept
Publication Date:
Language:
English
Physical Description:
xiii, 93 leaves : ill. ; 29 cm.

Subjects

Subjects / Keywords:
Dissertations, Academic -- Coastal and Oceanographic Engineering -- UF ( lcsh )
Genre:
bibliography ( marcgt )
non-fiction ( marcgt )

Notes

Thesis:
Thesis (M.S.)--University of Florida, 1995.
Bibliography:
Includes bibliographical references (leaves 90-92).
General Note:
Typescript.
General Note:
Vita.
Statement of Responsibility:
by Michael R. Krecic.

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
33662099 ( OCLC )

Full Text
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 ........................................................................................................ A
LIST OF SYM BOLS ...................................................................................................... ix
CHAPTER 1: INTRODUCTION .................................................................................... I
CHAPTER 2: PREVIOUS W ORK .................................................................................. 4
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 otion ................................................................................................. 27
Boundary Conditions ............................................................................................... 31
M ethod of Solution ................................................................................................. 32
Sensitivity Analysis .................................................................................................. 32
Random Trajectory M odel .......................................................................................... 38




D eterm ination of Initial Angles ................................................................................ 40
Boundary Conditions ............................................................................................... 42
Sensitivity Analysis .................................................................................................. 43
CH APTER 4: SAL TA TION M OD EL RESULTS ......................................................... 45
Single G rain M odel ..................................................................................................... 45
W iberg and Sm ith (1985) Shear Lift ........................................................................ 45
Saffm an (1965) Shear Lift ....................................................................................... 49
A dding M agnus Effect ............................................................................................. 53
Random Trajectory M odel .......................................................................................... 55
Com parison to Previous D ata Set ............................................................................ 55
Com parison of G rain Ejection M ethods ................................................................... 56
CH A PTER 5: CON CLU SION S .................................................................................... 68
APPENDIX A: FORMULATION OF FORCE COMPONENTS ................................... 71
APPENDIX B: BASSET HISTORY FORCE EVALUATION ...................................... 73
APPENDIX C: MATLAB PROGRAMS USED FOR MODEL DEVELOPMENT ........ 75 REFEREN CES .............................................................................................................. 90
BIO GRAPH ICAL SKETCH ......................................................................................... 93




LIST OF TABLES

Table Pa e
3.1 D istribution of Initial A ngles .................................................................................... 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 Sum m ary 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 T, 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 efinition of B ed Slope ............................................................................................ 29
3.3 The Saltating G rain 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 5u ............................................. 42
3.12 Particle Trajectories Of Same Initial Velocity And Angles Symmetric About The V ertical P lane ................................................................................................................ 4 3
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 Simulated 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 V elocities .................................................................................... 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, Ucr ...................................................... 58
4.18 Simulated Trajectories With Initial Velocity, 2uc ..................................................... 59
4.19 Simulated Trajectories With Initial Velocity, 3u ..................................................... 59




4.20 Simulated Trajectories With Initial Velocity, 411cr........................................... 60
4.21 Simulated Trajectories With Initial Velocity, 5u ....................................... 60
4.22 Simulated Trajectories With Initial Velocity, 6u1 ......cr.................................. 61
4.23 Simulated Trajectories With Initial Velocity, 711cr........................................... 61
4.24 Simulated Trajectories With Initial Velocity, 8u....cr ....................................... 62
4.25 Simulated Trajectories With Initial Velocity, 911cr ........................................... 62
4.26 Simulated Trajectories With Initial Velocity, 1011cr .........................................63
4.27 Simulated Trajectories With Initial Velocity, 10011cr.................................... 63
A. 1 Particle Relative Velocity ................................................................ 71




LIST OF SYMBOLS

A cross-sectional area of particle
Ar Archimedes number
b proportionality constant
cl 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
FMagnus Magnus lift force
Fshear shear lift force
g acceleration due to gravity
ge effective gravity
H saltation trajectory maximum height
I moment of inertia
k, 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
T, time particle travels to complete saltation
u saltating particle longitudinal velocity
Ucr critical velocity for incipient motion
Uf fluid velocity
Uffop 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
z0 zero fluid velocity level
P3 bed angle
0 horizontal plane angle
0C nondimensional Shields' parameter
v fluid kinematic viscosity
fluid dynamic viscosity particle angular velocity




Q0 initial particle angular velocity
TI1 flow angle of attack
T-2 related to airfoil centerline curvature
p density of fluid
PS density of saltating particle
Tb bed shear stress
Tcr critical bed shear stress
0o azimuth angle
01 arbitrary surface angle
02 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 grainto-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
p ( b-1g equation 2.1
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 0~, values versus grain Reynolds number (Re.),
R e U equation 2.2
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 p P g "2 d equation 2.3
where p, is the density of the particle.
00
,
ci
: 0~
U)
10-2 0 2
100 10' 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 (195 8) 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,
F, = (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
F= Vduf + (mduf _du"1 equation 2.5
1 dt (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 = -P du 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 FCDA1pV2 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 ufBot and down if ufBot is greater than UTop.
Uf (top) Uf (top)

uf (bot)

uf (bot)

Figure 2.2 Shear Lift Results From Velocity Gradient Across Grain




The shear lift force may be written as
F, (shear)= CLA p(U2op -2 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
1 du ~
C~j2 21 fI
FL (shear)= C oV d d- equation 2.9
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
67td 2 Tdu duf
F f ~i dt dt dt equation 2. 10
0 S
where i is the dynamic viscosity, T, is the total saltation time, and t is a dummy variable (Mci, 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.
or (top) ur (top)
grain N\I Lgrain } jF
of (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
FL (Ma gnus)= i Qd3pV Q-1du equation 2.11
(Mgu)=8 \ 2 dz
with Q as the angular velocity of the grain with units in radls (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
Momet =dQ =-7cd' Q- Idufequation 2.12 dt y2dz )
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,
ge = 2g(ps p) equation 2.13
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= 0 equation 2.14
2g,
for length
L -2uw equation 2.15
and for duration
02w equation 2.16
where uo is the initial horizontal velocity of the sphere and w0 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




AL I tanh

equation 2.17

b=240p(p, p)
(2p, + p)2

equation 2.18

C + gd )--2
21 2W 0+ w

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.5uln(ZI
.z0( )

equation 2.20

where z is the vertical position of the center of the sphere and z0 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

where




their obtained data set. After doing this, they believe the lengths of the trajectories are matched fairly well.

2.5

0 W
0
WO
W
0
* Experiment
* Model

2
E Z 1.5 PM a)
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 Rib (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




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, k~, 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).




E2.5 'I
,- 2
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




FL (Magnus)= 27c- i 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) = 4itpvc2 sin(N11 +- J2 ) equation 2.23
where c2 is proportional to the length of the centerline of the airfoil, T, 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 TP2 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).




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




C.
0)
-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
c
0
C
0.
z 0.8

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: l)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.

E1.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 mmn with a shear velocity of 3.6 cmls, CL equal to 20, and a k, 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 (FMagnus). 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:




du 1
pSV =-CDA2 u-u) +w2(u-u)
dd 2 du
-PVCM du 62 2 dt d
*dt o~l 0
nt3 2 1duw
+-d dp (ufu, + W 22 +z 4(u-u, ~+w 8 (2 f)) (U Uf) + W
+p(s 1)gV sin 3
and
dw 1 / 2
pVdt =-CDA P(u-u) + w2(w) dt 2
(ddw 6t dw
d 2 -dtdc
PVCM dt drXF
7td 3 l(u _uf + w U u
8 2 dz (u2)+w
+-~ACLu~TP -u~0)- ~s 1)~c u -13

+ AC, (Uirop -Uno,)- p(s-1)gV cos
where where

equation 3.1

equation 3.2

2 o(Up fiop )2 2 UATop "- W2

equation 3.3

2 ( f2 e2 UABo U ) U fBot W

equation 3.4




The variables, urrop and ufot, 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(she F1(Magnus)
Vr r F
tOI F d & ,
wO Fb
UO0.5d
Figure 3.1 Force Definition Sketch of Grain Saltation
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

equation 3.5

FL = -Cpuo.d2u f)2 +W2 duf U-u f dz )f (U f)2 + W

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 Re

22.73 0.0903
C + + 3.69
CD Re Re2
29.1667 3.8889
CD= +1.222
D Re Re2
46.5 116.67
CD = e- +0.6167
DRe Re2
98.33 2778
CD + 0.3644
CD Re Re2
148.62 4.75e04
CD = Re- 0.357
D Re Re2

for Re<0. 1

for 0.1
for 1000
equation 3.6

equation 3.7 equation 3.8 equation 3.9 equation 3.10

equation 3.11




D -490.546 57.87e04 0.46 for 5000 CD Re Re2
-1662.5 5.4167e06
CD + +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
z0 = 0.11k-- + 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




Figur 3.3The altaing Gain nitil Poitio
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 CT 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 I
0 05 1 15 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 1P = 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
1.2 .
0- /
Z/
0 1 2 3 4 5 6 7
Nondimensional Length, L/d
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 k,=2d, Cl=O.2, u.=4 cmls, u0=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.




-a)
05 1 6 7 8 9 1
NodmninlLntU
Fiur 36Pril rjcoisSniiiyt atceRtto
AssaeErvosy h nta odtos ftemdlgetyefc h eutn
satto rjcoy h nta nl n eoiyo h ri stenx rao neet Agi0h ri ieue s0 8c.Fgr ipasterslso ayn h ae
Fniurel.6tPartce Toghrajectorntvy ToesParetic rotio syfo h attah
offanglestand aintiinhecntatinta velocity. Thmtefudihalre nitl angered 30,ia 45,landt60




.
0.9
r
in
S0.8
E
o
Z

Nondimensional Length, L/d Figure 3.7 Particle Trajectory Sensitivity to Initial Trajectory Angle

I 0.9
LM
c 0.8
0
.n
c
E c 0.7
z

0.5 I I I
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/cm'. 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 =13. 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




- 1.2
0.9
. 01.18
0.7
ci)
0.6
0.9
C
09
Z 0.8
0.7 0.6
0.5
0 1 2 3 4 5 6
Nondimensional Length, L/d
Figure 3.9 Comparison Of Variable Bed Roughness Values
2.5
Ss=1.3
. . . . . .
- s = 2.65 2
72
O
- 1
.
ci)
E
:
0
z
0.5
0 5 10 15 20 25 30
Nondimensional Length, Lid
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, 0, is the azimuth angle when a three-dimensional axes is considered. Since this is a two-dimensional problem, however, only the 21 increments of 0 are considered. This angle ranges from 0 to 7t/2. The following equation is used for surface area:
AS = sinOA0 equation 3.16
where AS is the quasi-surface area because the angle, 0, is absent. This value should remain constant while changes. Or this equation can be written
sinodo = -d(coso) equation 3.17
which is also a valid expression. Since C0 ranges from 0 to 7/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 41 to 02. The following shows how to obtain this angle.
2
S f sin odo = cos o1 COS02 equation 3.18

Since S is a constant, it may be stated that




cosO1 -cos2 = sino(d2 -01)

or solving for 0.

equation 3.20

,arcsin cos1O -cos02
02-01

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 COSO

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

sinoAo

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

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 ur, 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
3
2.5
j 2.5
0
I
- c 2
.2.9
Cl) E
c 1.5
0 z
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




2 2.5 2)
2

E
cc 1.5
z

0.5 1 1
-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, 0 = 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 u0 = 8 cm/s, w0 = 8 cm/s, P3 = 0, CL = 0.85CD, and k, = 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, T. 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




4.5 .... Fernandez Luque
van Beek (1976) R
4 Model, d=0.18cm,k
CI=0.85Cd, u*=4c
3.5
S.
2.5
() E
0 2
Z
... ~ ~~............... . . .
1.5
0.5
0 5 10 15 20 25
Nondimensional Length, L/d
Figure 4.1 Model Comparison With Fernandez Luque and van Beek (1976)
300 .. ...

E 250
n
Z
20
c,200 (D
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

ange ks=2d, m/s

&
{ange
ks=2d, /s




0 1 2 3 4
Nondimensional Length, L/d Figure 4.3 Simulated Particle Velocities

C
0 (
0 .,
C
.2-2
E
"o o-

-6 :

2 3 4
Nondimensional Length, L/d

Figure 4.4 Simulated Particle Accelerations

- uflu*
- ug/u* .... v/u*

- Horizontal Acceleration
Vertical Acceleration




48
0.8 Lift Force, FI/Fg
0.7 .. Drag Force, Fd/Fg
L 0.6
LL
c6 0.5
U
U
S0.4
E2 0.1 ...........
E0 .1
0
z0
.o 0.3
.- 0.1
U
-0.2
II I I I I I I I I
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
0)
..
u.. 0.15
C
o
L 0.15
U/)
(
E
~0 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 Tim. 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 cmls. 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=O.18cm,ks=2d, C1=20, u*=4cm/s
4.5
f- 4
- 3.5
C
.)
E
V2.5 .
0
Z "
2 .
1 .. . . . . . . .
0.5 -"
0 5 10 15 20 25 30 35
Nondimensional Length, Lid Figure 4.7 Model Comparison With Fernandez Luque and van Beek (1976)




400
.6350
E
z
(D 250
200
0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

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
" o
75
C4
a) _0
0
/ 4

0 5 10 15 20
Nondimensional Length, L/d Figure 4.9 Simulated Particle Velocities

-




15
10
5
C 0
S-5
0
o
c -15O
Z
*a -10
-20
E
C -15
0
z
-20
-25
-30,
0 5 10 15 20
Nondimensional Length, L/d
Figure 4.10 Simulated Particle Accelerations
4
3.5
0)
LL 3-- 2.5
LL oi .
0 0
c 1.5
) E
0 C
0.5 .

25 30 35

u0 5 10 15 20 25 30 35
Nondimensional Length, L/d Figure 4.11 Simulated Particle Forces

Horizontal Acceleratio Vertical Acceleration

I r

n




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- van Beek (1976)
- Model, 33 rev/s
3.5- 0 Model, Shear lift
0 Only
"I /
/ \/
D)j /
3- /
/ 0
C2.5 /
0
.o /\
I1) k
E0 2- /
0/
o\
Z /\
1.5 9 /
0.
0 5 10 15 20 25 30
Nondimensional Length, L/d
Figure 4.13 Simulated Particle Trajectory With Magnus Effect




0) IL
.......
o 0 .5 .. .............. .. ..........
(D
" 0.
0E 0
0 z
-0.5
-1 I
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 Lugue And van Beek (1976)
How Grain Ejected Avg. Length (Lid) Avg. Height (Hid) Avg. Velocity (Uiu*)
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




-3.5

" 2.5
D
C 0
E
C- 2
0
z
1.5
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 2Ur through IUcr. 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
-
"3 0.9
C 0
= 0.8
E
0
Z r)7

-0.5 0 0.5 1 1.5 2 2.5
Nondimensional Length, Lb/d
Figure 4.17 Simulated Trajectories With Initial Velocity, Ucr




_0
S1.4
E 1
0
Z
0.8
0.6
-1 0
Figure 4.18 Simulated
2 1.2
C
c
ElI
C 0 z
0.8
-l.0
Figure 4.18 Simulated
2.2i

Z2
1.6
-D
-1.4
0
in
0 z 12
E
zi

1 2 3
Nondimensional Length, Lb/d Trajectories With Initial Velocity, 2ucr

0.8 E

-1 0 1 2 3 4 5
Nondimensional Length, Lb/d
Figure 4.19 Simulated Trajectories With Initial Velocity, 3ucr

6 7 8




0.5
-2 0 2 4 6
Nondimensional Length, Lb/d
Figure 4.20 Simulated Trajectories With Initial Velocity, 4Ucr

A. Fl.

0.5' 1
-2 0 2 4 6 8
Nondimensional Length, Lb/d
Figure 4.21 Simulated Trajectories With Initial Velocity, 5ucr

10 12 14
10 12 14

10

A




-o3
0)
"T 2.5
E 2
E
*0 Z C; 0-

-2 0 2 4 6 8 10 12 14
Nondimensional Length, Lb/d
Figure 4.22 Simulated Trajectories With Initial Velocity, 6Ucr
4.5


f3 a)
2.5
0
C
.E 2

0
z

-5 0 5 10
Nondimensional Length, Lb/d Figure 4.23 Simulated Trajectories With Initial Velocity, 7ucr

20




5
4.5
4
- 3.5
(D I 3
o
0
5 2.5
c2
E
Z
z

0.51 v N \ \1 \ \ \\A \\\ \\\\\
-5 0 5 10 15
Nondimensional Length, Lb/d
Figure 4.24 Simulated Trajectories With Initial Velocity, 8Ucr
5.5
5
4.5
E
-
.
._)
S.5

O- 3
E
l 5 2.5
C 0
z 2
1.5
1
0.5
-5 0 5 10 15
Nondimensional Length, Lb/d
Figure 4.25 Simulated Trajectories With Initial Velocity, 9ucr

25

I




63
7
6
14
0
Z U)
E3
0 z
2
1
-5 0 5 10 15 20 25 30
Nondimensional Length, Lb/d
Figure 4.26 Simulated Trajectories With Initial Velocity, lOUcr
100 90
80
7 70

.-2 60
"1
c 50
0
._)
E 40
-5
0
z 30
2010
-200 0 200 400 600 800 1000
Nondimensional Length, L/d
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 (Lid) Avg. Height (Hid) Avg. Velocity (uiu*)
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




20
18
_>16
0
a) > 14
12
C10
J .-
E 6
*0 Z A

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 0
18 + Average Height
* Average Velocity O
16 16
o 0
U)
> 14
-
0)0
D12
-c 0
c10
_ 0
Ca c8
E6 0 *K
.E 6- 0o
a X
0 0 *K

5 6
Initial Velocity, n*ucr

Figure 4.29 Average Saltation Length And Height From Particles Ejected With Flow

I I I i
0 Average Length 0
+ Average Height
* Average Velocity 0
0
0
0
0
0 X
X*
O *K +
X + +
+ +
++
+
+
+
S I I I II I I I




Table 4.3 Saltation Characteristics With Grain Ejected With Shear Flow Initial Vel. (fl*ucr) Avg. Length (Lid) Avg. Height (Hid) Avg. Velocity (ulu*)
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 1 58.19 149.72
ejected at velocities ranging from ur to 1 OUcr 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 1 OOucr 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 revs 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 I
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 (Uuyf +W2

equation A. 1

From figure A. 1

U -f
cose
(uu2 +W2

equation A.2




sinO =
(u-u1)2 +W2

equation A.3

The drag force is defined as

equation A.4

Fd =CDA p Vr2
2

The longitudinal component is
Fd, = -CDA p|V|VcosO
2
=-CDA p (u-u)2 +w2 (UU) +W2 U 2 2
2 (uu) +W2

=-CDA 2 p (u-uY)2 +w2 (u-uf)

equation A.5

The vertical component is
1
F, = -CDA -pVIVsinO
2
=-CDA'p (u -u)2 +W2 (Uu12 +W2
2 (u-u1) +w2

equation A.6

=-CDA p (u u)2 +W2 ()




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

dii
t
Fb = const dt dt
0 -,It-c

equation A.7

Let a,, a2, a3, etc. represent the particle's acceleration over two time steps with a as the initial acceleration. For the first time step, At, the Basset history integral may be written
t A
f(aa 2 At-t
0rAAt-t )2 t
=-(a, +t a2 X.\A -'At -- O)

(a, + a2}v'-

equation A.8

For the next time step, the Basset history integral may be written




a + a2
At d
rf 2 =
2At-t

a2 a3
+ 2At-T

= _(a,+ ,At (a +2aA 2A t- t
2 2,2At- -,r' 2 2i 2-c At

=[(a +a2 XI2 1)+(a02 +a3)]At

equation A.9

Therefore in summation notation, the term can be expressed as
S T dtJr
o
a +2a (2 At) n + 1- j + j)+ a, +a.+, (2At)
X~ 2 n+-+~2

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/(Re^2) + 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 C1 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'




Cli=0.85*Cd; %%or remark this out if C1=0.2 end
gravx=rho*(s-1)*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(1)-fltop)^2 + vel(2)^2); ubot = sqrt((vel(1)-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 gray
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 gray
veldot=[dudt;dwdt]; end




"random.m"
clear global dudt dwdt tf btu btw s d g flow fltop flbot C1 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=[];

%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-1)/g); lscale=(s+Cm)*(10/3)*d;
ufnon=uf/vscale; %%Nondimensional shear velocity

Ks=2*d;
zO=0.11*(nu/uf) + 0.033*Ks;

((d/2)^2))/sqrt(pi*nu);




%%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=1; %%This loop is for firing grains into and with flow
for loop=1l: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==1
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
elseif loop==2
if m==l
u0=-sqrt((init^2)/(l+(tan(ang(m)))^2));
wO=abs(u0)*tan(ang(m));
else
direct=sum(ang(l:m));
u0=-sqrt((init^2)/(+(tan(direct))^2));
w0=abs(u0)*tan(direct);

end

velO = [uO w0 2*pi*20]'; % initial conditions uvel(1)=vel0(1);
z=0.5*d;
%%Call fluidvel.m to determine fluid velocity fluidvel
fluid (1) =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 C1=0.85*Cd;
%%Initial gravity force/buoyancy force gravx=rho*(s-1)*g*V*sin(beta); gravz=-rho*(s-1)*g*V*cos(beta);

u0 = ? w0




%%Initial drag force dragx=-Cd*Ar*0.5*rho*vrel*(vel0(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((ve10(1)-flbot)^2 + velO(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)*(velO(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(1)=(bassetz +gravz + dragz + liftz...
+magz)/(rho*s*V*admass);
%%Initial grain position dep(1)=0.5*d; x=0;
len(1)=0;
time(1)=0;
%Calculation of velocities until grain hits bed while jones==0
[t,vel] = ode45('motion',tO,tf,velO);
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 velO=[vel(length(vel),l) vel(length(vel),2)...
vel(length(vel),3)] '; t0=tf;tf=tf+deltat; A=[A;B];
uvel(i)=velO0(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=1;
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 C1 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=1; %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)/(1+(tan(ang))^2)); else
ang=180-ang;
ang=pi/180*ang;
xdot=-sqrt((init^2)/(l+(tan(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(1)=flbot;
%%Calculation of velocity gradient across grain dudzs(1)=dudz;
%%Determination of relative velocity vrel=sqrt((velO(1)-flow)^2+vel0(2)^2); Vr(1)=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(1)=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'
%Cl=0.2;
Cl=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(1)=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-l)*g*V*cos(beta);
%%Initial drag force if drag=='Y'
dragx=-Cd*Ar*0.5*rho*vrel*(vel0(1)-flow);
dragz=-Cd*Ar*0.5*rho*vrel*vel0(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(1)-fltop)^2 + ve10(2)^2); ubot=sqrt((vel0(1)-flbot)^2 + vel0(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*(vel0(3)-0.5*dudz)*(vel0(2)/vrel);
magz=-(pi/8)*(d^3)*rho*vrel*(vel0(3)-0.5*dudz)*((vel0(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(1)=x;
time(1)=0;
%Calculation of velocities until grain hits bed while jones==0
[t,vel] = ode45('motion',t0,tf,vel0);
B=[t vel];
%%Calculation of grain position
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;
%%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))*sgrt(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),l) vel(length(vel),2) ...
vel(length(vel),3)] '; else
velO=[vel(length(vel),l) 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=1;
for q=1: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) ,l); 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)