Computational Modeling of Advanced Mechanisms for Bridge Failure in Response to Fluid Forces

MISSING IMAGE

Material Information

Title:
Computational Modeling of Advanced Mechanisms for Bridge Failure in Response to Fluid Forces
Physical Description:
1 online resource (121 p.)
Language:
english
Creator:
Robeck, Corbin A
Publisher:
University of Florida
Place of Publication:
Gainesville, Fla.
Publication Date:

Thesis/Dissertation Information

Degree:
Master's ( M.S.)
Degree Grantor:
University of Florida
Degree Disciplines:
Civil Engineering, Civil and Coastal Engineering
Committee Chair:
THIEKE,ROBERT J
Committee Co-Chair:
INGLEY,HERBERT A,III

Subjects

Subjects / Keywords:
bridge-failure -- computational-mechanics -- fluid-structure-interaction -- scour -- wave-forces
Civil and Coastal Engineering -- Dissertations, Academic -- UF
Genre:
Civil Engineering thesis, M.S.
bibliography   ( marcgt )
theses   ( marcgt )
government publication (state, provincial, terriorial, dependent)   ( marcgt )
born-digital   ( sobekcm )
Electronic Thesis or Dissertation

Notes

Abstract:
Computational fluid dynamics was used to analyze the factors of three failure mechanisms for bridges. Two mechanisms involved the scouring of bed material around foundation and support structures, while the third involved hurricane wave forcing on the bridges superstructures. Analysis was primarily numerical and was meant as a complement to previous experimental and analytical work. Computational models appeared to reveal previously-overlooked factors associated with experiments. Additionally, computational wave models appeared to be capable of reproducing experimental wave signals and associated forces to a high level of accuracy.
General Note:
In the series University of Florida Digital Collections.
General Note:
Includes vita.
Bibliography:
Includes bibliographical references.
Source of Description:
Description based on online resource; title from PDF title page.
Source of Description:
This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Statement of Responsibility:
by Corbin A Robeck.
Thesis:
Thesis (M.S.)--University of Florida, 2014.
Local:
Adviser: THIEKE,ROBERT J.
Local:
Co-adviser: INGLEY,HERBERT A,III.

Record Information

Source Institution:
UFRGP
Rights Management:
Applicable rights reserved.
Classification:
lcc - LD1780 2014
System ID:
UFE0046723:00001


This item is only available as the following downloads:


Full Text

PAGE 1

1 COMPUT ATIONAL MODELING OF ADVANCED MECHANISMS FOR BRIDGE FAILURE IN RESPONSE TO FLUID FORCES By CORBIN ROBECK 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 2014

PAGE 2

2 2014 Corbin Robeck

PAGE 3

3 To those who pushed me to never give u p

PAGE 4

4 ACKNOWLEDGMENT S The author of this thesis would like to thank the researchers at the Argonne National Laboratory (ANL) Transportation Research and Analysis Computing Center (TRACC), particularly Steve Lottes and Cezary Bojanowski for t heir help and ongoing support Additional thanks to the member s of the Turner Fairbank Highway Research Center Hydraulics Laboratory, particularly Kornel Kerenyi, for his enco uragement throughout the studies My greatest level of appreciation is extended t o the members of my committee : Dr. Skip Ingley has been a con stant source of support and has always expressed a limitless amount of excitement for all of the work I have conduct ed Dr. Robert Thieke has been a perpetual guide throughout my academic career. His role has shifted over my time at the university from fluid mechanics professor, to academic advisor, coauthor, and most recently as the chair of my c ommittee. Finally to my advisor Dr. Raphael Crowley for his persistent pursuit of excellence and refusal to accept anything less than my best.

PAGE 5

5 TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................ ................................ ................................ .. 4 LIST OF TABLES ................................ ................................ ................................ ............ 7 LIST OF VARI A BLES ................................ ................................ ................................ .... 11 CHAPTER 1 IN TRODUCTION ................................ ................................ ................................ .... 1 5 2 COMPU TA T IONAL MODELS AND NUMERICAL ALGORITHMS .......................... 18 2.1 Governing Equations Formulation ................................ ................................ 18 2.2 Turbulence Models ................................ ................................ ........................ 21 2.2.1 Direct Numerical Simulation ................................ ................................ 25 2.2.2 Large Eddy Simulation ................................ ................................ ........ 27 2.2.3 k ................................ ................................ ........................... 31 2.2.4 k dels ................................ ................................ .......................... 36 2.2.5 Spalart Allmaras ................................ ................................ .................. 38 2.3 SIMPLE Algorithm ................................ ................................ ......................... 39 3 COMPUTATIONAL MODELING OF BED MATERIAL SHEAR STRESSES IN PISTON TYPE EROSION TESTING DEVICES ................................ ..................... 41 3.1 SERF Modeling Overview ................................ ................................ ............. 41 3.2 Background ................................ ................................ ................................ ... 42 3.2.1 Rotating Devices ................................ ................................ ................. 45 3.2.2 Submerged Jet Devices ................................ ................................ ...... 47 3.2.3 Piston Devices ................................ ................................ .................... 49 3.3 Approach ................................ ................................ ................................ ....... 53 3.4 Computationa l Model ................................ ................................ .................... 55 3.4.1 Mesh Geometry ................................ ................................ ................... 55 3.4.2 Inlet Formulation ................................ ................................ ................. 55 3.4. 3 Outlet Formulation ................................ ................................ ............... 56 3.4.4 Wall Treatment Formulation ................................ ................................ 57 3.4.5 Model Computation Scheme ................................ ............................... 59 3.4.6 Turbulence Model Formulation ................................ ........................... 60 3.4.7 Justification ................................ ................................ ......................... 60 3.5 Results ................................ ................................ ................................ .......... 61 3.5.1 Configurations ................................ ................................ ..................... 61 3.5.1.1 Conical protrusion ................................ ................................ ........ 67 3.5.1.2 One mm protrusion ................................ ................................ ...... 68 3.5.1.3 Differential erosion ................................ ................................ ....... 68 3.5.1.4 Wavy ................................ ................................ ............................ 68 3.5.2 Analysis ................................ ................................ ............................... 74

PAGE 6

6 4 AN OVERLOOKED LOCAL SEDIMENT SCOUR MECHANISM ............................ 78 4.1 Governing Equations for Flow Around a Cylinder ................................ ......... 78 4.2 Computational Modeling ................................ ................................ ............... 82 4.2.1 Model Dimensions ................................ ................................ ............... 83 4.2.2 Mesh Geometry ................................ ................................ ................... 84 4.2.3 Turbulence Model Formulation ................................ ........................... 87 4.2.4 Wall Treatment Formulation ................................ ................................ 87 4.2.5 Inlet and Outlet Formulation ................................ ................................ 87 4.2.6 Model Computational Scheme ................................ ............................ 88 4.3 Results ................................ ................................ ................................ .......... 88 5 WAVE IMPACTS ON BRIDGE DECKS ................................ ................................ .. 91 5.1 Background ................................ ................................ ................................ ... 91 5.2 Experimental Data Set ................................ ................................ .................. 93 5. 3 Modeling Overview ................................ ................................ ........................ 97 5.4 Computational Model ................................ ................................ .................... 98 5.4.1 Model Geometry ................................ ................................ .................. 98 5.4.2 Po lyhedral Meshing Scheme ................................ .............................. 99 5.4.3 Volume of Fluid (VOF) Model ................................ ............................ 100 5.4.4 Moving Boundary Method ................................ ................................ 102 5.4.5 Stokes Wave Theory Method ................................ ............................ 105 5.4.5.1 F irst order (linear) wave f ormulation ................................ .......... 106 5.4.5.2 Fifth o rder ................................ ................................ .................. 106 5.5 Anal ysis and Results ................................ ................................ ................... 107 5.6 Further Work To Be Conducted ................................ ................................ .. 114 6 SUMMARY AND CONCLUSIONS ................................ ................................ ........ 115 LI ST OF REFERENCES ................................ ................................ ............................. 116 BIOGRAPHICAL SKETCH. ................................ ................................ ......................... 121

PAGE 7

7 LIST OF TABLES Table page 3 1 Average amplification factors for each SERF configuration ................................ 69 5 1 Reported valu es for run BSXX014 ................................ ................................ ...... 94 5 2 Values for moving boundary method equation ................................ ................. 107 5 3 Values for linear wave theory model ................................ ................................ 107 5 4 Values for fifth order wave theory model ................................ .......................... 108 5 5 Wave forces modelled vs experimental ................................ .......................... 114

PAGE 8

8 LIST OF FIGURES Figure page 2 1 Visualization of density gradients from large eddy simulation of exhaust flow from a 2 inch conic nozzle conducted by researchers at Argonne National Laboratory ................................ ................................ ................................ .......... 30 3 1 Photograph of SERF ................................ ................................ .......................... 54 3 2 Schematic Layout of SERF showing three dimensional grid .............................. 54 3 3 Modeled shear stress versus measured shear stress ................................ ........ 62 3 4 Correlation between wall friction factor, f and grain size ................................ ..... 64 3 5 Computed shear stress using pre ssure drop from CFD model for varying roughnesses ................................ ................................ ................................ ....... 65 3 6 Grid sensitivity study results ................................ ................................ ............... 66 3 7 SERF sample section configurations ................................ ................................ .. 67 3 8 Amplification factors for conical protrusion configuration. Flow directi on is along the x axis. ................................ ................................ ................................ 70 3 9 Amplification factors for 1 mm protrusion configuration. Flow direction is along the x axis. ................................ ................................ ................................ 70 3 10 Amplification factors for differential erosion configuration. Flow direction is along the x axis. ................................ ................................ ................................ 71 3 11 Amplification factors for wavy configuration when top edge is held flush with flume bottom. Flow direction is along the x axis. ................................ ............... 71 3 12 Amplification factors for wavy configuration when average elevation is held flush with flume bottom. Flow direction is along the x ax is. ............................... 72 3 13 Local shear stress variations for uniformly rough samples ................................ 73 3 14 Centerline shear stress variations for uniformly rough samples ......................... 74 4 1 Depiction of flow around a cylinder with associated variable representation in polar coordinates ................................ ................................ ................................ 83 4 2 Top view of pile with associated refinement region ................................ ............ 85 4 3 Top view of pile with separation inducing notches ................................ .............. 86

PAGE 9

9 4 4 View of prism layers along bottom boundary ................................ ...................... 86 4 5 Side view of model with flow ................................ ................................ ............... 89 4 6 Pressure contour an d gradient vectors for pile radius 0.02 m at bed .................. 90 4 7 Pressure gradient magnitude contours for pile radius 0.02 m at bed .................. 90 5 1 Physical wave force experiments conducted at the University of Florida ........... 92 5 2 BSXX014 Experimental wave height signal ................................ ........................ 94 5 3 BSXX014 Experimental vertical force signal ................................ ...................... 95 5 4 Spectral analysis of wave signal ................................ ................................ ......... 96 5 5 Schematic of deck configuration from p hysical experiment ................................ 97 5 6 Model bridge deck before meshing ................................ ................................ ..... 98 5 7 A typical polyhedral cell used in the wave tank mesh ................................ ....... 100 5 8 D epiction of invalid and valid use of the VOF model ................................ ........ 101 5 9 Volume of fluid depiction on polyhedral mesh around deck .............................. 102 5 10 Depiction of vertex motion ................................ ................................ ................ 103 5 11 Comparison of wave generated by moving boundary method with experimental wave ................................ ................................ ............................ 109 5 12 One to one comparison of wave generated by moving boundary method with experimental wave ................................ ................................ ............................ 109 5 13 Comparison of vertical force generated by moving boundary method (green) with experimental wave ................................ ................................ .................... 110 5 14 Comparison of wave generated by linear wave theory method (red) with experimental wave ................................ ................................ ............................ 110 5 15 One to one comparison of wave generated by linear wave theory method with experimental wave ................................ ................................ .................... 111 5 16 Comparison of vertical force generated by linear wave theory method (red) with experimental wave ................................ ................................ .................... 111 5 17 Comparison of wave generated by fifth order wave theory method (red) with experimental wave ................................ ................................ ............................ 112

PAGE 10

10 5 18 One to one comparison of wave generated by fifth order wave theory method with experimental wave ................................ ................................ .................... 112 5 19 Comparison of vertical force generated by fifth order wave theory method with experimental wave ................................ ................................ .................... 113

PAGE 11

11 LIST OF VARI A BLES (Unless Otherwise Noted) c Speed of sound D Intersection point between logarithmic layer and viscous sublayer D h SERF hydraulic diameter E Friction constant f Friction factor G b Turbulent production due to buoyancy term G k Turbulent production term g Acceleration due to gravity g Two layer wall blending function I Turbulence intensity ratio k Turbulent kinetic energy k s Equivalent roughness p Change in SERF pressure L Distance between SERF pressure ports Mass flow rate r Roughness parameter Re Reynolds number Re+ Wall roughness number Re y Wall Reynolds num b er S Modulus of the mean rate of strain tensor S Mean rate of strain tensor

PAGE 12

1 2 T Temperature U Mean velocity in SERF u + Non dimensional wall parallel velocity u Friction velocity V Grid cell volume v Fluid velocity v g Grid velocity W Rotation tensor x Cell centroids y Distance from wall y + Non dimensional wall coordinate Variables (Greek) Turbulent dissipation rate Ambient source value in source terms of TKE eqn. Reconstruction gradient data Dilation coefficient Von Karmann constant Wall (bed) shear stress Dynamic viscosity of water Turbulent viscosity Turbulent viscosity Kinematic viscosity of fluid Density of fluid

PAGE 13

13 Turbulent Schmidt numbers Mathematics Gradient of a scalar quantity S: T Dot product of two tensors, S and T Dot product of two vectors, b and c Cross product of d and e

PAGE 14

14 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 COMPUTATIONAL MODELING OF ADVANCED MECHANISMS FOR BRIDGE FAILURE IN RESPONSE TO FLUID FORC ES By Corbin R obeck May 2014 Chair: Robert J. Thieke Major s : Civil Engineering Computational fluid dynamics was used to analyze the factors of three failure mechanisms for bridges. Two mechanisms involved the scouring of bed material around foundation and support structures while t he third involved hurricane wave forcing on the Analysis was primarily numerical and was meant as a complement to previous experimental and analytical work. Computational models appeared to reveal p reviously overlooked factors associated with experiments Additionally, computational wave models appeared to be capable of reproducing experimental wave signals and associated forces to a high level of accuracy

PAGE 15

15 CHAPTER 1 INTRODUCTION T he beginning of the 21 st century saw a great rise in the versatility and adaptability of computational mechanics code s Specifically, in recent years computational fluid dynamics codes have grown at a pace matching machine numerical power. Due to the nonlinear nature of the Navier Stokes equations finding analytical solutions to anything but highly specific classes of fluid mechanics situations has proven difficult. Furthermore added complexity arises when analyzing real world applications in which multiple factors ( complex structural boundaries, moving boundaries, and fluid boundary interaction) are present. Therefore t he ability to numerically model analytically unsolvable problems has been a major milestone for the fluid dynamics and engineering community as a whol e ( Anderson 1995 ) Computationally inclined scientists have seen this increase in computational affluence as a path to replace expensive wide scale physical testing. Experimentalists however remain skeptical of numerical scheme reliability and its inheren t disconnect from the physical nature of engineering design (Anders on 1995) This paper is an effort to bridge the gap between the physical and numerical community by providing numerical model validation in three distinct cases where fluid boundary interac tions are directly responsible for the structural failure Included are two numerical models of physical case stud ies and a numerically verified new application of potential flow theory. T he physical experiments model led were the Sediment Erosion Rate Flume (SERF) tests at the University of Florida and the w ave forcing experiments conducted by Sheppard and Marin (2009) in the University of Florida wave tank between 200 5 and 2009 The SERF is a physical flume device used in the meas urement and prediction of

PAGE 16

16 river bed material erosion rates (Crowley et al. 2012 ) In addition to producing a numerically accurate numerical model an emphasis was placed on classifying the fluid dynamic mechanisms of scour and initiating a constructive rhet oric for the strengths and limitations of this, and all, empirical scour measurement devices. The wave experiments modelled consisted of numerous wave tanks tests conducted to analyze the forcing associated with the impact of waves upon bridge superstruct ures as a result of structural damage sustained during hurricanes in the early 2000s. A specific goal of the experiment, and therefore of the numerical model, was to accurate ly quantify the distinct vertical and horizontal forces associated with the wave i mpacts. Perhaps most importantly however was to address bridge deck structural integrity and failure prediction based on the wave deck interaction. The new application of potential flow theory is an attempt to classify a previously over looked mechanism f or hydraulically induced bridge scour. The analysis looked to further of the understanding of scour depth and g rain size (D 50 ) dependence. Much previous research has found a dependency between equilibrium scour depth ( y s ) and ratio between sediment grain s ize ( D 50 ) and pier width ( b ). Examples include Ettma (1980), Froehlich (1988), Melville and Sutherland (1989), Garde et al. (1993), Sheppard et al. (1995), Sheppard (2004), Sheppard et al. (2004), and Sheppard and Miller (2006). Sheppard (2004) provided an analytical basis for the b/D 50 dependency on equilibrium scour depth based upon potential flow theory. However, the fundamental assumption behind the Sheppard (2004) analytical expressions was that potential flow theory, which is known to be applicable in higher elevations in a water column, is also applicable near the bed. Investigators in this analysis revisit the Sheppard (2004)

PAGE 17

17 equations and provide more supporting evidence for their applicability near the b ed through the us e of high end computational modeling. The common theme present in the studies presented for discussion was the effort taken to use numerical methods to determine fluid structure interaction effects on possible cause s for coastal structure vulnerability.

PAGE 18

18 CHAPTER 2 COMPU TA T IONAL MODEL S AND NUMERICAL ALGORITHMS 2.1 Governing Equations Formulation The set of non linear partial differential equations governing fluid flow is known as the Navier Stokes equations. When the fluid is taken as incompressible the three dimensional Navier Stokes equations are a set of four equations: three momentum statements (for each flow direction), and one of mass continuity. The momentum equations are presented in vector form, for brevity, as: i n which V is the velocity vector, t is time, is the fluid density, g is the gravity vector, and is the fluid dynamic viscosity. While assuming the fluid flow to be incompressible Additionally, a c ontinuity relationship may be established: i n which u is the velocity component in the x direction, v is the velocity component in the y direction and w is the velocity component in the z direction. If the flow is laminar these equations completely describe all components of the flow. To evaluate these partial differential equations numerically a continuous integral must be obtained that applies across the entire control volume continuum. This contin uous integral form of the governing equations is presented as follows:

PAGE 19

19 In which the fluid density, is the porosity, V is th e continuum volume, v is the velocity, v g is the grid velocity a is the face area vector, S u is the control surface, I is the identity matrix, T is the viscous stress tensor, and f is any additional body force placed on the fluid. Gravity would be an example of one such additional fluid body force. The next step in preparing the governing equation for numerical analysis is to convert the continuous equations to one of finite sums. This translation from a continuous integral to a finite set of equations allows physical problems to be broken into many parts, solved as a set of multiple smaller problems and pieced together. T he domain is discretized into many control volumes for which the equation s are solved and summed together The discretized momentum formulation for control volume i when summed across the face f of each volume is as follows: Correspondingly the dis cretized mass flow equation is: In which is the uncorrected face mass flow rate and is the mass flow correction. The uncorrected mass flow rate is the flow rate which is solved directly from the discretized equations However because the governing equation is a continuous integral approximate d by sums it is sometimes necessary to add a mass flow correction to ensure continuity is enforced through out the control volume.

PAGE 20

20 Although this control volume formation is discr ete there is no guarantee that the system will behave linearly or that there will be a unique solution to any given set of conditions. It is therefor e necessary to adopt a finite volume scheme that relates the discretized equations to a n implicit linear ized system of equations based on a representative grid. Within the finite volume approach this grid corresponds to a physical domain split into many finite control volumes known as cells. This finite number of cells is an indication of the number of degr ees of freedom with in the linear system. Although the equations have now been approximated as a linear system solving the m simultaneously as a coupled system can prove challenging Additionally, numerical instabilities can arise when solving large coupl ed systems leading to model convergence issues A segregate d approach is therefore taken to split the equations into transient, source, convective, and diffusive terms These separated terms are then solve d separate ly within each cell. The discrete form of the control volume equations when related to the finite volume method is: In which the fluid density, is the por osity is an arbitrary cell scalar value, V is the cell volume, v is the velocity a is the face area vector, G is the grid flux, and S is the volumetric source term.

PAGE 21

21 2.2 Turbulence Models Taylor and von Karman (1937) proposed a qualitative definition of turbulence as irregular motion that makes its appearance in fluids, gasses or liquid, when they flow past solid surfaces (Biswas and Eswaran 2002) A more modern quantitative, definition characterizes turbulence by the random fluctuation of flow parameters away from the mean in both space and time (Pope 200 2 ) These random fluctuations even when small, can grow in a manner described mathematically as ch aos in which seemingly negligible difference s in input parameters can produce diverging solutions in a few number of iterations. The chaoti c nature of turbulence can be seen analytically in the governing equations as instantaneous velocity and pressure term s which become difficult to quantify if the flow extends beyond the laminar range Even if the flow is considered incompressib le, homogen ous, and isothermal velocity term s appear in the form of a nonlinear divergence located in the convective portion of the Navier Stokes momentum equations as: For which v is the velocity v ecto r in each dimension This equation can be made more accessible if presented in Cartesian coordinates as:

PAGE 22

22 i n which u is the velocity compo nent in the x direction, v is the velocity component in the y direction, and w is the velocity component in the z direction. These velocity terms describe the transfer of kinetic energy that commonly dominates turbulent fluid flows. While turbulent, incomp ressible flows contain pressure terms that are linear, they often interact with velocity terms which are highly nonlinear. This nonlinear interaction makes both terms difficult to quantify. Moreover t hese nonlinear parameters frequently span many orders o f magnitude and occur across a wide range of time and length scales (Moin and Mahesh 1998). As such there is no analytical solution for solving general turbulent fluid flow problems. As with many phenomena prone to chaotic behavior mathematically complete models have been difficult to formulate. Therefore t he most common way of accounting for these flow fluxes is to use a method known as decomposition. The decomposition method is based o n the mathematica l concept that any variable which is a funct ion in both space and time can be broken into th e sum of two parts: the mean value and a fluctuating component. This is represented by: In which is an arbitrary variable, is the time independent component of and is t he fluctuations from the mean value as a function of both space and time. Specifically this method is implemented by decomposing velocity and pressure components as such:

PAGE 23

23 In which barred variables represent mean values and primes as the fluctuations away from the mean as a function of space and time. The difficulty of solving these equations arises when deriving an analytic method for predicting these fluctuating components. This inability to account for these term s leads to a lack of closure in the governing equations because implementation of this technique will result in six unknowns and only four governing equations It is therefore necessary to empirically approximate some of the terms. Usually, in the context of turbul ent flow, the viscous fluxes are approximated empirically. Although the complex details of turbulence is thought to be beyond current fluid mechanics models it is generally possible to captur e the mean flow quantities related to turbulence. It is therefore necessary to invoke a turbulence model that provide s some method of averaging erratic flow quantities in an effort to quantify the mean flow characteristics (Celik 1999). A common way of approximating this mean flow field is by a stochastic time averaging method known as the Reynolds Averaged Navier Stokes (RANS) method which utilizes the aforementioned decomposition technique Within the RANS formulation th e nonlinear velocity terms in the governing equations are replaced with a turbulence tenso r referred to as the Reynolds stress tensor. I n this ten sor scaling arguments are used to label specific parameters as n egligible and mean values are considered constant From there the remaining terms, Reynolds S tresses ominate the flow. The Reynolds Averaged Navier Stokes equations, presented in Einstein notation for clarity, then take the form:

PAGE 24

24 w here the only terms which require an additional model are the apparent n stresses defined as: RANS model formulations are typically presented with one or two equations for determining the Reynolds Stress. These equations often provide independent mean transport equations for two parameters : turbulent length scales and turbul ent kinetic energy. In the case of single equation models one of the above parameters is chosen and an attempt is then made to approximate the other from the modeled quantity. This formulation often leads modern RANS equations to be referred to simply as equation models. Within the Reynolds time averaging nomenclature there are two major assumptions that must be made: local turbulence isotropy and local equilibrium within the flow (Celik 1999). Because of these basic assumptions the RANS mod els have inherent limitations when describing anything more than the macroscopic effects of turbulence. In an effort to more accurately describe the physical phenomenon that takes place in turbulent fluid flow complex fluid models have been proposed tha t address the problem of calculating turbulence at vastly different scales. Modern theory suggests that turbulence is not isotopic but rather dependent on the scale at which it occurs (Pope 2002) Small scale turbulence tends to more closely exhibit close to isotropic behavior wh ile large scale fluctuations are much more erratic (Celik 1999). More advanced

PAGE 25

25 models therefore calculate turbulent eddies on multiple, vastly varied, scales within the continuum. Although a great deal of turbulence closure models have been formulated four have been chosen, due to their popularity and robustness, for discussion. Of the models chosen one is of the advanced eddy formulation and three from the Reynolds Averaging family. For ease of navigation they are presented in decreasing order of complexity as: Large Edd y Simulation (Advanced Hybrid) K Epsilon ( Two Equation RANS) K Omega (Two Equation RANS), and Spalart Allmaras (One Equation RAN S) Of the RANS models chosen all are of the linear eddy viscosity formulation and therefore adhere to the Boussinesq hypothesis of a linear relationship between the Reynolds Stress Tensor and the mean Strain Rate Tensor ( Pope 2002 ). 2.2.1 Direct Numerical Simul ation Although not a turbulence model a moment will be taken to discuss Direct Numerical Simulation (DNS) as a starting point for turbulence in computational fluid mechanics. The DNS method works on the assumption that the governing equations are correct i n their current form No turbulence model is used in the DNS approach and as such all velocity terms linear and nonlinear, must be calculated directly across all dimensional scales from the governing equations. Accordingly, there is no reliance on empiric ally based mean flow parameters such as t he Reynolds Stress tensor since the velocity fluxes are directly numerically calculated Within the DNS theory, turbulence is dictated by eddy formation and decay described by the governing differential equations T urbulence itself is then described as consisting entirely of different scale eddies. F ormations of large scale eddies tend to become unstable and release energy in the form of turbulence. This release of energy is

PAGE 26

26 classified by the creation of multiple, sm aller eddies. This process continues until the scale at which eddies can no longer be sustained is reached (Kolmogrov 1941 a ). At this minimum eddy scale turbulent energy is no longer converted into smaller eddies but rather released as heat. This release o f heat is known as the turbulent kinetic energy dissipation rate and is generally given the symbol This minimum turbulent scale, below which eddies can not exist is derived in a set of classical turbulence papers ce (Kolmogorov 1941 a ). Contained in this theory is which dictate this minimum scale beyond which eddy formation ceases to occur These scales are given as: i n which is the minimum length scale, is the kinematic viscosity, is the turbulent kinetic energy is the minimum velocity scale, and n is the m inimum time scale. Turbulence is the n shown to be contributed to by eddies that span widely varied time and le ngth scales. As such capturing these phenomena numerically requires grid size and time steps to be on the order of these minimum scales ( Pope 2002 ). Consequently analysis has shown that computation time for DN simulations grow cubically as Reynolds numbers increase s Due to the small time step and grid size constraint DN simulations are regarded as the most computationally intensive approach to solving a flow condition currently in implementation ( Lee et al. 2013)

PAGE 27

27 Although DN S simulations are generally considered to be the closest to physically correct there are still some assumptions which must be made with smaller ones (Pope 200 2) While this is considered to dominate the process experiments and simulations have shown the existence of upward energy flowing from smaller eddies to larger ones. This reversal of energy dissipation progression is known as turbulent backscatter Currently the resea rch community is divided on whether it can be captured by the governing equations alone (Mo i n 1998) Another assumption which must be made to implement a DNS scheme is that turbulence occurrences are entirely random. While this assumption generally holds t here are situations where stable, large scale turbulent structures can form ( Pope 2002). Due to the extremely high computation cost and existence of numerical instabilities DNS is not readily used by industry or significant portion s of the engineering comm unity. Its current application has remained limited to small, low Reynolds number flows of academic and scientific interest regarding the micromechanics of turbulence ( Anderso n 1995 ). Currently the largest DN simulation to date was conducted by Argonne Nat ional Laboratory and University of Texas at Austin researchers in November of 2013. The simulation spanned 242 billon degre es of freedom and r equired 260 million processor core hours. Even so the simulation was limited to a Reynolds number of ~ 5000 and the academic problem of channel flow between two infinite parallel plates (Lee e t al 2013 ) 2.2.2 Large Eddy Simulation The Large Eddy Simulation (LES) is a hybrid theory first presented by Joseph Smagorinsky ( 1963 ) Within the LES formulation an attempt is ma de to bridge the gap

PAGE 28

28 between highly expensive direct simulations, which calculate all scales of turbulent eddies, and Reynolds Averaged methods, which calculate only the mean The assumption within the LES formulation is that turbulent flow characteristics are primarily dominated by the physics of large scale eddies within the continuum. Large scale eddy formation and decay is therefore solved directly in a manner similar to full scale DNS. S maller scale eddies are considered less important to the overall f low and modeled using an averaged approach consistent with the RANS formulations ( Celik 1999). This approach provides more accurate results than those of RANS calculations however computation costs can still be considerable ( Cd a dapco 2012). To formula te t he LES equations the discretized Navier Stokes equations are put through a filtering process, (in lieu of averaging ) This filtering process determines the smallest eddy to be calculated directly, beyond which t he RANS approach will take over. Through alg ebraic manipulation the filter equations can be presented in a manner similar to the segregated flow equations with the exception of the Reynolds Stress tensor. Th at term is then replaced by the Boussinesq approximation tensor as: In which T is the turbulent stress tensor, t is the subgrid scale viscosity, S is the strain rate tensor, v is the vector velocity, is the fluid density, k is the turbulent kinetic energy, and I is the identity matrix. From there the strain rate tensor, S is resolved as: Because the LES formation makes no attempt to calculate small scale eddy structures a model is needed to close the governing equations and define turbulent viscosity beneath a specifi ed resolution at which eddy calculate is deemed

PAGE 29

29 unnecessary This is appropriately referred to as a subgrid scale model. Although multiple subgrid scale models have been put forward the most widely used is the Wall Adapting Lo cal Eddy Viscosity (WALE) mode l d ue to its relative ease of use Additionally, the majority of the LES computational effort has been placed in optimizing the WALE m ethod for high performance computing (CD adapco 2012) The WALE method is therefore presented here, in place of other subg rid methods, as representative of the LES formulation The WALE model presents a mixing length argument similar to that of Prandtl T urbulence T heory in which subgrid scale viscosity is defined as: w is a deformation parameter. The grid filter width is then defined as a delta function as: In which C w constant (Generally taken to be 0.41 ) and d is the wall distance. The deform ation parameter is defined as: Where S w is the deformation parameter, S is the strain rate tensor defined above, and S d is a modified strain tensor defined as: for which v is the velocity vecto r.

PAGE 30

30 The Large Eddy formulation is considered a significant step up from formulations which only consider mean flow parameters across the entire field ( Figure 2 1) Since large eddy configurations are calculated from the governing equations solutions more a ccurately reflect the physics of turbulent fluid flow. A consequence of calculating these additional eddy structures is a significant increase in computation time over two equation models. Although LES models do not approach the machine power required for a DN simulation studies have shown an order of magnitude increase in computation time over comparable two equation schemes. Another inherent limitation of LE simulations is their disregard of small scale eddy structures and the related assumption of isotr opic turbulence in those regions. Figure 2 1 Visualization of density gradients from large eddy simulation of exhaust flow from a 2 inch conic nozz le conducted by researchers at Argonne N ational L aboratory ( Joseph Insley ANL 2013)

PAGE 31

31 Despite the relatively high computational power needed to simulate flows of engineering interest there has been a large number of scientific publications and industrial reports for use of the method. As computational power continues on the path to exascale computing ( Argonne National Laboratory 201 3 ) use of LES and other highly computationally expensive methods is expected to increase as access to supercomputing capability become s more widespread. 2.2.3 k Models The k model is a two equation turbulence closure model in which empirical parameters for mean flow values are solve d Like other RANS model s the turbulent kinetic energy, k is used to define the turbulence tensor. The model differs from other two equations models in its use of a turbulent energy dissipation rat e, as the second parameter to close the governing equations. The k approximation 1 (1941a) theory describing turbulent kinetic and the rate at which the process progresses. However, the first formal theory was described by Jones and Launder (1972 ) as an empirical method for obtaining turbulent flow parameters to numerically evaluate the governing equations. This original formulation is c ommonly referred to as the standard k model. And while it was considered a large breakthrough at the time the model has inherent limitations ( beyond the limits of mean flow models ) and is generally not used in modern CFD analysis (Bardina et al. 1997). One such limitation with the standard model is its lack of realizability. Realizability places a constraint on the growth rate for turbulent kinetic energy with the grid. A shortcoming of initial formulation is the propagation of turbulent kinetic energy f rom stagnation points contaminating the surrounding solution. A second

PAGE 32

32 limitation was the inherent problem of resolving boundary layers at the edge of the domain. Standard formulations exhibit significant instabilities and error around wall boundaries. As a consequence many modified k models, which account for these limitations, are used currently. Popular variations of the original model currently in implementation include: Abe Kondoh Nagano (AKN), V2F, and the Realizable Two Layer models. The AKN and V 2F models are typically suited for low Reynolds number flow in which head exchange through complex geometry is of interest. The V2F model specifically is known to capture skin friction and near wall turbulence quite well. The Realizable Two Layer model is generally considered the most flexible model and it implemented in a wide range of applications. Most significant ly the Realizable Two Layer formulation allows for high wall distance values required when calculating shear stress due to surface roughness. A realizable, two layer, k model with two layer all y+ wall treatment was chosen for all studies presented in this thesis ( Crowley et al. 2014). This realizable model was developed by Shih et al. (1994), and it contains a different transport equation for the turbulent dissipation rate, than the traditional k approach. It also parameterizes the model coefficient, C as a function of mean flo w and turbulence. In the standard k approach, this term is assumed to be constant. This parameterization appears to be consistent with experimental observations in boundary layers (Cd adapco 2012).

PAGE 33

33 Specifically, in the realizable k model, the standard transport equations are: ( 2 1 ) ( 2 2 ) in which is the density of the fluid; k is the turbulent kinetic energy; V is the cell volume; v is the velocity; v g is the grid velocity; a is the face area vector; is the dynamic viscosity of the fluid; t is the t urbulent viscosity; k and are the turbulent Schmidt numbers; is the turbulent dissipation rate; 0 is the ambient turbulence value in the source terms that counteracts turbulence decay; M is the dilation dissipation coefficient; is the kinematic viscosity of the fluid; and S k and S are user specified source terms (CD adapco 2012). The turbulent terms, G k and G b for production and production due to buoyancy are given by Equation 2 4 and Equation 2 5 respectively: ( 2 3 ) ( 2 4 ) In these equations, T is the temperature (although flow was assumed to be isothermal); is the coefficient of thermal expansion; g is the acceleration due to gravity; and S is the modulus of the mean rate of strain tensor, S given by: ( 2 5 )

PAGE 34

34 ( 2 6 ) The dilation coefficient M is expected to be small since the SERF was modeled as a single phase volume of water. However, for completeness, its formulation is presented below: ( 2 7 ) where c is the speed of sound and C M = 2 (CD adapco 2012). The turbulent visco sity is defined as: (2 8) where, as mentioned, unlike the standard k approach, C is variable and is given by: (2 9 ) (2 10 ) ( 2 11 ) (2 1 2 ) (2 1 3 ) (2 1 4 ) W is the rotation tensor: (2 1 5 ) The model coefficient, C is defined as: (2 1 6 )

PAGE 35

35 where is defined as: (2 1 7 ) And the balance of the model coefficients are given as C 2 = 1.9, = 1.0, and = 1.2 (CD adapco 2012). This realizable model was combined with a two layer approach (Rodi 1991) which allows the k model to be applied in the viscous sublayer. In this approach (as implied), the computation is divided into t wo layers. In the layer near the walls, and t are specified as functions of wall distance. Values for are blended smoothly with values computed by solving Equation 2 2 far from the wall using blending functions described in detail by Jongen (1998). Meanwhile, Equation 2 1 is solved throughout the entire flow regime (CD adapco 2012). Specifically, the near wall model is parameterized as a length scale function and a turbulent viscosity ratio function (Wolftstein 1969): (2 18 ) ( 2 19 ) ( 2 20 ) ( 2 21 ) ( 2 22 ) where C = 0.09, = 0.42, A m = 70.

PAGE 36

36 The dissipation rate is simply: (2 23 ) The use of this two layer approach allows for the flexibility of a two layer all y+ wall treatment. 2.2.4 k Like the k model t he k model is a two equation turbulence model that belongs to the RANS family. Within the formulation two, independent, mean flow characteristics are empirically defined as the turbulent kinetic energy, k and the specific dissipation rate, The formation of the k in its original form was presented by Wilc ox in 1988. The model has had a somewhat tumultuous history to reach its current form. The original model presented significant iss ues with numerical stability. Additionally the model shows a high sensitiv ity to the empirical dissipation parameters in free stream flow. This poses an issue when estimating initial flow conditions. This parameter is often difficult to measure experimentally and the input of an expected value often caused numerical divergence One advantage of the k model that was deemed important was its robustness in boundary layer s close to walls. In an effort to capture this boundary layer advantage Menter (1996) rearranged the formulation in 1996 to create a two equation hybrid model. This new model evoked the k equations in areas of close relation to the wall and called upon k assembly. This hybrid RANS model was labeled as the Shear Stress Transport (SST) k Omega model. Wilcox (1998) released a set of corrections to his original k model; although the theory is thought to be sound there has yet to be wide scale use of the corrected model. Consequently model validation has not been

PAGE 37

37 established and care must be taken when implementing the corrected expressions. Due to its wide use a nd large presence of validation studies the SST variation of the model will be presented here as representative of modern k expressions. The SST transport equations are formulated as: I n which is the density of the fluid; k is the turbulent kinetic energy; V is the cell volume; v is the velocity; v g is the grid velocity; a is the face area vector; is the dynamic viscosity of the fluid; t is the turbulent viscosity eff is the effective intermittency is a compressibility correction constant, f is a free stream dissipation limiter coefficient, S k is a user defined source term, and D is the cross derivate term The SST k model is currently a popular two equation model within the aerospace and wind engineering industry due to its ability for rel atively, accurately descriptions of the mean characteristics of compressible flow and boundary layer separation. The model is generally used in place of more complex models due to its computational efficiency and avoidance of numerical instabilities associ ate d with detailed turbulent flux calculations.

PAGE 38

38 2.2.5 Spalart Allmaras The Spalart Allmaras is a one equation model, meaning that it aims to resolve the entire relationship between the Reynolds stress and mean flow velocity field with a single transport expressi on. This expression comes in the form of a model for the parameter known as the eddy viscosity. Previous one equation models chose to combine an equation for turbulent kinetic energy and an algebraic description of length scale to define eddy viscosity whe reas the Spalart Allmaras seeks to solve directly for the eddy viscosity. Although this method is a step more complex than algebraic (No Equation) models critiques often note its lack of a turbulent length scale parameter (Wilcox 1998) The model was origi nally developed in 199 2 for use in the aerospace industry by NASA engineers (Spalart and Allmaras 1992) specific formulation with the direct intention of computational implementation. The equations are presented in a fo rm that is intentionally written in a manner to allow for straightforward translations to unstructured grids within computational codes. Although unstructured grids were fairly new at the time of publication they have become increasingly prevalent in recen t years. This attention to bridging the gap between the theoretical model and the struggles of implementation has led the model to become one of the more popular approaches within the commercial CFD and aerospace industry (Moin 1998). Another si gnificant a dvantage of the Spal art Allmaras approach is the reduced computational cost when compared to other widely used models. Since the model only makes use of a single additional transport equation it can yield design caliber results at significantly lower compu tation cost than higher complexity models.

PAGE 39

39 Much like many turbulence models there exists a standard form which adheres to the original publication and various modified forms which have been developed over the years to overcome hurdles the formulation has encountered as use became more mainstream. The standard Spalart Allmaras is generally only suitable for flows of low Reynolds number flow which are free of shocks and have little if any, separation ( Spalart and Allmaras 199 2 ). A disadvantage of the origin al formulation is that it is implement ed wit hout the use of wall functions. Due to this omission if the model is to be used to acquire even reasonable accuracy it requires a fairly fine grid at boundaries to capture wall turbulence in the boundary layer (C d adapco 2012) It is also noted that although the model provides reasonable results for wake, mixing layer, and radial jet problems it does not provide accurate results for plane and round jet problems (Wilcox 1998) primarily due to its lack of an equatio n to describe the eddy length scales involved with the problem formulation. The model equations are generally considered less complex than those of the other RANS models presented. Therefore the presentation of empirical constants and derivation of model v ariables is left to the reader. 2.3 SIMPLE Algorithm Once appropriate models have been selected the segregated fluid approach discussed above was implemented using the Semi Implicit Method for Pressure Linked Equations (SIMPLE) method (Patankar and Spalding 1972) internally within a variety of CFD solvers (Andersen 1995) In a general sense the SIMPLE algorithm is used to set up the workflow when solving impli cit unsteady transient problems. The formulation is broken into twelve steps; the proce ss is then repeated for each time iteration. For completeness the workflow is presented below as:

PAGE 40

40 1. Set the boundary conditions 2. Compute the reconstruction gradients of velocity and pressure 3. Compute the velocity and pressure gradients 4. Solve the discretized mo mentum equation 5. Compute the uncorrected mass fluxes, at cell faces 6. Solve the pressure correction equation 7. Update the pressure field by : 8. Update the boundary pressure corrections 9. Correct the face mass fluxes by: 10. Correct the cell velocities by: In which v is the uncorrected cell velocity, v* is the uncorrected cell velocity, is the vector of the central coeffic ient for the discretized linear system representing the velocity equation, and V is the cell volume 11. Update density due to pressure changes (If Flow is Compressible) 12. Free all temporary storage

PAGE 41

41 CHAPTER 3 COMPUTATIONAL MODELING OF BED MATERIAL SHEAR STRESSE S IN PISTON TYPE EROSION TESTING DEVICES 3.1 SERF Modeling Overview The Sediment Erosion Rate Flume (SERF) device was computationally modeled using CD CCM+ at varying flow rates and sample roughnesses so that wall shear stresses could be evaluated during a piston style erosion test. Shear stress data were matched between the model and data from previous physical tests. Pressure differential upstream and downstream from an eroding specimen displayed similar behavior during both p hysical and modeled tests in that as eroding sample roughness increased, pressure differential did not appear to change. A series of complicated bed configurations were added to the computational model to quantify sample over advancement during an erosion test. Results appeared to indicate that small deviations in sample geometry may have large effects on localized shear stresses. Another series of models was run to shear stress development over a rough sample. Results showed that relatively large localized shear stresses may develop simply because of the introduction of roughness. Finally, results also indicated that the most co nservative method for future testing is to keep the bottom edge of an eroding sample flush with a flume bottom.

PAGE 42

42 3.2 Background From the outset, the major point of emphasis within the numerical adaptation was placed on providing a clearer picture of the hydraulically induce d cohesive scour that occurs in empirical scour testing devices Although laboratory erosion rate measurements have been taken in a wide variety of conditions their relation to full scale physical scour amounts has yet to be fully e stablished (Crowley et al. 2012). Erosion rates within these devices have not been analytically described. As with any entirely empirically based system care must be taken when evaluating their results. Many of these scour devices rely on accurate shear s tress measurements and the governing assumption that shear stress is constant across the entire sample. This constant shear stress assumption causes the most pressing concerns. Localized shear stresses across a sample can cause laboratory tests to either s ignificantly over and under estimate material scour resistance. Underestimating the strength of a material could lead to the overdesign of coastal structure adding unnecessary expensive. Overestimation of material strength is of the highest concern as it w ould lead to the direct lowering of design safety factors, which could ultimately lead to structural failures. Although it can take several forms and s cour is broadly defined by t he Hydraulic Engineering Circular (HEC) as erosion of streambed or bank material due to flowing water (Crowley et al. 2013). Each form of scour has its own mechanisms however it is generally grouped into to three major categories: General, Contraction, and Local. General scour is defined as bed elevation changes that result from the later al instability of the waterway (Florida Bridge Scour Manual 2005). This change in elevation occurs due to the cyclic nature of tidal and riverine cycles. Occasionally river flows can reverse direction cau sing slow changing oscillations in the river veloc ity profile. These

PAGE 43

43 oscillations over time tend to shift bed material and rearrange river profiles. This change in global flow patterns occurs gradually and rarely causes rapid loss of bridge foundation mate rial. As a consequence of the large time scales at which global scour occurs s cientific work on the mechanisms has been limited. Because the current design life of coastal structures is generally limited to fractions of centuries global scour is not commo nly accounted for in the design process. Contraction scour occurs do the rapid decrease in channel cross sectional area due to an obstruction, natural or man made (Florida Bridge Scour Manual 2005). This effect can be caused by a variety of factors with th e most common being the construction of sizable bridge pier foundations. As the river cross section narrows the cross sectional velocity must increase ( assuming a constant flow rate ) This increase in mean flow velocity has a twofold effect. First, as the law indicates that a higher percentile of sediment sizes will be transported downstream with the moving flow. Secondly, over t ime the amount of sediment removed from the contracted area will increasing in a compounding ma nner This increase can cause an unbalance of sediment entering the contracted area versus leaving. Contraction scour is g enerally considered predictable. As a result foundations are designed to extend past t he maximum calculated scour point (Trammell 2004 ) Local scour can be considered similar to contraction scour but i n a much more specific, localized location. Because of the stagnation pressure (and associated local accelerations) present at the leading edge of a bluff body, high fluid forces develop at the base of coastal foundations. This increase in forces coupled with the increased contraction velocity is the main mechanism for local scour. Local scour is generally

PAGE 44

44 determined around specific supports rather than the entire structure all at once. The integrity of the structure is most highly dependent on local scour. Local scour generally recognized as the major factor in removal of material around open channel bridge piers. Due to its highly specific nature local scour can be the most difficult to pre dict and therefore design around. Consequently coastal design manuals have placed a large emphasis on local scour mitigation and have deemed it the largest contributor to modern bridge collapse in the United States ( Arneson et al. 2012 ). The effects of loc al scour were the primary focus of this study. Currently the physics of local bridge pier scour is not well understood and as a consequence mitigation is often overly conservative (Annandale 2004 ). This approach tends to compound construction costs and le ave designers with a wide possible range regarding the ultimate capacity of the structure. Bridge designs have reported accounts of d oubling the initial design pile depth embedment due to scour concerns ( Trammell 200 4). In response to the growing concern o ver bridge integrity the Federal Highway Administration (FHWA) released the first Hydraulic Engineering Circular No. 18 Evaluating Scour at Bridges in 1993. The current, fifth edition, was published in April of 2012 (Arneson et al. 2012) Within the manual, bridge pier scour is placed into two broad categories regarding the foundation material type: sandy cohesionless materials and clay like materials exhibiting strong cohesive properties. For cohesionless materials empirical formulas are presented based on laboratory tests from various institutions. One set of such equations was developed by Max Sheppard at the University of Florida

PAGE 45

45 ( Florida Bridge Scour Manual 2005). In the case of materials exhibit strong cohesion, however, such as clay and rock a physical e rosion test must be performed. The latest standard stipulates that when cohesive bed material erodes, a method similar to the Briaud et al. Erosion Function Apparatus Scour Rate in Cohesive Soils (EFA SRICOS) method (Briaud et al. 1999, 2001, 2004a, 2004b) may be used to estimate scour depth. The crux of this scour erosion function, or erosion rate versus shear stress relationship, which should exist for any bed material (Einstein and Kro ne 1962, Partheniades 1965) must be measured using laboratory equipment. The advantage to this semi empirical method is that site specific and soil specific conditions are taken into account during design. But its disadvantage is its reliance on an accur ate erosion test. Three erosion function measuring devices are discussed in HEC 18: rotating type devices, submerged jet type devices, and piston type devices 3.2.1 Rotating Devices Rotating scour devices began produc tion in 2001 with the development of the Ro tating Erosion Test Apparatus (RETA) at the University of Florida ( Crowley et al. 2012 ). The RETA testing scheme is only applicable to samples exhibiting str ong cohesive properties such as stiff clays, sandstone, and limestone. The reason for this is twofo ld. First, the sample is place d on a support shaft running through the center of the rotating device. Samples must therefore be partially supported vertically by their own cohesive forces. Secondly, the shear stress calculation is based upon a decrease in sample radius. For samples exhibiting low cohesion this reduction in sample radius represents a nonphysical property.

PAGE 46

46 To conduct a RETA erosion test field samples are collected and formed to a 4 aced through the sample material and positioned in the RETA device. An outer Plexiglas annulus is then placed around the sample material and filled with water. The outer cylinder is then made to rotate by a motor driven chain. The sample rod is then attach ed to a torque cell. As the outer cylinder rotates shear stresses accumulate on the sample causing material to erode. At the end of a RETA test samples are removed from the device and dried in an oven to remove saturated water. Erosion rates are then deter mined from: change in sample mass, change in sample diameter, and change in sample shape ( Crowley et al. 2012b ). The average shear stress, is then give by: In which T is the torque measure d at the torque cell, R is the radius of the sample, and L is the length of sample From t his the erosion rate can be calculated by: is the change in mass of the sample, R is the radius of the sample, L is the length of sample, and D is the duration of the test ( Crowley et al. 2012b ) Because t he torque on the inner cylinder is used to calculate the shear stress on the sample the RETA loop within the system mandating a constant shear stress. There has however been considerable criticism on the RETA device. Because the shear stresses in a river bed occur at the soil fluid interface ero sion will occur in a top to bottom fashion. Within the

PAGE 47

47 RETA erosion occurs along a radial plane which could cause inconsistencies when relating measured shear stress within the RETA to river bed shear st resses. Additionally sample shear stresses are consid ered to be constant across the entire sample. This constant shear stress assumption holds only if the cylinder is uniform and free of defects. This assumption causes issues when relating the shear stress of an uneven sample. Although care is taken to ensur e that the sample is sufficiently uniform before test ing some degree of error is expected during preparation. Therefore although the stresses on the rotated outer cylinder are approximately constant those on the sample can be quite localized ( Crowley et al 2012b) 3.2.2 Submerged Jet Devices The jet device operates on the concept that the shear stress in an eroding sample is related the maximum scour depth. An assumption is then made that as a scour hole grows the shear stress will decrease until equilibrium is achieved after which scour will no longer occur. This concept is shown through jetting experiments as a result of increasing jet dissipation energy (Weider 2012). The J et Erosion Test (JET) apparatus was developed by Hanson in 1997. The device consists of a jet tube, an adjustable head tank, and a jet submergence tank (Weidner 2012). The JET device is filled with water and placed so that its orifice is submerged in soil. Water is pushed through the JET orifice and the sample is made to erode. As the sample erodes depth of erosion as a function of time is recorded. From the erosion test, center line velocity of the jet can be found by:

PAGE 48

48 In which U is the center line velocity, U 0 is the jet nozzle velocity, Cd is the diffusion coeffic ient (Generally taken as 6.2), d is the nozzle diameter, and H is the distance from the jet origin. From there an assumption is made that the jet flow is parallel to the planar boundary and that the velocity parallel to bed is equal to the jet velocity at distance, H, from the jet nozzle. From this assumption it is t hen possible to calculate the shear stresses caused by the jet as follows: In which C f is the coefficient of friction. According to Hanson et a l (1997) river flows are dominated mostly by surface roughness Therefore, the coeffici ent of friction should be independent of Reynolds number give n a sufficiently fast moving flow. From this shear stress erosion is then measured until a maximum depth of erosion has occurred and equilibrium is said to be established. The maximum scour depth is then related to the critical shear stress by: The critical shear stress is then integrated through dimensional analysis and a theoretical time of erosion can be calculated for comparison with the experimental time. An upside of the submerged jetting device is its ease of use in the field Additionally, typical results are obtain ed in a relatively short amount of time; most tests typically run for around an hour (Hanson 2004). On the other hand a JET test require s a large number of assumptions which may affect its accuracy. For example, the test itself can cause hole formation of differen t shapes. Because scour holes are known to come in a variety of shapes (Crowley et al. 2012) this may cause issues when scour

PAGE 49

49 occurs unevenly. Additional ly when testing cohesive soils jet tests can take a s long as 80 hours. To mitigate large run time s, results are typically extrapolated outward from shorter time experiments. This extrapolation has been shown to overestimate scou r hole depth s found in full time length experiment s A theory for why this overestimation occurs is documented by Crowley et al. (2012). Crowley et al. proposed that soil density and stiffness increase the further down it is located. This increase in stiff ness can lead to reduced erosion rates and overly conservative scour estimates (Crowley et al. 2014) This extrapolation error compounds upon already exist ing concerns that studies have yet to show that full time length erosions tests accurately predict co hesive scour (Weider 2012). 3.2.3 Piston Devices This study focused entirely on piston type devices, which have been in use since the mid McNeil et al. 1996), the Adjustable Shear Stress Er osion and Transport Flume (ASSET, Roberts et al. 1998), the Erosion Function Apparatus (EFA, Briaud et al. 1999, 2001, 2004a, 2004b), and the Sediment Erosion Rate Flume (Crowley et al. 2012b, i.e. the SERF). To use these piston style devices, first an i n situ Shelby tube or rock core is obtained. The specimen is extracted into a testing tube, a piston is inserted into one end of the testing tube, and the assemblage is attached to a lead screw. This assembly is inserted into a rectangular duct with a ci rcular cutout so that piston advancement causes the sample to protrude into the duct. Water is made to flow over the sample,

PAGE 50

50 bottom. Thus, erosion is approximately c aused only by shear stress. Erosion rate, or y/ t ( y is change in piston position; t is elapsed time) is directly measured. However, it is often difficult to obtain accurate shear stress estimates in these instruments because no known method exists t hat is capable of measuring erosion rate and shear stress simultaneously for a Shelby tube specimen. Earlier devices such as stress by using a smooth wall approximation given by: ( 3 1) ( 3 2) (3 3 ) where f is the friction factor defined by the Darcy Weisbach Equation given in Eq. ( 3 3); U the average velocity; D h the hydraulic diameter defined by Eq. ( 3 2) ; the kinematic viscosity of water; h the height of the duct; w the width of the duct; and the shear stress on the four walls of the duct. The EFA SRICOS method recommends using the Colebrook White Equation (Eq. 3 4) with a roughness height, k s of one half the eroding sediment diameter to estimate f (Briaud et al. 1999 2001, 2004a, 2004b): ( 3 4)

PAGE 51

51 Annandale (2006) criticized this approach, and he recommended using a pressure drop measurement (White 1986) to estimate shear stress as per Eq. 5: ( 3 5) Where, p is the pressure drop across length, L ; and w and h height respectively. Crowley et al. (2012 sensor to directly measure shear stress at different roughnesses (smooth, 0.125 mm, 0.25 mm, 0.5 mm, 1.0 mm, and 2.0 mm). The sensor may be described as a spring type apparatus. It is inserted into a hole in the flume bottom, and as water flows over a 50 mm test disc, the sensor slightly deflects. Solenoids are used to close the system so that deflection is minimal. This study concluded that using a pressure drop to estimate shear stress across an e roding sample did not produce accurate results because the small. Results showed that as roughness on a relatively small portion of the duct increased, p remained n early constant for a given velocity (shear stress). Investigators also found that using k s equal to half the sediment diameter, as recommended by Briaud et al. (1999, 2001, 2004a, 2004b), appeared to be incorrect because back solving Eq. 3 4 for k s relies upon instrument specific hydraulic diameters. Crowley et al. (2012) also discuss forcing on a sample that erodes non uniformly. Annandale (2006) identified a situation in the EFA where the downstream portion of a specimen eroded much more quickly than its upstream edge. Similarly, during SERF testing, differential

PAGE 52

52 whereby erosion of smaller particles is interspersed by events where large chunks of material erode (Crowley et al. 2012b). Under these conditions, using uniform surface shear stress assumptions would appear to be incorrect because localized shear stresses should govern the erosion rate. In concept, the distribution of shear stress over a flat, rough sample should show some similarity to development of the boundary layer over a flat plate. At the leading edge of the flat plate, both boundary layer thickness and shear stress begin at zero and grow with distance downstream. In the case of a rough sample in a piston style devi ce the shear on a smooth duct boundary. And, presumably the shear over the remaining portion of the sample should increase as the flow moves downstream (at a faster r ate than shear stress growth for a smooth sample). This increase in shear development However, Crowley et al. (2012b) appear to show that from a design perspective, assuming nearly uniform surface conditions in the testing device may be conservative because estimating a lower shear stress for a given erosion rate forces the design erosion function to translate toward the y axis and to become steeper. Therefore, the most conse rvative design approximation would be to use smooth wall approximations do not significantly reduce the localized shear stress. Investigators sought to determine: (1) i f computational results verify experimental significant effect on shear stresses. Question (2) bears some further discussion beyond

PAGE 53

53 non uniform erosion. The Briaud e t al. (1999, 2001, 2004a, 2004b) EFA testing procedure stipulates that during testing, samples should protrude 1.0 mm into their duct. Presumably, this is done to ensure conservative results. Causing the sample to protrude into the flume will induce a no rmal stress upon its leading edge, which in turn should increase erosion rate for a given estimated shear stress. On the other hand, Crowley et al. (2012b), Roberts et al. (1998) and McNeil et al. (1996) recommend keeping the sample approximately level wi th the flume bottom so that field shear stress alone is approximated. 3.3 Approach To evaluate the questions presented in the previous section a computational model of the SERF (Crowley et al. 2012b) was generated using CD a CCM+. Details of the model are discussed below: Dimensions of the SERF are 2.0 in. (5.08 cm) by 8.0 in. (20.32 cm) by 9.0 ft. (2.74 m). The 1.0 ft. (30.48 cm) duct inlet and 1.0 ft. (30.48 cm) duct outlet are n. (10.16 cm) polyvinylchloride (PVC) connective plumbing. The inlet of the device is equipped with a flow straightener that is designed to aid in transition to a hydraulically smooth, fully turbulent flow. A photograph of the device is presented in Figu re 3 1 while a schematic diagram (and associated mesh approximation) is presented in Figure 3 2.

PAGE 54

54 Figure 3 1 Photograph of SERF (Crowley et al 2014) Figure 3 2 Schematic Layout of SERF showing three dimensional grid; (a) shows a three dimensional isometric view; (b) shows a top view. The location of the shear sensor is highlighted in the figure; the dark circles show the location of the sample section. (Crowley et al 2014)

PAGE 55

55 3.4 Computational Model 3.4.1 Mesh Geometry A detailed, three dimensional drawing of the SERF was prepared using AutoCAD by Autodesk, Inc. All volume meshing was performed within Star CCM+. First investigators attempted to adapt a polyhedral volume meshing scheme to the geometry. However, due to co nvergence issues with the model, researchers eventually were forced to abandon the polyhedral mesh and revert to a tetrahedral scheme. Between 1.5 and 2 million cells were used depending on the simulation. Cell resolution was approximately 0.5 mm per cel l near the shear stress sensor and sample sections. Resolution was expanded to approximately 10.0 mm further away from these areas of interest in improve computational efficiency. 3.4.2 Inlet Formulation were specified explicitly in the normal direction. The boundary face pressure was extrapolated from cells adjacent to the inlet boundary using a hybrid Gauss least square method (LSQ) reconstruction gradient. The unlimited (superscript u ) reconstruction (subscript r ) of data value is given by: ( 3 6 ) (3 7 ) (3 8 ) (3 9 ) (3 10 )

PAGE 56

56 where x 0 and x n are the centroids of cell 0 and its neighboring cell n through face f ; 0 and n are data values in cell 0 and its neighbor; A f is the face area vector; V 0 and V f are cell volumes; and is the geometric Gauss LSQ gradient blending factor field function. The blending factor is user specified, and it determines the weight given to the Gauss/LSQ computed gradients (CD adapco 2012). Initial turbulent kinetic energy, k was specified throug h introduction of a turbulence intensity parameter, I such that: (3 11 ) where v is local velocity magnitude and I was manually set at 5%. Inlet turbulent dissipation rate was derived using an initial turbulent viscosity ratio, of 30 such that: ( 3 12 ) 3.4.3 Outlet Formulation CCM+ defines the velocity at a flow split outlet as: ( 3 13 ) Where is the velocity that is extrapolated from the adjacent cell value using reconstruction gradients; a/|a \ is the outward normal vector; and x i is a scale factor that is computed for outlet boundary i as per Eq. 48 through Eq. 50 (CD adapco 2012): ( 3 14 ) where f i is the specified fraction of the flow leaving outlet i (in the case of the SERF, 100%); and is the total inlet flow defined as:

PAGE 57

57 (3 15 ) And is the unscaled outlet mass flow rate through outlet i : (3 16 ) where is the maximum density of the fluid. 3.4.4 Wall Treatment Formulation A two layer all y+ wall treatment implies that no assumptions are made about how well the viscous sublayer is resolved. Instead, a blended wall law was used to estimate shear stress. Blend ing is achieved between high y+ wall treatment and low y+ wall treatment. High y+ wall treatment is similar to a wall function type approach in that near wall cells are assumed to lie within the logarithmic region of the boundary layer. Low y+ wall treat ment is similar to the traditional low Reynolds number approach where no modeling beyond the assumption of laminar flow is needed in the wall cells. In other words, it assumes that the viscous sublayer is properly resolved (Cd adapco 2012). All y+ wall t reatment attempts to mimic high y+ wall treatment when a mesh is coarse and low y+ wall treatment when a mesh is fine. Specifically, a blending function, g is defined as a function of Re y : (3 17 ) Next, a reference (friction) velocity is d efined: ( 3 18 ) in which y is the normal distance from the wall.

PAGE 58

58 Wall cell turbulence production is defined as: ( 3 19 ) where u + is the wall parallel velocity non dimensionalized with u and y + is the non dimensional wall coordinate given by: ( 3 20 ) Wall cell dissipation is given the same formulation as Equation 2 23 to provide two layer approach consistency. At the walls, a Neumann boundary condition is used for k such that In the viscous sublayer, velocity distribution is modeled as: (3 21 ) In the logarithmic layer, velocity distribution is modeled as: (3 22 ) where E is a constant (9.0) and f is the roughness function (friction factor). Obviously, a discontinuity exists between the viscous sublayer and the logarithmic layer i.e. the buffer region. A blended wall law (Reichardt 1951) is used to model the velocity distribution in this region : (3 23 ) where D is the intersection point between the logarithmic layer and the viscous sublayer, and: (3 24 )

PAGE 59

59 (3 25 ) Wall friction factor is defined through introduction of a wall roughness number: ( 3 26 ) where r is a roughness parameter (analogous to k s ); and is the kinematic viscosity of water. Friction factor f is related to Re+ via: ( 3 27 ) where the exponent, a is given by: (3 28 ) where C and B are calibration coefficients and Re + smooth and Re + rough are constant conditions. By default, B = 0; C = 0.253; Re + smooth = 2.25; and Re + rough = 90 (Cd adapco 2012). 3.4.5 Model Computation Scheme Initial conditions of the model were e stablished such that the SERF was filled v 0 was specified at the velocity downstream. Average wall shear stress across t section and sensor section was used as an approximate indicator of a fully developed flow condition. Once these average shear stresses approached a constant value, or leveled off, the flow was

PAGE 60

60 assumed to be fully developed. Then, a mode conditions were typically achieved in approximately 4 5 seconds of modeled time. This time approximately corresponded to observations during physical tests in previous studies. Modeled readings were compared with data from Crowley et al. 2012a. Typically, average y+ values near the sample section and sensor section were between 75 and 80, which would tend to shift all y+ wall treatment toward the high y+, wall function style of wall boundary. This technique was deemed acceptable because physical data. 3.4.6 Turbulence Model Formulation As stated above a realizable, two layer, k model with two layer all y+ wall treatment was chosen for this simulation. Section 2.2.3 provides an in depth discussion of the model and other relevant information. 3.4.7 Justification The models that were used in this study were chosen for a number of reasons. CCM+ appea red to indicate that the combination of models described here would provide accurate results. Investigators (ANL) Transportation Research and Computing Center (TRACC) and the J. Sterling Jones Hydraulics Laboratory at the Turner Fairbank Highway Research Center experiences, they agreed that the aforementioned modeling choices appeared to be appro priate. Review of the Star similar to the one described here, the models used in this paper were appropriate for this particular commercially available software (CD adapco 2012).

PAGE 61

61 Finally, models dis cussed here have been compared with other modeling options in literature. For example, Davis et al. (2012) compared Star k model with its k model and its V2F k variant for flows over a wall mounted cylinder. The traditional k a pproach was not studied because of its known deficiencies in resolving flows in the viscous sublayer. This study concluded that the realizable k model consistently performed better than the other two models in almost every measurable aspect. 3.5 Results 3.5.1 C onfigurations Several simulations were conducted using this setup at varying flow rates. First, computation, researchers wanted to be sure that these wall functions accur ately replicated real data for the selected mesh geometry and associated y + values. Once smooth data had been verified, the roughness parameter, r from Eq. 326 was varied to match data for grain sizes of 0.125 mm, 0.25 mm, and 0.5 mm at each velocity (as reported in Crowley et al. 2012a). Averaged modeled shear stress across the sensor section was plotted as a function of measured shear stress from Crowley at el. (2012a) to demonstrate quality of fit (Figure 3 3 ). Note that the same roughness parameter was used at each velocity for a given grain size to preserve physical significance. For rougher sediments, increases in r led to increases in y + This, in turn led to increases in the near wall cell centroids near the areas of interest in the flume (reca ll that resolution r value was defined as approximately 0.3 mm the approximate distance from the wall

PAGE 62

62 B from Eq. 3 27 was adjuste d to fit the measured data ( C in Eq. 38 was left untouched so that modeling conditions were changed as minimally as possible). As an unintended consequence of this procedure, very rough data (grain sizes 1.0 mm and 2.0 mm) at lower velocities (less than 4 .0 m/s) failed to match data; however, at higher velocities, matching was achieved. Data that failed to match was omitted from further analysis. Figure 3 3 Modeled shear stress versus measured shear stress (Crowley et al 2014) The grain diameters reported in Crowley et al. 2012a appear to be on the order of magnitude of sand grains, and the SERF (and other piston style devices) are designed to measure the erosion rates of cohesive material. However, previous

PAGE 63

63 research with the d often occurred during a cohesive erosion test (Crowley et al. 2012b). Use of a uniform sample (or shear stress) section with a high roughness could approximate an infinitely chunked speci men. In other words, it is not difficult to imagine a scenario where many surface approximates a sample with the roughness of a sand specimen surface. Likewise, previous rese arch has shown that the introduction of a small amount of cohesive material to a sand specimen may cause the specimen to erode like a cohesive material (Mitchener and Torfs 1996 for example). These specimens may have surface roughnesses on the order of ma gnitude of a typical, uniform sand specimen, but they may erode slowly like a cohesive soil and exhibit higher, cohesive style critical shear stresses. Next, mean values of r section. These data appear to correspond to previous research in that r appears to be approximately an order one multiplier of d (Einstein and El Samni 1949, Kamphuis 1974, Bayazit 1976, Dancey et al. 2000, Rahman and Webster 2005, Camanen et al. 2006, etc.). In fact, when data fr om this study was plotted alongside data in a similar range (Kamphuis 1974, Rahman and Webster 2005) results followed a very similar trend (Figure 3 4). Additionally, modeled pressure differentials were used to compute shear stresses via Equation 3 5 (Fig ure 3 5). Pressure differentials were computed just upstream and just downstream from the sensor (2012a) previous experiments. Total distance between pressure ports was approximately 4.0 in. (10.16 cm).

PAGE 64

64 Figure 3 4 Correlation between wall friction factor, f and grain size (Crowley et al 2014)

PAGE 65

65 Figure 3 5 Computed shear stress using pressure drop from CFD model for varying roughnesses (Crowley et al 2014) To verify these results, a convergence study was conduct ed using a representative flow rate (5.0 m/s), a smooth sample section, and a smooth sensor section. The refined portions of the mesh (across the sample section and across the sensor section), were replaced with several finer and coarser meshes. Average shear stress was computed across the sensor section and the sample section, and results were plotted as a function of cell size. A best fit, linear regression curve was established, and shear stresses for a cell size of zero were computed (i.e. a Richards on extrapolation). Results (Figure 3 6) appear to indicate that using average cell sizes of

PAGE 66

66 1.0 mm produces cell size computational errors less than 5.0%, which is why the ~0.5 mm grid size was deemed to be acceptable. Figure 3 6 Grid sensitivity study results (Crowley et al 2014) 3 2) was replaced with four configurations (Figure 3 7) to approximate differential erosion rates observed by Crowley et al. (2012b), Annandale et al. (2006), a nd recommendations from Briaud et shear stresses. Configurations were dubbed conical protrusion configuration, 1 mm

PAGE 67

67 protrusion configuration, differential erosion co nfiguration, and wavy configuration. These configurations and the rationale for their selection are described below. Figure 3 7 SERF sample section configurations showing (a) conical protrusion configuration; (b) 1 mm protrusion configuration; (c) dif ferential erosion configuration; and (d) wavy configuration. Flow is from left to right for all configurations. (Crowley et al 2014) 3.5.1.1 Conical p rotrusion As discussed in Crowley et al. (2012a), piston advancement in the SERF is controlled through a feedb ack loop between an ultrasonic depth sensor array and a Servo stepper motor. The ultrasonic array consists of twelve crystals that measure distance from the top of the flume to the top of the sample. When the average of these differences is within the to lerance of the stepper motor, an advancement signal is triggered. Sometimes, a sample may over advance because of an errant ultrasonic

PAGE 68

68 signal. When sand is tested and an error occurs, a small conical shaped slope as illustrated in Figure 3 5a is often ge nerated because sand grains tend to roll from the center of the eroding sample toward its edges. For this configuration, the conical sample protrudes 1.50 in. (3.81 cm) into the flume. 3.5.1.2 One mm p rotrusion As discussed, Briaud et al. (1999, 2001, 2004a, 2004b) recommend a 1.0 mm protrusion into their EFA device when cohesive samples are tested. While this will face, it was unclear how this affected Figure 3 7b shows this 1.0 mm protrusion. 3.5.1.3 Differential e rosion During previous SERF testing with cohesive sediment and rock, investigators often found that the upstream portions of the samples eroded m uch more slowly than their downstream portions. A photograph in Annandale (2006, p. 287) appears to show similar behavior during an EFA test. Figure 3 5c is meant to approximate this phenomenon. In this configuration, the front half of the sample is hel d flush with the flume bottom, while its back edge is recessed 0.50 in. (1.27 cm) below the flume floor and a linear plane is used to connect the two points. 3.5.1.4 Wavy scenario often observed during cohesive and rock tests with flume style devices. A random pattern of recesses was prepared to illustrate forcing on a very rough sample during an ero sion test. This configuration was tested in two different ways. First, its top

PAGE 69

69 edge was held level with the flume bottom. Then, average elevation was held flush with the flume bottom. Smooth bed assumptions were made for each configuration, and result s were compared with smooth, flush results. Amplification factors were computed by dividing non uniform configuration shear stresses by flush configuration shear stresses. Because the grid points did not align exactly with one another, an interpolation a lgorithm was used. Results are presented from Figure 3 8 through Figure 3 12. A table of average amplification factors was also prepared (Table 3 1). Table 3 1 Average amplification factors for each SERF configuration Configu ration Average Amplification Factor Conical protrusion 1.01 1 mm protrusion 0.98 Differential erosion 0.72 Wavy top portion level with flume bottom 0.37 Wavy average elevation level with flume bottom 1.16

PAGE 70

70 Figure 3 8 Amplification factors for conical protrusion configuration. Flow direction is along the x axis. Figure 3 9 Amplification factors for 1 mm protrusion configuration. Flow direction is along the x axis.

PAGE 71

71 Figure 3 10 Amplification factors for differential erosion configuration. Flow direction is along the x axis. Figure 3 11 Amplification factors for wavy configuration when top edge is held flush with flume bottom. Flow direction is along the x axis.

PAGE 72

72 Figure 3 12 Amplification factors for wavy configuration when average elevation is held flush with flume bottom. Flow direction is along the x axis. Finally, the sample portion of the SERF was roughened using the same roughness coefficients developed in Figure 3 2. Simulated flow was applied at 5.0 m/s to illustrate a re presentative worst case flow condition. Contour maps were generated of bed shear stress as a function of downstream flume distance to quantify shear stress development along the sample length as a function of roughness (Figure 3 13). Similarly, the cente rlines of each contour plot were compared with one another (Figure 3

PAGE 73

73 Figure 3 13 Local shear stress variations for uniformly rou gh samples (contour units are in Pa)

PAGE 74

74 Figure 3 14 Centerline shear stress variations for uniformly rough samples 3.5.2 Analysis Figure 3 5 appears to confirm previous results in that as roughness increases, pressure differential does not significantly increas e. This is believed to be caused by appears to solidify previous arguments that using p to estimate shear stress in piston style erosion equipment is ineffective. Results for variable testing configurations appeared to show that small changes in sample geometry may have large effects on localized shear stresses. For the case of protruding samples (conical protrusion configuration and 1 mm protrusion configurati on), average surface shear stress was similar to shear stresses for a flush specimen. However, as demonstrated in Figure 3 8, localized amplification factors for a

PAGE 75

75 protruding specimen may become quite large as much as 6.0 times that for a flush specimen This may affect localized erosion rates, and brings into question the validity of assuming a constant shear stress for the entire sample surface. The recessed specimens exhibit more significant problems. As samples become more recessed, average shear stress appears to become significantly affected. For the differential erosion configuration, average shear stress was on average only approximately 0.72 times the shear stress for a flush sample on average. However, as demonstrated in Figure 3 edge may reach values as high as 2.5 times the stress for a flush sample. Likewise, the shear stress amplification factor appears to reach a similar maximum at the point where the sample begins to slope downward. This would appear to indicate that if differential erosion was to occur over time (and as discussed, it has been repeatedly observed during testing), the issue will exacerbate itself. In other words, because differential erosion occurred, a higher downstream edge shear stress was produced in some locations. Ultimately, this may lead to higher downstream edge erosion rates, which in turn would produce higher downstream stresses at certain localized positions. Meanwhile, the presence of th ese high stress concentration regions appears to be somewhat balanced by other localized areas where shear stress is reduced. Thus, in sample geometry and turbulent flow fl uctuations. Figure 3 11, random recesses may produce large amplification factors, although on average, shear stress appears to decrease. Conversely, in Figure 3 12 where the

PAGE 76

76 s ample was on average level with the flume bottom the average shear stress on the sample surface appears to be much closer to smooth wall results. The presence of f localized stress amplification. Meanwhile, areas that failed to erode will continue to fail to erode because the stresses upon these regions of the sample are too low. Results in Figure 3 13 and Figure 3 14 appear to provide some As illustrated, when roughness is added to the flume (as would be the case when a sample is eroded), shear stress increases as a function of distance downstream. And, as roughness increase s, the average shear stress upon the sample becomes an increasingly poor estimator of actual stress conditions. All samples appeared to exhibit approximately equal localized shear stresses for 2.0 mm. Very rough samples (samples with equivalent grain siz es of 1.0 mm and 2.0 mm) appeared to exhibit a local shear stress maximum just beyond this upstream 2.0 mm edge; while intermediate roughness specimens (specimens with roughnesses of 0.125 mm, 0.25 mm, and 0.5 mm) showed maximum shear stresses closer the s shown in Figure 3 14, for very rough samples, once the localized maximum is reached, shear stress appears to level off and slightly decrease. This would appear to indicate that the sudden introduction of high roughness acts a flow system. Thus, adding high roughness generates enormous localized shear

PAGE 77

77 stresses that require some distance to resemble a more developed flow condition. This result appears to show, from a stress perspective, why and persist throughout an erosion event. And, it calls into question the accuracy of piston style erosion test. For conservative testing re sults, estimated average shear stress must be lower than actual sample shear stress, or estimated erosion rate must be higher than actual sample erosion rate for a given shear stress. If erosion rate and shear stress are correlated to one another, then th ese results provide interesting clues as to how to properly implement flume style erosion testing. First, it appears to be incorrect to keep the sample recessed relative to the flume bottom when differential erosion, blocking, or chunking occur because d oing so reduces shear stress in regions of the sample. This in turn may reduce localized erosion rates, which would lead to a non conservative erosion function. A better solution is to keep an eroding sample in a piston style device level with the flume bo ttom on average. Briaud et al. (1999, 2001, 2004a, 2004b) recommendation for a 1.0 mm protrusion appears to be close to correct, although even this recommendation may not go far enough. The Briaud et al. (1999, 2001, 2004a, 2004b) recommendation does no t take into account a situation where a localized portion of the sample erodes more than 1.0 mm. From the results in this paper, for a situation quantify average sample el evation, the most conservative method for an erosion test regardless of any protrusion into the device.

PAGE 78

78 CHAPTER 4 AN OVERLOOKED LOCAL SEDIMENT SCOUR MECHANISM 4.1 Governing Equations for Flow Around a Cylinder For flow around a circular cylinder, potential flow theory is only valid for flow scour applications, cylinder separation should be expec ted to occur at 90 degrees to the Following Sheppard (2004), for the definition sketch presented in Figure 4 1 for flow past a circular cylinder, the stream function may be written in polar coordinates as: (4 1) where is the stream function; is the depth mean velocity the approach flow; r is the radial coordinate; is the angular coordinate; and a is the radius of the circular pile (i.e. b/2). The radial and angular compone nts of velocity in terms of the stream function are then: (4 2) (4 3) where is the radial velocity component and is the angular velocity component. The magnitude of velocity is then: (4 4)

PAGE 79

79 pressure and velocity (assuming a constant elevation): (4 5) in which is the density of water; is the pressure at mid elevation of the approach flow; and v is the fluid velocity. Using Equation 4 5, the stagnation pressure ( ), or pressure along the leading edge of the pile may b e established: (4 6) Rewriting Equation 4 6 yields: (4 7) Likewise, the pressure at any point in the flow field may be expressed as: (4 8) This may be rewritten as: (4 9) Equation 4 4 and Equation 4 9 may be normalized by using the following nondimensional variables: (4 10) (4 11) (4 12 )

PAGE 80

80 Substituting results in: (4 13) (4 14) Next, the pressure gradient ( ) is defined as a vector with a magnitude equal to the maximum spatial change in pressure and a direction defined by that maximum change: (4 15) in which and are unit vectors in the radial an d angular directions respectively. Applying Equation 4 15 to Equation 4 14 yields: (4 16) The magnitude of the pressure gradient at any point in the flow field is then: (4 17) where N is defined as or the normalized coordinate in the direction of the maximum change of pressure in which n is the dimensional coordinate in the direction of the maximum change of pressure. Now, the forces acting on a sediment particle at the surface of the bed in the vicinity of the pile must be examined. The magnitude of the drag force ( F D ) upon this particle is: (4 18)

PAGE 81

81 where is the drag coefficient; A is the projected area of the particle; and u is the approac h velocity at the level of the particle. For the case of a spherical particle, A may be defined as where c is half the particle diameter ( D 50 /2). The force on the particle due to the pressure gradient ( ) may be expressed in differential form as: (4 19) where is the angular coordinate measured from the positive n direction (positive counterclockwise) and: (4 20) The component of this force in the n direction ( ) is then: (4 21) Integrating the component of the pressure force in the n direction over the particle surface results in: (4 22) Equation 4 22 represents an estimate of the effect of the structure induced pressure gradient on the particle. After some algebra, Equation 4 22 reduces to: (4 23) To find the relative magnitudes of the force on the particle due to the structure imposed pressure gradient by the force on the particle, must be divided by (i.e. dividing Equaiton 23 by Equation 18):

PAGE 82

82 (4 24) After some further algebraic manipulation: (4 25) Note the apparent dependency term between D 50 and b 4.2 Computational Modeling As discussed above, this method relies upon application of potential flow theory near the bed. As a pseudo verification of this analysis, a computational model of flow past two circular piles was prepared using CD CCM+ whereby one pile was prototype scale (3.0 m diameter) and the other pile was laboratory scale (2.0 cm diameter). Based upon Sheppard et al. (2004) and Sheppard and Miller (2006), the nondimensional parameters that govern equilibrium scour depth (beyond b/D 50 ) were assumed to be the ratios between water depth and pile width ( h/b ) and free stream velocity to critical velocity ( V/V c ). As such, water depth was varied to preserve the h/b ratio between one model and another. Note that the implication of this analysis is that V does not change from one model to an other because V c from laboratory to prototype scale cannot change. In other words, sand cannot be scaled because doing so would introduce cohesive forces into the laboratory model. Note that this assumption appears to be in conflict with several scour eq uations which rely upon the Froude Number to predict equilibrium scour depth. As discussed in Sheppard (2003) however, this apparent discrepancy is not at all unusual; after almost 40 years of scour research, engineers still cannot agree upon which non di mensional parameters actually govern scour. Details of the computational model are given below:

PAGE 83

83 4.2.1 Model Dimensions As discussed above, two piles were modeled one at prototype scale (3.0 m), and one at laboratory scale (2.0 cm). For the purposes of the m odel, a rectangular flow domain was prepared around the piles. Dimensions of the flow domain were 30 m wide by 15 m high for the larger pile and 0.2 m wide by 0.1 m high for the smaller pile. These dimensions were chosen to minimize side wall effects on the flow. Because of known deficiencies associated with Reynolds Averaged Navier Stokes (RANS) modeling associated with flow separation, two small perturbations (2 cm for the large pile and 0.16 mm for the small pile) were included in each model to induce separation (Figure 4 2 ). These were placed at on the pile at 90 degrees to trip separation exactly where potential flow theory suggests it will occur. It should be noted that th is is a separation modelling technique and nonphysical. Zdravkovich (1997) sho ws actual separation at locations closer to 110 degrees. Figure 4 1 D epiction of flow around a cylinder with associated variable representation in polar coordinates (Zdravkovich 1997)

PAGE 84

84 4.2.2 Mesh Geometry All meshing was performed within Star CCM+ using its hexagonal meshing a prism layer was included to improve computation accuracy near the bed. For the large model the layer consisted of 10 prismatic cells having a total thickness of 2 mm .For the small model the layer consisted of 5 prismatic cells having a total thickness of 0.05 mm. Additionally, two refinement planes were prepared bottom bou level. Approximately 5 Million and 10.87 Million cells were used for the larger and smaller model respectively. Along the refinement planes cell resolution as approximately 0.1 m per cell fo r the larger model and 0.5 mm per cell for the smaller model. These dimensions were significantly smaller near the boundaries due to the prism layer. Outside of the refinement planes, cell resolution was approximately 0.4 m per cell for the larger model an d 1 mm per cell for the smaller model. Illustrations of the associated smaller pile mesh is shown in Figure 4 3. The mesh depiction for the larger pile appears nearly identical with the except of a scaling factor Cell resolution yielded typical y values (near the pile and the bed) of 50 mm for the larger model and 0.5 mm for th e small model This corresponded to y + values less than one in both cases.

PAGE 85

85 Figure 4 2 Top view of pile with associated refinement region

PAGE 86

86 Figure 4 3 Top view of pile w ith separation inducing notches Figure 4 4 View of prism layers along bottom boundary

PAGE 87

87 4.2.3 Turbulence Model Formulation As stated above a realizable, two layer, k model with two layer all y+ wall treatment was chosen for this and all other, simulation s Section 2.2.3 provides an in depth discussion of the model choices and other relevant information. 4.2.4 Wall Treatment Formulation All model choices regarding the wall boundary formulation were the same as those chosen in Crowl e y et al. (2014). Section 3.4.4 provides an in depth discussion of the model choices and other relevant information. 4.2.5 Inlet and Outlet Formulation As discussed above, a rectangular flow regime was prepared around the modeled piles. Thus, the simulations wer e similar to flume style flow experiments around a circular pile. And, therefore, an inlet and an outlet were defined. CCM) is defined such that its velocity vectors are specified explicitly in the nor mal direction. Boundary face pressures are extracted from cells adjacent to the velocity inlet using a hybrid Gauss least square method (LSQ) reconstruction gradient. CCM+ reconstruction gradients. Outlet flow velocity is computed adding a nondimensionalized scale factor to the inlet velocity whereby the scale factor is a function of the mas s into the most downstream cells within the flow domain.

PAGE 88

88 4.2.6 Model Computational Scheme Initial conditions of each model were established such that the flow domain was filled with a stagnant mixture of air and water using a volume of fluid (VOF) model (A detailed describe of the model is laid out in Section 5.4.3.1) velocity was specified at the velocity inlet within the water portion of the fluid domain. Over time, water flow became fully developed nstream. Once average pressure at both the bottom on the modeled 4.3 Results Although analysis is still being conducted current results show a strong verification of the theory application by the numerical model. Simulation illustrations and associated pressure gradients for the small pile are presented below.

PAGE 89

89 Figure 4 5 Side view of model with flow (left to right)

PAGE 90

90 Figure 4 6 Pressure contour and gradient vectors for pile radius 0.02 m at bed Figure 4 7 Pressure gradient magnitude contours for pile radius 0.02 m at bed

PAGE 91

91 CHA P TER 5 WAVE IMPACTS ON BRIDGE DECKS 5.1 Background The third case study was that of quantifying vertical forces associated with wave impacts on bridge superstructures. Current research in the coastal engineering community is highly focused on the structural integrity of coastal structures and br idges (ASCE 2013 ). In parallel with structural integrity a major focus was placed on better understanding the overall mechanics of forces that bridge support structures encounter. A large percentage of this increase in investigating coastal forces has been allocated to examining the existence of slamming forces caused by high frequency wave impacts. Slamming forces are defined as the portion of a wave impact force that is not quantifiable as a quasi static loading. The frequency at which that force is measu red is regarded as the major factor in filtering the slamming force from that of the hydrostatic wave forces (Sheppard et al 2006). Throughout the hurricane seasons of 2004 and 2005 multiple bridge failed due to damage from wave impact s and storm surge. Examples include the I 10 Escambia Bay Bridge in Pensacola, FL; the Lake Pontchartrain Causeway just outside of New Orleans, LA, the Biloxi Bay Bridge in Biloxi, MS, the Mobile Bay onramp in Mobile, AL, and the St. Louis Bay Bridge just outside of Bay St. Louis, MS. Each of these structures were closed for extended periods as a result of bridge superstructure damage (Sheppard and Marin 2009) Additionally, as a consequence of cyclic economic situations and rising energy costs infrastru cture inspections ha ve become vu lnerable portions of local and national budgets. Bridges especially have encountered increased negligence. Consequently

PAGE 92

92 failures of aging bridges ar e at an all time high (ASCE 2013 ). As the overall age of the ncreas ed number of factors become relevant when predicting and analyzing failures. In 200 5 the Florida Department of Transportation (FDOT) responded to the mounting concerns over bridge integrity by commissioning the University of Florida Civil and Coast al Engineering Department to conduct wave forcing studies in their coastal engineering laboratory. Wave tank experiments (Figure 5 1) were conducted through 200 8 and a final report was submitted in December of 2009 ( Sheppard and Marin 2009 ) Figure 5 1 Physical wave force experiments conducted at t he University of Florida (Sheppard and Marin 2009 )

PAGE 93

93 Their findings indicated the presence of a significant slamming force and pointed to deficiencies in current bridge design regarding their wave impact integri ty. It was determined that although the physical experiments gave a significantly clearer picture of the physics involved further research was required to get a complete picture of the full nature of the problem 5.2 Experimental Data Set In an effort to addre ss the concerns regarding the unanswered questions left after the physical wave tank experiments concerning the slamming forces, a data set consisting of experimental wave height and vertical deck forcing ( as a function of time ) was obtained from the exper iments discussed above. The wave height data was collected with a wave gauge located 32 feet upstream of the forward most portion of the deck. The force data was collect with four three axis load /torque cells mounted on the deck structure. These data w ere the basis used for all numerical model comparison s Although the physical experiments analyzed multiple deck configurations of varied complexity, investigators determined that modeling should initially be conducted on one complex configuration for the sak e of efficiency. The experiment data set used is that of Transportation ( Sheppard and Marin 200 9 ) Plots of the data obtained are shown in figure 5 2 and 5 3 below. The values for run BSXX014 reported in the FDOT report are also listed in table.

PAGE 94

94 Table 5 1 Re ported values for run BSXX014 (Sheppard et al 2009) Report ed Values Value Water Depth (ft) 1.58 Wave Height (ft) 0.55 Period (s) 3 Vertical Force Max (lbf) 112.81 Figure 5 2 BSXX014 Experimental wave height s ignal

PAGE 95

95 Figure 5 3 BSXX014 Experimental vertical force s ignal During the investigation researchers noticed the presence of frequency interference effects located in the wave. Therefore a spectral analysis was conducted to (Figure 5 Although large scale interference was determined to be minimal investigators determined that althou gh the report states that the waves had a period of 3 seconds the spectral analysis appears to suggest that the period was actually 2.7 sec onds As a result the period used for the numerical analyses was change d to 2.7 seconds. All other reported values were taken as valid.

PAGE 96

96 Figure 5 4 Spectral a nalysis of w ave s ignal s consisted of a 1 inch thick polypropylene slab with seven, evenly spaced, girders attached the underside. The deck was constructed of three smaller slabs. Therefore the total transverse width was 72 inches from the three 24 inch sections All other dimensions and configuration layout is presented in figure 5 5.

PAGE 97

97 Figure 5 5 Schematic of deck configuration from p hysical e xperiment (Sheppard et al 2008) 5.3 Modeling Overview To better understand the mechanisms involved in the w ave forcing the previously mentioned experiments, conducted in the University of Florida wave tank were computationally modeling using CD adapco's Star CCM+. Data taken from the experiments was then used to calibrate the computational model. First, the wave signal itself was extensively studied using three wave generation methods by using mesh morphing to model a moving wall boundary sim ilar to a typical piston style wavemaker; implementation of linear wave theory; and implementation of fifth order wave theory. Ultimately, the oscillating fluid motion method was discarded because the latter three methods appeared to reproduce the wave si gnal more accurately. Finally, the bridge was rotated to study the effects of attack angle for bri dges subjected to wave attack.

PAGE 98

98 5.4 Computational Model 5.4.1 Model Geometry Because the simulation was to reproduce a physically conducted experiment every effort wa s made to match the dimensions of the actual wave tank in the University of Florida coastal laboratory as closely as possible. Initially a three dimensional rectangular domain was developed entirely with Star CCM+. The fluid domain was a rectangular box 150 feet (45.72 m) long by 6 feet (1.8288 m) wide by 10 feet (3.048 m) tall. 68 feet downstream of the front tank boundary the model bridge deck was placed into the domain in a manner which did not allow fluid to flow around the deck sides. Fluid would then only be permitted to flow over or under the modelled deck. The final geometric version of the deck and associated boundaries is shown in f igure 5 6. Figure 5 6 Model bridge deck before m eshing

PAGE 99

99 5.4.2 Polyhedral Meshing Scheme After the initial geometry was assembled meshing was performed exclusively within Star CCM+ using its polyhedral meshing scheme. A polyhedral mesh consists of cells containing an arbitrary number of faces that can be oriented in any direction The advantage to this meshing scheme is that in principle, it should allow for flow calculations normal to each cell face Due to the omni directional nature of most practical i nvestigations this is not always possible when using standard cell formulations (Davis et al 2012 ). The increased number of faces polyhedral meshes exhibit, especially compared with hexahedral and tetrahedral, increase s the probability that a particular ce ll face will be in or close to, the ideal flow orientation. Within Star CCM+ polyhedral cells average 14 cell faces, compared to 6 faces in a hexahedral scheme and 4 cell faces in a tetrahedral formulation CD a d apco states that in cases of complex flows and geometry polyhedral meshing schemes generally exhibit great levels of accuracy when compared to hexahedral and tetrahedral meshes of the same size (CD adapco 2012). It is noted that r Research and Analysis Computing Center (TRACC) have shown that in most cases the difference between hexahedral and polyhedral meshes are small (Lottes et al. 2013) Therefore although a polyhedral mesh was chosen for this simulation investigators believe a hexahedral mesh could have been used to produce similar results The mesh for all simulations consisted of ~6.3 million polyhedral cells. The average cell siz e used for the air region was 50 m m. In an effort to increase accuracy in areas of interest, a ref inement region was developed around the bridge deck and wave region which brought the average cell size to 5 mm. It should be noted that due to the arbitrary nature of polyhedral meshes cell sizes are estimates and varied somewhat through the mesh.

PAGE 100

100 A typi cal, non boundary, polyhedral cell used within the analysis is shown below in figure 5 7 Figure 5 7 A typical polyhedral cell u sed in the wave tank m esh 5.4.3 Volume of Fluid (VOF) Model T o simulate the effects of a multicomponent multiphase flow investiga tors invoked a volume of fluid (VOF) model to capture the significant difference in fluid properties between the two fluids being simulation (water and air). The VOF model is generally consider ed to be a somewhat simple approach to the multiphase problem a s it only accounts for the difference in fluid parameters and says nothing of the interaction between them. The basic assumption of the model is that both fluids share the same field values when evaluating velocity, pressure, and temperature terms. This as sumption generally only cause s there is significant fluid entrainment. It is noted that the VOF model can still produce

PAGE 101

101 reasonable results in situations where there is significant mixing of the fluids howev er this requires a fairly fine mesh as the errors associated with the fluid interactions are that of the discretization error. Figure 5 8 illustrates examples of when the VOF might produce significant errors on coarse meshes. Figure 5 8 Depiction of in valid (a) and valid (b) use of the VOF m odel (CD adapco 2012) Taking note of the assumption s la id out above the equations governing the volume of fluid model of the i th phase are: And

PAGE 102

102 The modified parameters are then used to create the governing transport equation for m ultiphase flow as Where is the phase source or sink term, and is the material derivative of the phase densities (for compressible flow) Figure 5 9 Volume of fluid depict ion on polyhedral mesh around deck 5.4.4 Moving Boundary Method The wave signal for the experimental data described was generated using a piston style wave maker device based on linear wave maker theory The theory put simply is that a sinusoidal motion pattern of the piston should result in the production of a linear wave moving downstream from the wave maker. For a more detailed discuss ion

PAGE 103

103 and derivation of corresponding governing equations the read er is referred to Water Wave Mechanics for Engineers and Scient ists ( Dean and Dalrymple 1984 ) and Waves Generated By a Piston Type Wavemaker (Madsen 1970) From the theory investigators determined that the piston motion should be governed by an equation of the form : (5 1 ) Where a b c and d are arbitrary empirically determined constants and t is simulation time. Within the CFD model investigators employed a numerical techni que known as mesh morphing to allow the boundary designated as the modelled piston to move relativ e to the other portions of the mesh. In the mesh morphing algorithm cells in the mesh are deformed after every timestep based on the prescribed boundary motion. This is done internally in Star CCM+ by moving the mesh vertices profile After each time step cells are either stretched or compressed depending on their A displacement vector is then determine d based the distance between the vertices old time step and new time step location. This is illustrated in the Star CCM+ manual as : Figure 5 10 D epiction of vertex motion ( Cd adapco 2012 ) i n which x 1 is the vertex location at the old time step, x 2 is the vertex location at the new time step, and d i is the displace ment of the specified vertex. A multiquadratic

PAGE 104

104 interpolation algorithm then interpolates flow properties, such as velocity and pressure The entire morphing process is the broken d own into five steps describe by CD a dapco as : 1. Collection of vertex control points and their prescribed displacements 2. A process known as thinning is applied to reduce the number of vertices needed to calculate the interpolation field. This is done to decrea se the time required by the morpher. 3. The calculated displacements are then used to generate the interpolation field base on the following system of equations (One equation for each vertex moved): (5 2) Where (5 3) (5 4 )

PAGE 105

105 4. The system of equations described in step 3 is then solved and the interpolation field for every vertex in the mesh can be found based on its displacement equation, from (5 5) Where can be applied to every vertex in the mesh 5. Final adjustments then are made to cells close to or at boundaries to ensure a nonintersecting mesh In the wave tank simulation the using the function described by E quation 5 1 while the mesh used to resolve the flow (i e excluding the bridge deck) was stretched and compressed in response to the moving wall boundary. This type of mesh motion is designated internally within Star CCM+ as 5.4.5 Stokes Wave Theory Method The other approach chosen for wave generation was to use the internal wave model within Star CCM+. This mode l follows Stokes Wave T heory and enables the user to generate either first order (linear) wave approximations or fifth order approximations. In this approach the boundary condition s we re modified from those stated in the moving boundary method. The front w (described in section 3.4.2) and the back wall of the tank was made a The pressure outlet boundary is similar to the however a pressure is inst ead defined as a known quantity. The formulation of the fluid velocity and pressure terms is then defined based on the Stoke Theory.

PAGE 106

106 5.4.5.1 First order (linear) wave f ormulation The first order theory used within Star CCM+ presented as Wave Horizontal Velocity: Wave Vertical Velocity: Wave Surface Elevation: Where a is the wave amplitude, is the wave frequency, is the vector, K is the magnitude of the wave vector, and z is vertical distance from the mean water level. The wave period, T is then defined by And the wavelength, is defined as 5.4.5.2 F ifth o rder The fifth order stoke wave is given as the expansion of the Stokes W aves series out to the fifth order terms. The full solution to the series expansion out to five terms can be cumbersome and is somewhat out of the scope of this paper. The read er is therefore referred to A Fifth Order Stokes theory For Steady Waves (Fenton 1985) for a full description of the model parameters.

PAGE 107

107 5.5 Analysis and Results Once an adequate wave signal had been generated, a comparison was made between experimentally obtained force data and data from the computational analysis. Additionally, downward force was not properly simulated because the computationally This was accomplished by adding an outlet boundary along the water washed onto the top portion of the bridge, it was allowed to drain immediately. The parameters found to gi ve the best results for the piston style wavemaker the linear theory method, and the fifth order method are shown in Table 5 2, Table 5 3, and Table 5 4 respectively. Table 5 2 Values for moving boundary method equation Variable for Eqn. Value a 0.45 b 1 c 1 d 20 The parameters found to give the best results for the first order wa ve approximation are give n Table 5 3 Values for linear wave theory model Wave Parameters Value Water Height 0.4816 m

PAGE 108

108 Wave Amplitude 0.08382 m Wave Period 2.7 s The parameters found to give the best results for the fifth order wave approximation we re found to be Table 5 4 Values for fifth order wave theory model Wave Parameters Value Water Height 0.4816 m Wave Height 0.16764 m Wave Period 2.7 s Although parameters used within the Stokes wave model were slightly different than those reported by Sheppard et al. (2009) they did fall within the range of acceptable values based upon the sensitivity of the instrumentation used. After each simulation results f ro m the simulated wave signal was compared to the physical signal. Because values upon from the simulation did not temporally line up with the experimental values an interpolation algorithm was used to allow for a direct comparison.

PAGE 109

109 Figure 5 1 1 Compariso n of wave generated by moving boundary method (red) with experimental (blue) wave Figure 5 1 2 One to one comparison of wave generated by moving boundary method with experimental wave

PAGE 110

110 Figure 5 1 3 C omparison of vertical force generated by moving boundary method ( green ) with experimental (blue) wave Figure 5 1 4 C omparison of wave generated by linear wave theory method (red) with experimental (blue) wave

PAGE 111

111 Figure 5 1 5 O ne to one comparison of wave generated by linear wave th eory method with experimental wave Figure 5 1 6 C omparison of vertical force generated by linear wave theory method ( red ) with experimental (blue) wave

PAGE 112

112 Figure 5 1 7 C omparison of wave generated by fifth order wave theory method (red) with experimental (blue) wave Figure 5 1 8 O ne to one comparison of wave generated by fifth order wave theory method with experimental wave

PAGE 113

113 Figure 5 1 9 C omparison of vertical force generated by fifth order wave theory method (red) with e xperimental (blue) wave Associated regression coefficients of the one to one line discussed above ( between modeled and physical wave signal data ) were 0.81, 0. 88 and 0.8 3 for the piston method, the linear method and the fifth order method respectively (Shown in figures 5 13, 5 15, and 5 17) The best method in terms of both wave signal and force magnitude was determined to be the linear wave theory method which predicted the wave force magnitude within a factor of 0.96 6 ( Table 5 5) However, although the forcing pattern was reproduced with high levels of accuracy both in terms of amplitude and period, the physical model appeared to behave more non linearly in time than any of the modeled results.

PAGE 114

114 Table 5 5 Wave forces modelled vs experimental Fo rces Value Experimental 113 lbf Simulation 117 lbf Amplification Factor 0.966 5.6 Further Work To Be Conducted The current work represents only a portion of the analysis that should be conducted to gain a full picture of this phenomena. Further will therefore continue i n following areas: 1. A signal analysis should be conducted on the simulation results and optimized to ensure the wave being produced within the model is as numerically close to the physical wave as possible 2. The bridge deck should be rotated to be investigate the effects of wave angle of attack upon forces 3. Once the fluids model has been appropriately opt imized and deemed a close to the physical experiments as possible a structure analysis should be conducted to examine the effects of the fluid forces upon the structural integrity of the deck slab.

PAGE 115

115 CHAPTER 6 SUMMARY AND CONCLUSIONS The Sediment Er osion Rate Flume (SERF) at the University of Florida was Computationally Modelled Computational model were able to reproduce laboratory experiments with statistically significant accuracy and precision The computational model revealed an otherwise unobse rved mechanism for the blocking and chucking of a soil specimen during an erosion test Investigators believe that care must be taken to determine conservative estimate of scour when a piston style device is being use for design. Additionally, the most con servative method for an erosion test would be to keep protrusion into the device. A new application of Potential Flow Theory was used to explain an overlook ed mechanism fo r bridge scour. The use of Potential Flow Theory draws attention to the dependency to sand grain size when determine the maximum scour hole size The computational model reinforce d these find ing s and showed that pressure gradients around the base of a pile are adequate to invoke the new theory application. Wave tank test s conducted by D. Max Sheppard and Justin Marin (2009) were computationally mode l l ed in an attempt to reproduce wave signal and forcing data Three method s were used to generate the waves The most effective was to invoke a linear wave theory model Although there is more work to be conducted wave signal comparison show a strong correlation (R 2 =0.9) and wave forcing magnitude comparison are within the statistically significant range (0.96 6 ) Continued work will be conducted in an effort to determine forcing patterns as a function of wave angle of attack.

PAGE 116

116 LIST OF REFERENCES American Society of Civil Engineers (ASCE) (2013). Infrastructure Washington D.C. Anderson, J. (1995). Computational Fluid Dynamics McGraw Hill, New York. Annandale, G.W. (2006). Scour Technology D. Ing, Ed. McGraw Hill, New York, NY. Argonne National Laboratory (2013). Tools and Technology for Solving Critical Scientific Problems Chicago. Arneson, L. A., Zevenbergen, L. W., Lagasse, P. F., and Clopper, P. E. (2012). Evaluating Scour at Bridges HEC 18 Fifth Edition U. S. Department of Transportation, Federal Highway Administration, Washington D. C. Bardina, J., Huang, P., C oakley, T. (1997) Turbulence Modeling Validation, Testing, and Development. National Aeronautics and Space Exploration Administration. Moffett Field, CA. surface J. of Hydraulic Re s. 14, pp. 115 126. Biswas, G., Eswaran, V. (2002) Turbulent Flows Fundamentals, Experiments, and Modeling. CRC Press Boca Raton. Bloomquist, D., and Crowley, R. of Erosion Rates of Rock, Sand, and Clay Mi SERF Equipment. Final Report No. BDK75 977 09, Florida Department of Transportation, Tallahassee, FL. EFA method for complex piers in fine gr J. of Geotechnical and Geoenvironmental Engineering 130(11), 1180 1191. Briaud, J. L., Ting, F.C.K., Chen, H. C., Gudavalli, R., Perugu, S., and Wei, G. (1999). SRICOS J. of Geotechnical and Geoenvironmental Engineering 125(4), 237 246. Briaud, J.L., Chen, H. C., Li, Y., Nurtjahyo, P., and Wang, J. (2004a). Pier and contraction scour in cohesive soils NCHRP Report 516 Transportation Research Board, Washington, D.C. Briaud, J.L., Ting, F.C.K., Chen, H. C., Cao, Y., Han, S. W., Kwak, K.W. (2001). J. of Geotechnical and Geoenvironmental Engineering, 127 (2), 105 113.

PAGE 117

117 Camenen, B., Bayram, A., and Larson, M. J. Hydraulic Engineering 132(11), pp. 1146 1158. Cd adapco (2012). User Guide Star CCM+ Version 7.06 Cd adapco, Melville, NY. Celik, I. (1999) Introductory Turbulence Modeling Lecture Notes West Virginia University Department of Mechanical and Aerospace Engineering. Morgantown. Crowley, R. W., Bloomquist, D. B., Hayne, J. R., Holst, C. M., and Shah, F. D. (2013). osion rate J. of Hydraulic Engr. 138(11), 990 994. Crowley, R. W., Bloomquist, D. B., Shah, F. D., and Holst, C. M. (2012 The sediment erosion rate flume (SERF) : a new testing device for measuring erosion rate and shear stress Geotechnical Testing Journal 35 (4) 649 659. Crowley, R., A Description of Erosion Rate Testing Devices and Correlations Between Rock and Erosion Rate and Cohesion International Conference on Scour and Erosion P aris. Crowley, R., Robeck, C. Thieke, R. (2014). Computational Modeling of Bed Material Shear Stresses in Piston Type Erosion Rate Testing Devices J. Hydraul. Eng ., 140(1), 24 34 Dancey, C. L., Balakrishnan, M., Diplas, P., and Papanicolaou, A. N. (200 spatial inhomogeneity of turbulence above a fully rough, packed bed in open Exp. Fluids 29, pp. 402 410. based turbulence modeling for flow over a wall m Proc. 20 th Annual Conference of the CFD Society of Canada Canmore, AB, May 9 12. Dean, R. G, Dalrymple (2000). Water Wave Mechanics for Engineers and Scientists World Scientific Publishing Company. London. Rev. Modern Physics 21, pp. 520 524. J. Geophysic al Research 67(4), 1451 1461. Ettema, R. (1980). Scour at bridges. Report No. 216, University of Auckland, School of Engineering, Auckland, New Zealand. Fenton, John D. ( 1985 ) A Fifth Order Stokes Theory for Steady Waves J. Waterway, Port, Coastal and Ocean Eng ., 111 (2), pp. 216 234.

PAGE 118

118 Froehlich, D. C. (1988). Analysis of on site scour at bridge piers. Proc. ASCE National Hydraulic Engineering Conference American Society of Civil Engineers, Reston, VA. Garde, R.J., Ranga Raju, K.G., and Kothyar i, U.C. (1993). Effect of Unsteadiness and Stratification on Local Scour International Science Publisher, New York. Hanson, G. J., and Cook, K. R. (2004). Apparatus, test procedures, and analytical methods to measure soil erodibility in situ." Applie d Engineering in Agriculture 20(4), 455 462. Jones, W.P., Launder, B.E. (1972) The Prediction of Laminarization with a Two Equation Model of Turbulence Int. J. Heat and Mass Transfer 15 pp. 301 314. Thesis, Lausanne EPFL. Kamphuis, J. W. (1974). Determination of sand roughness for fixed beds. J. Hydraulic Res., 12, pp. 193 203. Kolmogorov, AN. (1941a). The local structure of turbulence in incompressible vis cous fluid for very large Reynolds number SSSR 30: 9 13. doi:10.1007/978 94 011 3030 1_45. Lottes, S. (2012). Principle Mechanical Engineer, Energy Systems Division, Argonne National Laboratories Transportation Research and Analysis Computing Center. Pe rsonal conversations. Madsen, O. (2011). W aves generated by a piston type wavemaker Coastal Engineering Proceedings 1 (12). doi:10.9753/icce.v12. J. of Hydraulic Engineering 122(6), pp. 316 324. Melville, B.W. and Sutherland, A.J. (1989). Design method for local scour at bridge pi ers. Journal of Hydraulic Engineering, 123 (2), 125 136. Menter, F.R. ( 1994 ) Two equation eddy viscosity turbulence modeling for engineering applications AIAA Journal, 32(8), pp. 1598 1605. Coastal Engineering 29, pp. 1 25. Moin Mahesh. ( 1998 ) Direct Numerical Simulation: A Tool in Turbulence Research Annual Review Fluid Mechanics 1998 30:539 78 J. of the H ydraulics Division, ASCE 91(1), 105 138.

PAGE 119

119 Patankar, S. V. and Spalding, D.B. (1972) A calculation procedure for heat, mass and momentum transfer in three dimensional parabolic flows Int. J. of Heat and Mass Transfer Volume 15, Issue 10, October 1972 pp. 1787 1806. Pope, S. (2000). Turbulent Flows Cambridge University Press. Experiments in Fluids 28, pp. 372 384. Reichardt, Z. Andrew. Math. Mech 31(7), pp. 208 219. bulk density o J. of Hydraulic Engineering 124 (4), pp. 1261 1267. layer models combining the k e model with a one 29 th Aerospace Sciences Meeting January 7 10, Reno, NV, AIAA 91 0216. Sheppard, D.M, and Marin, J. (2009). Forces on Bridge Decks Final Report No. BD545 58, Florida Department of Transportation, Tallahassee, FL. Sheppard, D.M. (2003). Scour at complex piers. Final Report, Florid a Department of Transportation, Tallahassee, FL. Sheppard, D.M. (2004). Overlooked local sediment scour mechanism. Transportation Research Record: Journal of the Transportation Research Board 1890, TRB, National Research Council, Washington, D.C., 107 111. Journal of Hydraulic Engineering 132 (7), 635 642. water local pier Journal of Hydraulic Engineering 130 (10), 957 963. Proc. 1 st International Conference on Water Resources Engineering Sand Antonio, TX, August 14 18, vol 2, 1809 1813. Shih, T. eddy viscosity model for high Reynolds number turbulent flows model development Smagorinsky, J ., (1963). General circula tion experiments with the primitive equations. Mon. Wea. Rev ., 91, pp. 99 164.

PAGE 120

120 Spalart, A one equation turbulence model for aerodynamic flows AIAA 92 0439. Trammell, M. (2004). Laboratory Apparatus and Methodolgy for Determi ng Water Erosion Rate of Erodible Rock and Cohesive Sediments. MS Thesis. University of Florida. Gainesville, FL. Versteeg, H., Malalasekera, W. (2007). An Introduction to Computational Fluid Dynamics Pearson Educatio n Limited. Weidner, K. (2012). Evaluation for the Jet Test Method for Determining the Erosional Proerties of Cohesive Soils: A Numerical Approach. MS Thesis. Virginia Polytechnic Institute and State University. Blacksburg, VA. White, F. (1986). Fluid mechanics, McGraw Hill, New York, 3 02 304. Wilcox, D.C. ( 199 8). Turbulence Modeling for CFD 2nd edition, DCW Industries, Inc. dimensional flow International Journ al of Heat and Mass Transfer 12, pp. 301 318. Zdravkovich, M.M (1997). Flow Around Circular Cylinders Volume I, Fundamentals Oxford University Press. New York.

PAGE 121

121 BIOGRAPHICAL SKETCH Corbin Robeck was born in McLean, Va in 1989, T he second son of a medical doctor and a history teacher. He came to the University of Florida in 2007 initially to pursue a biochemistry degree but transferred to engineering midway through his sophomore year Th r ough his undergraduate studies he was a membe r of multiple student clubs on campus including: the Gator Marching Band, the University Jazz Orchestra, and the UF Steel Bridge Team He graduated with a B. S. in c ivil e ngineering in 2012. When presented with the opportunity to join the industry he inst ead chose to remain at UF and enroll in concurrent graduate degree programs (civil and m echanical engineering). As a graduate student in the D epartment of C ivil and Coastal Engineering and Department of M echanical and Aerospace Engineering he w as exposed to a variety of new experiences. One of the more instructive experiences was his work as a collaboratory Computing Center (TRACC) in Chicago. During his time at Argonne Co rbin was under the tutelage of Dr. Steve n Lott es and Dr. Cezary Bojanowski. His work at Argonne consisted primarily of computational fluid structure interaction related to infrastructure problems Scientifically, h e publish ed his first conference paper in mid 2012 and his first journal paper in the summer of 2013 He hopes to publish two additional papers by the time of his graduation in May of 2014. Upon completion of his M.S. degrees Corbin hopes to enroll in a PhD program to continue his work in the fiel d of fluid structure interaction.