UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 1 NUMERICAL INVESTIGATION OF PARTICLE LADEN FLOW IN OSCILLATORY CHANNEL AND ITS IMPLICATION TO WAVE INDUCED FINE SEDIMENT TRANSPORT By CELALETTIN EMRE OZDEMIR A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2010 PAGE 2 2 2010 Celalettin Emre Ozdemir PAGE 3 3 To my mother light in my eyes when it is dark my mind when I cannot think, my legs when I cannot walk; and to my father my source of strength and greatest teacher, even though he is not with me, may his soul rest in peace. PAGE 4 4 ACKNOWLEDGMENTS I would like to thank Dr. Tian Jian H su not only for his supervision of this work but a lso for being a great mentor, a friend and also for providing a good research environment through my doctoral studies. I also would like to present my heartfelt and sincere thanks to Dr. S. Balachandar for his advice the discipline he instilled and his contributions to my rapid intellectual development for the last two years. I am also grateful to Dr. Ashish Mehta and Dr. Mamta Jain for their favors which cannot even be paid by my lifelong loyalty to them, and reminding me there are good people in the world I also would like to present my thanks to my committee members Dr. Kirk Hatfiled, Dr. Louis Motz, Dr. Ashish Mehta, and Dr. John Jaeger for their support. My friends in computational multi physics group CMG, deserve special thanks not only for their friendship but also for their constructive comments through the course of this study Especially I would like to present my thanks to Dr. Mariano I. Cantero, who is now assistant professor in University of Bariloce, Argentina, and Dr. Thomas Bonometti, who is now assistant professor in Institut National Polytechnique de Toulouse, France, for their help and support at the initial stage of my dissertation. Needless to say, I am thankful to my family: my mother Zeliha Ozdemir, my father Cumhur Ozdemir may his soul rest in peace and my sister Dr. Nilufer Aslihan Ozdemir. I also thank my dear friend Dr. Ahmet Emr e Tekeli and anyone who does not lack my well being and success in their prayers. This study is supp orted by the U.S. Office of Naval Research (N00014091 0134) and National Science Foundation (OCE 0913283) to University of Delaware and by U.S. Office of Naval Research grant (N00014 071 0494) to University of Florida. This work is partially supported by the National Center for Supercomputing Applications under OCE70005N and PAGE 5 5 OCE80005P utilized the NCSA Cobalt and PSC Pople. Also, the author acknowledges the University of Florida High Performance Computing Center for providing computational resources and support that have contributed to the research results reported in this thesis. PAGE 6 6 TABLE OF CONTENTS page ACKNOWLEDGMENTS .................................................................................................................... 4 LIST OF TABLES ................................................................................................................................ 8 LIST OF FIGURES .............................................................................................................................. 9 ABSTRACT ........................................................................................................................................ 12 CHAPTER 1 INTRODUCTION ....................................................................................................................... 14 1.1 Overview ............................................................................................................................... 14 1.2 Literature Survey................................................................................................................... 18 1.2.1 Turb ulence Modulation due to Fine Particles ........................................................... 19 1.2.2 Boundary Layer Turbulence in Oscillatory Flow ..................................................... 20 1.3 Scope of the Thesis ............................................................................................................... 23 2 PROBLEM STATEMENT AND MATHEMATICAL FORMULATION ............................ 25 2.1 Problem Statement ................................................................................................................ 25 2.2 Mathematical Formulation ................................................................................................... 26 3 NUMERICAL METHODOLOGY ............................................................................................ 32 3.1 Overview ............................................................................................................................... 32 3.2 Pseudo Spectral Methods ..................................................................................................... 34 3.3 Solution of Governing Equations ......................................................................................... 38 3.4 Solution of Helmholtz and Poisson Equations .................................................................... 41 4 CLEAR FLUID SIMULA TIONS: PASSIVE CASE ............................................................... 44 4.1 Grid Resolution and Domain Size ....................................................................................... 45 4.2 Comparison of Mean and Turbulent Quantities .................................................................. 47 5 TWO WAY COUPLED SIMULATIONS: EFFECT OF PARTICLE VOLUME CONCENTRATION ................................................................................................................... 55 5.1 Laminar Solution................................................................................................................... 55 5.2 Turbulent Structures ............................................................................................................. 57 5.3 Mean Quantities and Turbulent Intensity ............................................................................ 60 5.4 Flow Instability ..................................................................................................................... 64 5.4.1 Stability of Stokes Layer ............................................................................................ 64 5.4.2 Effect of Density Stratification .................................................................................. 66 PAGE 7 7 5.5 Flux and Gradient Richardson Numbers ............................................................................. 68 5.6 TKE Budget ........................................................................................................................... 71 5.7 Flow Regimes a nd Field Interpretations .............................................................................. 75 5.7.1 Regime 1 ..................................................................................................................... 75 5.7.2 Regime 2 ..................................................................................................................... 76 5.7.3 Regime 3 ..................................................................................................................... 76 5.7.4 Regime 4 ..................................................................................................................... 78 6 TWO WAY COUPLED SIMULATIONS: EFFECT OF PARTICLE SETTLING VELOCITY ................................................................................................................................. 94 6.1 Qualitative Observations: Turbulent Structures .................................................................. 95 6.2 Mean and Turbulent Quantities ............................................................................................ 98 6.3 Mixing and Suspension ...................................................................................................... 100 6.4 TKE Budget ......................................................................................................................... 103 7 CONCLUSIONS ....................................................................................................................... 120 LIST OF REFERENCES ................................................................................................................. 127 BIOGRAPHICAL SKETCH ........................................................................................................... 133 PAGE 8 8 LIST OF TABLES Table page 2 1 List of simulations performed in this study. ......................................................................... 31 PAGE 9 9 LIST OF FIGURES Figure page 1 1 Illustrative sketch for the processes involved in sediment source to sink. ......................... 23 1 2 Temporal and spatial scales for different sediment transport processes ............................. 24 2 1 Two dimensional descriptive sketch of the problem. The velocity profile is the laminar velocity profile at t .................................................................................... 31 3 1 Time series of friction velocity for (A) Ri = 0.0, (B) Ri = 4110 and (C) Ri = 4310 .................................................................................................................................... 43 4 1 Two point correlation functions in spanwise direction at y+ = 0.2 and y+ = 78 for phases (A) ................................................................. 50 4 2 Two point correlation functions in streamwise direction at y+ = 0.2 and y+ = 78 for phases (A) C) .................................................................. 50 4 3 Power spectrum of fluctuating velocities in streamwise and spanwise direction for phases at y+ = 0.2 A) ..................................................... 51 4 4 Power spectrum of fluctuating velocities in streamwise and spanwise direction for phases at y+ = 78 A) ..................................................... 51 4 5 Phase and planar average of streamwise velocity ( U ) compared with that of Spalart & Baldwin (1989) (SB1989). ................................................................................................ 52 4 6 Phase and planar average of streamwise turbulent intensity ( Urms) compared with that of Sp alart & Baldwin (1989) (SB1989). ............................................................................... 52 4 7 Phase and planar average of vertical turbulent intensity ( Vrms) compared with that of Spalart & Baldwin (1989) (SB1989) .................................................................................... 53 4 8 Phase and planar average of spanwise turbulent intensity ( Wrms) compared with that of Spal art & Baldwin (1989) (SB1989) ................................................................................ 53 4 9 Vortex structures identified by swirl intensity ( ci)2 (Zhou, et al. 1999) at different phases ..................................................................................................................................... 5 4 4 10 Iso contours of concentration at C = 1.3 at different phases. .............................................. 54 5 1 Turbulent vortical structures extracted in terms of iso surfaces of swirling strength (ci).for cases (A) Ri = 0, (B) Ri = 104, and (C) Ri = 3 x104. ............................................. 80 5 2 Turbulent vortical structures extracted in terms of iso surfaces of particle concentration for case 1, case.2, and case3. ......................................................................... 81 PAGE 10 10 5 3 Mean concentration profiles for Ri = 0.0, 104 and 3104.. ................................................. 82 5 4 Mean streamwise velocity profiles for Ri = 0, 104 and 3104. .......................................... 83 5 5 Turbulent rms velocity fluctuation along (A) streamwise, (B) wall normal and (C) spanwise directions. ............................................................................................................... 84 5 6 Iso surfaces of concentration at six different phases from case 2.. ..................................................................................................................................... 85 5 7 Time variation of extremum velocity for Ri = 104 and Ri = 3 x 104. .................................. 86 5 8 Spanwise vorticity contours and time series of maximum and minimum perturbation vortcity for case 3. .................................................................................................................. 86 5 9 Plot of gradient Richardson number at six different phases for (A) Ri = 104 and (B) Ri = 304 in log scale. ......................................................................................................... 87 5 10 Plot of flux Richardson number at six different phases for (A) Ri = 104 and (B) Ri = 3104 in log scale. ................................................................................................................. 88 5 11 The ratio of turbulent particle suspension flux over settling flux for (A) Ri = 104 and (B) Ri = 304. ...................................................................................................................... 89 5 12 Profiles of TKE budget terms for Ri = 0.0. .......................................................................... 90 5 13 Profiles of TKE budget terms for Ri = 0.0 and Ri = 104 at y = 15~25. .............................. 91 5 14 Profiles of TKE budget terms for Ri = 104. ......................................................................... 92 5 15 Profiles of TKE budget terms for Ri =3 04. ..................................................................... 93 6 1 Representation of simulation parameters in parametric space.. ........................................ 107 6 2 Iso contours of concentration for case 5 and case 6. ......................................................... 108 6 3 Iso contours of swirl strength for case 5 and case 6. ......................................................... 109 6 4 max min) for case 6.. ............................................................................................................................. 110 6 5 Ensemble averaged velocity profiles for case 2, case 3, case 5, and case 6.. .................. 111 6 6 Ensemble averaged concentration profiles of case 2, case 3, case 5, and case 6.. ........... 112 6 7 r.m.s of velocity fluctuation profiles for case 2, case 3, case 5, and case 6 at six different phases. .................................................................................................................... 113 PAGE 11 11 6 8 Flux Richardson number and gradient Richardson number profiles for A) case 5 ( Vs = 3 104 and Ri = 1 104 ) and B) case 6 ( Vs = 27 104 and Ri = 1 104 ).. ............ 114 6 9 RSTS profiles for (A) case 5 (Vs = 3 104 and Ri = 1 104 ) and (B) case 6 ( Vs = 27 104 and Ri = 1 104 ) at six different phases.. .............................................................. 115 6 10 TKE budget terms of case 5 ( Vs = 3 104 and Ri = 1 104 ).. ....................................... 116 6 11 Profiles of TKE budget terms for case 2 ( Vs = 9 104 and Ri = 1 104 ) and case 5 (Vs = 3 104 and Ri = 1 104 )at y = 15~25 : .................................................................. 117 6 12 Profiles of TKE budget terms for case 2 ( Vs = 9 104 and Ri = 1 104 ) and case 5 (Vs = 3 104 and Ri = 1 104 )at y = 25~30. ................................................................... 118 6 13 TKE budget terms of case 6.. .............................................................................................. 119 7 1 Illustrative sketch that shows the vertical structure of fluid mud transport. .................... 125 7 2 The location of the simulations in twodimensional parametric space of Vs and Ri ....... 125 7 3 A complete set of simulations in three dimensional parametric space. ............................ 126 PAGE 12 12 Abstract of Dissertation Presented to the Graduate S chool of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy NUMERICAL INVESTIGATION OF PARTICLE LADEN FLOW IN OSCILLATORY CHANNEL AND ITS IMPLICATION TO WAVE INDUCED FINE SEDIMENT TR ANSPORT By Celalettin Emre Ozdemir August 2010 Chair: T. Hsu Cochair: S Balachandar Major: Civil Engineering Studying particle laden oscillatory channel flow constitutes an important step towards understanding practical applications such as sediment transport in coastal environments. This study aims to take a step forward in our understanding of the role of turbulence on fine particle transport in an oscillatory channel and the back effect of particles on turbulence modulation using an Eulerian Euler ian framework. The governing equations are solved by a 3D pseudo spectral Navier Stokes solver previously used for direct numerical simulation (DNS) of turbulent flows. The EulerianEulerian modeling framework chosen in this study, makes the two way coupl ed problem computationally feasible without compromising accuracy as long as particles are of small response time (Shotorban & Balachandar, 2006). As a first step, the instantaneous particle velocity is calculated as the superposition of local fluid veloci ty and particle settling velocity with higher order particle inertia effect is neglected. Correspondingly, the only modulation of carrier flow is due to particle induced density stratification. In the present investigation at moderately energetic condition s the volume averaged concentration and particle settling velocity are varied over a realistic range. The simulation results reveal critical processes PAGE 13 13 due to different degrees of particle turbulence interaction. Essentially, four different regimes of part icle transport for the given Reynolds number are observed within the range of settling velocity ( Vs), which represents the par ticle size, and volume averaged concentration : i) virtually no turbulence modulation in the case of very dilute condition, ii) sl ightly modified regime where slight turbulence attenuation is observed near the top of oscillatory boundary layer. However, in this regime significant change can be observed in the concentration profile with the formation of a sharp gradient in ensemble av eraged concentration profile, lutocline; iii) regime where flow laminarization occurs during peak flow, followed by shear instability during flow reversal. Significant reduction of oscillatory boundary layer thickness is also observed; iv) complete lamina rization due to strong particle induced stable density stratification. These regimes revealed by the numerical simulations are consistent with the limited field observations. However, some of the flow features, especially regarding regimes iii and regime i v,cannot be predicted by the conventional Reynolds averaged models. It is concluded in this study the turbulence resolving simulation tool developed in this study is useful to understand wave induced fine sediment transport processes. Apart from particle s ettling velocity and particle concentration, the developed model can be further enhanced to include particle inertia effects and inter particle influence on the flow. PAGE 14 14 CHAPTER 1 INTRODUCTION 1.1 Overview Coastal zones have been a cradle for many civilizations and they are the forefront of the modern civilization deve lopment. As a result 44 % of the Worlds population lives near the coastlines (U.N atlas of the oceans), which also demonstrates the vitality of the coastal zones to humanity. There is also increasing concern regarding the importance of sustainability in the coastal zones, i.e. the principle of keeping all activities in balance withou t overexploiting the resources. Due t o its impact on geomorphologic evolution of the coastline, local and global environmental complications, and public health and safety, the fate of terrestrial sediment transport in the coastal environment is an importan t phenomenon to be quantitatively assessed. This is strongly demonstrated in the 40 km2 of land lost annually in Louisiana While directly impacting humans living on the coast, this massive amount of cohesive sediment contains high amount s of nutrient s du e to modern agricultural practices, and has reportedly caused hypoxia in the Gulf of Mexico (Goolsby et al., 1999). In a global perspective, massive storm induced runoff is claim ed to be responsible for the de position of carbon, so called carbon sequestrat ion, as the eroded particles containing carbon are deposited in the continental margin and which are removed, at least temporarily or permenantly, out of the carbon cycle (Goldsmith et al. 2008). All these intriguing examples are closely related to the fat e of terrestrial sediment sta rting from its introduction to the river mouth to the final deposition in the continental margin. The terrestrial sediment conveyed through river plumes undergoes several processes until their final deposition at the continental margin ( Wright & Nittrouer 1995; see also Figure 1 1). Depending on the sediment content in the river discharge, sediment particles may be deposited via episodic hyperpycnal plume s (Bates 1953; Mulder & Syvitski 1995) or settle through PAGE 15 15 floccula tion in hypopycnal plume s (e.g. Hill, Milligan, & Geyer 2000) These p articles are then resuspended on the continental shelf through waves, tides and currents. Resuspension processes are very complex to quantify due to complicated fluid sediment interactio n, wave current boundary layer processes (Ross & Mehta 1988; Traykovski et al. 2000) and biological effects (biological packing of fine grained sediment) Hence, fine sediment transport suspension due to waves is the focus of the current study. The continental shelf is comparatively shallow compared to the deep ocean and therefore the bottom boundary layer is dynamic and susceptible to the effects of waves, currents and tides. These effects are of very different time and spatial scales (Larson & Krauss 1995) (Figure 1 2). These variations in scale can range from seconds and centimeters for boundary layer proc esses and up to decades and hundreds of kilometers for st orm events and sea level rise. Failure to understand transport at a specific scale o f transport could lead to orders of magnitude of error in calculating other scales of transport The suspension of particles in turbulent boundary layer involves processes that take place at the smallest time scale but can have considerable long term effe cts ( Mehta 1987; Best 2005; Winterwerp 2001). Unsteady turbulent boundary layer physics itself and turbulence modulation due to particles constitute a major challenge to understand the wave driven fine sediment transport as both these aspects have not been fully understood. Transport and mixing due to current s wave s and tide s are very likely to be in effect concurrently. Though the effect of waves on suspension of fluidmud is known as early as Alishahi & Crone (1964), Jackson (1973), its significance is r ecognized by Wells & Coleman (1981), Maa & Mehta (1987), and Sternberg et al. (1996). Initial attempts to evaluate these transport mechanisms with large scale models failed to predict sediment transport quantitatively. While qualitatively good agreement is obtained by PAGE 16 16 including a wave boundary layer (WBL) module ( e.g. Scully 2003) model results are often sensitive to small changes in the input parameters (Harris Traykovski, & Geyer 2005). Large scale fiel d observations can reveal critical processes that ar e missing in our current understandi ng of the problem. F ield experiments of sediment transport are under the influence of multiple effects with further uncertainties such as measurement errors in turbulence, particle size and data processing technique s for the measurement devices. This study aims to accurately describe the smallest scale of turbulence in the wave boundary layer in an effort to understand the physics behind large scale sediment transport. The aspects that have been emphasized heretofore have also been pointed out in the studies of Hsu et al. (2004) and Hsu & Hanes (2004). In these two studies wave resolving two phase models with two equat ion turbulence closure models a re proposed to better understand the medium to coarse sediment suspension and deposition in the WBL. The two phase model in these two studies is based on EulerianEulerian framework incorporating particle particle interaction for coarse particles via kin etic theory of granular flow. For fine sediment with small particle respons e time the problem gets more complicated and the available closure for particle particle interaction can only be incorporated through empirical relations primarily due to flocculati on. In the later study of Hsu, Traykovski, & Kineke (2007) and Hsu, Ozdemi r, & Traykovski (2009) particle particle interaction has also been incorporated highly viscous rheological properties of the sediment fluid m ixture (Hsu et al. 2009) and its impact on the fluid mud transport is investigated. The details of these two studies are discussed in section 1.2.1. The turbulence models incorporated in both studies are two equation k closure models. Standard k closure model uses empirical parameters in the transport of the turbulent kinetic energy (TKE) equation and are based on the isotropic turbulence assumption. The transport terms of pressure PAGE 17 17 and turbulent convection are modeled by eddy diffusivity assumption, i.e. they are assumed as diffusive processes. The dissipation rate trans port equatio n is evaluated empirically. Moreover, in these studies the loss due to particle density stratification is again modeled by eddy diffusivity assumption. Although the aforementioned studies provide a more detailed intra wave picture for coastal sediment transport applications, a turbulence resolving approach is required to fully understand fine sediment dynamics near bottom boundary and to verify findings obtained from Reynolds averaged approach. The advantage of resolving the turbulence at all t he scales i.e, from the largest scale to dissipative Kolmogorov scale will yield the response of turbulence to the particle density stratification due to settling without empirical closure. Of relevance is the recent work of Cantero et al. (2009) who study fine sediment transport at steady state in a fully turbul ent turbidity current. They use a turbidity current with a roof model to consider the statistically stationary state. Their simulation results reveal a decrease in turbulent intensity with increasin g stratification and the variation of the effective von Krmn constant was revisited. Most importantly, it was observed that at a specified dilute particle concentration (correspondingly a fixed Richardson number) when particle settling velocity increased above a threshold the stratification effect became strong enough to kill turbulence and the flow becomes laminarized. The effect of turbulence on the suspension of fine particles and the back effect of stable density stratification on turbulence in an unsteady oscillatory cha nnel is not well understood This study aims to investigate th e influence of particle induced density stratification on the flow turbulence which dictates suspension and deposition of particles in the WBL for moderate field conditions with different availability fine sediment. While the turbulence is fully resolved at all the scales, the back effect of fine particles are taken into account in a Eulerian Eulerian PAGE 18 18 framework that makes the problem computationally feasible without compromi sing the accuracy (Ferry & Balachandar 2001). As a result of this study, the impact of particle density stratification on flow turbulence manifest s itself in four different flow regimes observed with varying volume averaged concentration. i) virtually no t urbulence modulation in the case of very dilute condition, i.e., ~ 0; ii) slightly modified regime where slight turbulence attenuation is observed near the top of oscillatory boundary layer. However, in this regime significant change can be observed in the concentration profile with the formation of a lutocline; iii) regime where flow laminarization occurs during peak flow, followed by shear instability during flow reversal. Significant reduction of oscillatory boundary layer thickness is also observed; iv) complete laminarization due to strong particle induced stable density stratification. These results are in coherence with many field and laboratory observations (See Winterverp & Van Kesteren 2004 and the references therein). However the current study pr ovide s a better understanding of the particle laden WBL processes and benchmark results that can be further used to test models with turbulence closure. This study will also be useful to further model the fate of terrestrial sediment in the continental she lf, and enhance large scale modeling studies by providing a tool for parameterizing these small scale processes. 1.2 Literature Survey Suspension and deposition i n the WBL are depende nt upon the unsteady turbulence as well as the back effect of the particl es on the flow turbulence. Unfortunately, there is no turbulence resolving study which considers both effects concurrently. Either particles are treated as passive tracers in oscillatory flow or the effect of particles is studied in steady channel. Theref ore, the relevant literature reviewed here falls into two categories: two phase turbulence modulation in simple flow condition and single phase turbulent flow characteristics of oscillatory flow. PAGE 19 19 1.2.1 Turbulence Modulation due to Fine Particles Despite pi oneering studies of turbulence modulation by experiments (Rogers & Eaton 1991, Tanaka & Eaton 2008), point particle simulations (Maxey, 1987; McLaughlin, 1989; Squires & Eaton, 1991; Elghobashi, 1991; Elghobashi & Trusdell, 1992; Elghobashi & Truesdell, 19 93), and fullyresolved simulations (Bagchi & Balachandar 2003, 2004, Burton & Eaton 2005, Zeng et al. 2008, Zeng, Balachandar, & Najjar 2009), our understanding on twoway coupled turbulence modulation remains incomplete. In general, small neutrally buoya nt particles at low concentrations tend to suppress turbulence, while larger particles contribute to turbulence enhancement. There are several mechanisms that contribute to modulation of carrier phase turbulence. Added inertia of the suspended particles, i ncreased dissipation due to hydrodynamic drag on the particles and enhanced effective viscosity of the suspension are three mechanisms that contribute to suppression of carrier phase turbulence. On the other hand, if the particle Reynolds number is suffici ently large, carrier phase turbulence can be augmented through wake oscillation and vortex shedding. Finally, local vertical variations in particle concentration can introduce stable or unstable stratification and contribute to turbulence modulation (see Balachandar & Eaton 2010). Here the scope of this study is limited to a sufficiently dilute transport of fine particles that the particle inertial, enhanced viscosity and vortex shedding effects are unimportant. Our focus is on the last mechanism where the effect of particle induced density stratification on carrier flow turbulen ce is of critical interest. In the case of fine particles (or particle flocs) in water their re lative negative buoyancy results in settling towards the bottom surface. Turbulence near the bottom boundary will contribute to the resuspension of particles and an effective diffusion of particles away from the wall. In a statistically stationary flow, as in a turbidity current, the above two mechanisms will statistically balance and determine the vertical variation of mean particle concentration and the PAGE 20 20 degree of stable density stratification. In the context of wave induced oscillatory turbulence, the above balance is still active and resu lts in a stably stratified particle concentration. The earliest study which addressed this problem was carried out by Vanoni (1946), which was followed by Einstein & Chien (1955). Both of these experimental works concluded that the change in Reynolds avera ged mean velocity is due to suspended particles and the modified mean velocity profile was parameterized in terms of a reduced value of von Krmn constant. More importantly, Vanoni (1946) hypothesized that stratification is responsible for the decrease i n turbulence, which in turn contributed to the modification of the mean flow velocity. Later studies account for the effect of particleinduced stable density stratification using different parameterizations. Gelfenbaum & Smith (1986) proposed a reduction factor for eddy viscosity as a function of Richardson number. Hino (1963) and Zhou & Ni (1995) evaluated changes to the von Krmn constant in terms of particle concentration. Turner (1973) and Barenblatt (1953) adopted damping functions for eddy viscosity as a function of MoninObukhov length scale (for further details refer to Winterwerp 2001 and the references therein). The numerical study by Noh & Fernando (1991) utilized one equation closure turbulence model incorporating the effect of particles on tur bulent flow due to density stratification to model particle suspension by an oscillating grid and the formation of lutocline. Winterwerp (2001, 2002) and Hsu et al. (2007, 2009) have used twoequation turbulence closure models that further incorporate the back effect of particles on the turbulent flow. 1.2.2 Boundary Layer Turbulence in Oscillatory Flow The unsteady nature of oscillatory particulate flow presents an important challenge. Under ideal sinusoidal forcing, the flow velocity varies between positi ve and negative peaks with flow reversal in between. If the instantaneous flow velocity becomes sufficiently large, then perturbations in the flow start to grow during the acceleration phase and depending upon the PAGE 21 21 extent of the acceleration phase the flow might become fully turbulent. During deceleration the perturbations might fully decay and the flow might relaminarize. Or, depending on the maximum velocity that the flow reaches, the perturbations might only partially decay, and the flow might remain full y turbulent over the entire oscillatory cycle. In flow of this type, the Reynolds Stokes strongly controls the physical features of the flow. The unsteady nature of the flow motivated early researches to identify the onset of turbulence theoretically (Kerc zek & Davis, 1974; Hall, 1978; Akhavan, Kamm, & Shapiro 1991b) experimentally (Akhavan Kamm, & Shapiro 1991a; Sarp kaya, 1993) or both (Hino, Sawamoto, & Takasu 1976) in these kinds of complex flows. In all these studies, the Stokes Reynolds number is sel ected : 0RefU (1 1) w here 0U is the dimensional maximum free stream velocity, is the dimensional Stokes boundary layer thickness defined as 2/f being the dimensional angular frequency of oscillatory forcing, and f is the kinematic viscosity of the fluid. Although the boundary values of Re remain uncertain, these studies have classified four flow regimes based on the Stokes Reynolds number: (1) laminar regime (Re < 100 ); (2) disturbed laminar regime (100< Re < 800; Hino et al., 1983); (3) intermittently turbulent regime (800< Re < ~1200; Hino et al. 1983; Sarpkaya, 1993); and (4) fully turbulent (Re ~1 200; Sarpkaya 1993; Jensen, Sumer, & Fredse 1989). In the laminar regime, the flow stays laminar over the entire wave cycle and in the fully turbulent regime it stays turbulent throughout the wave cycle. In the perturbed laminar regime, t he flow is perturbed during acceleration phase and laminarizes during deceleration without ever becoming fully turbulent. Finally, in the intermittently turbulent regime the flow becomes turbulent during acceleration and tends to decay during deceleration without complete PAGE 22 22 laminarization. Later studies by Jensen et al. (1989) and Chen et al. (2007) are at higher Re and they present detailed turbulent quantities as well as discuss evolution of flow structures. Direct numerical simulations, the selected mean s to solve the flow field in this study, of oscillatory flow in a clear fluid started with the study of Spalar t & Baldwin (1989). They carried out simulations in the intermittently turbulent regime for Re = 800, 1000, and 1200 using highly accurate pseudo spectral method (Canuto et al. 1987) Their main objective was to provide comprehensive turbulent flow statistics. To identify the nature of instability, Akhavan et al. (1991b) simulated oscillatory pipe flow in the intermittently turbulent regime. To bet ter understand the phenomenon of onset of instability and evolution of vortex structures, Vittori & Verzi cco (1998) and Costamagna, Vittori, & Blondeaux (2003) carried out simulations using a secondorder finite difference scheme. Salon, Armenio, & Crise (2007) have performed large eddy simulations of the oscillatory turbulent flow and tested subgrid scale LES closure models by using the DNS data provided by Costamagna et al. (2003). Recently, Ozdemir, Hsu, & Balachandar (2009) performed numerical simulati ons in the intermittently turbulent regime for very dilute particle concentration, where particles can be assumed to be passive and not influence the carrier flow. They adopt a highly accurate pseudospectral flow solver and implement the same domain size and resolution as Spalart & Baldwin (1988), whose simulations were of the highest fidelity since they employed the largest computational domain, the finest resolution and spectrally accurate numerical methodology. Ozdemir et al. (2009) simulated the case R e = 1000 for very dilute concentration where the particle effect can be ignored. They also evaluated the suspension and settling of fine particles both qualitatively and quantitatively. However, all these simulations ignored the back particle effect which is the focus of this study. PAGE 23 23 1.3 Scope of the Thesis This thesis is organized as 7 chapters. C hapter 1 presents the motivation together with the review of the relevant literature. In C hapter 2 the statement of the problem and the proposed methodology to s olve the stated problem are given. The numerical formu lation as well as methodology are presented in Chapter 3. Chapter 4 presents the very dilute particulate flow where the effect of particle density stratification can be ignored. The same chapter can als o be regarded as validation of the numerical methodology and formulation as the turbulent quantities are validated with the available numerical simulations of Spalart & Baldwin (1989). In C hapter 5 the effect of stable particle density stratification at d ifferent v olume averaged concentration is presented with detailed analyses. The effect of particle settling velocity is evaluated in Chapter 6. Finally, C hapter 7 summarizes the findings related to the current study together with recommendations for future studies. Figure 1 1. Illustrative sketch for the processes involved in sediment source to sink. PAGE 24 24 Figure 1 2 Temporal and spatial scales for different sediment transport processes (Adopted from Larson and Kraus 1995) PAGE 25 25 CHAPTER 2 PROBLEM STATEMENT AN D MATHEMATICAL FORMULATION 2.1 Problem Statement In this study, sediment transport in the WBL problem is idealized as particulate flow in an oscillatory channel. The flow is statistically homogeneous in streamwise and spanwise directions, i.e., 1DV. If the flow becomes nonturbulent without particles, the problem is simplified into Stokes second problem (Figure 21). As d iscussed in the C hapter 1 the WBL keeps the particles suspended. If there is a small gravitational force, considerable sediment transport might occur as a result. On the other hand, the turbulent characteristics of the flow keeps the fine particles suspen ded which is concurrently under the influence of a back particle effect in the form of particle density stratification. Hence, the two way coupled turbulence particle interaction and its effect on particle suspension and deposition events is the main focu s of this study. Therefore, in the context of understanding suspension and deposition events in the WBL which is directly affected by modified or unmodified turbulence, the stated problem is sufficient to mimic the actual conditions in a moderately energe tic field on the continental shelf. Although the mathematical formulation for fine sediment transport in the WBL is idealized, the range of flow parameters is realistic if compared to that of field observations (e.g., Traykovski et al. 2000; Traykovski, Wi eberg, & Geyer 2007). Although a wide range of dimensional variables satisfy the nondimensional parameters that are discussed in Section 2.2, the following can be given as an example for comparison with the realistic field conditions. I n this study, surfa ce waves at the continental shelf of 10 sec wave period are considered and the corresponding Stokes boundary layer thickness is 1.8 mm. Here, Reynolds number using Stokes boundary layer thick ness as length scale ( Eq. 1 1 ) is 1000, which corresponds to an oscillatory flow velocity amplitude of 0.56 m/s. This value PAGE 26 26 falls within the range of velocity magnitude that observed at continental shelf where fluid mud is suspended in the wave boundary layer (e.g., Traykovski et al. 2000; 2007) A suspension of fine specific gravity of 2.65, which results in a still fluid settling velocity of about 0.5 mm/s. The particle Reynolds number, Re()/psfVd becomes 0.02 which is sufficiently smaller than unity. In the definition of particle Reynolds number given d stands for the particle diameter, sV is the dimensional settling velocity of the particle. Therefore, t he corre sponding particle time scale to the carrier fluid motion defined as 2/18pfd is 8.5105 s, where p is particle densi ty, d is particle diameter and f is dynamic viscosity of water. 2.2 Mathematical Formulation For the current case of Re1000 the flow is turbulent provided the back effect from the suspended particles is small. The flow thus presents a range of time scales from the integral to the Kolmogorov scale. The time scale of the Stokes layer ( 0/ U ) is more than an order of magnitude larger than the particle time scale. From the numerical results we have verified that the particle time scale is also smaller than the Kolmogorov time scale. It can be established that the particles are sufficiently small that their Stokes number, defined as the ratio of particle to fluid time scale, is less than 1. As a result, here the inertial effect of particles is ignored (Ferry & Balachandar, 2001) the nond imensional particle velocity P iU is defined as the sum of fluid velocity f iU and particle settling velocity, SV as: 2 pf iisiUUV (2 1) where superscripts f and s stand for the fluid and particle phases, respectively. The purpose of this simplification is the reduction in the number of partial differential equations to be solved, PAGE 27 27 which results in computational efficiency. The accuracy of t his simplification has been tested in the study of Shotorban & Balachandar (2006) with the LagrangianEulerian counterpart. The assumption is valid as long as the ratio of the particle time scale to the time scale of the flow is less than unity, and very c lose results are obtained. By substituting the above algebraic relationship in the EulerianEulerian two phase equations and applying the Boussinesq approximation, valid for relatively small concentrations, the resulting governing equations can be greatly simplified (e.g., Cantero, Balachandar, & Garcia 2008). The continuity equation ( Eq. 2 2 ) an d momentum equations ( Eq. 2 3 ) of the fluid phase can be written as 0f U (2 2) '2 1222 1 sin() () ReRe Re U eeUf fD tPRiC Dt (2 3) These equations are nondimensionalized by length, velocity, and time scales selected to be the Stokes boundary layer thickness ( ), maximum free stream velocity ( 0U ), and 0/ U respectively In the present problem of an oscillatory channel flow, sufficiently away from the boundary layer, the time variation of the dimensional free stream velocity is given as follow 01cos() Ue Ut (2 4) w here 1e and 2e correspond to unit vectors along the streamwise and wall normal directions, respectively. Far from the bottom boundary, well outside the oscillatory boundary layer, the momentum equation can be simplified for inviscid fl ow to yield the mean streamwise pressure gradient to be 0sin()fPUt x (2 5) PAGE 28 28 Eq. 2 5 upon non dimensionalization, along with the definition of Stokes layer thickness and Reynolds number, reduces to the first term on t he right hand side of Eq. 2 3 The second term on t he right hand side of Eq. 2 3 corresponds to the fluctuating pressure gradient. The third right hand side term in Eq. 2 3 represents the back effect of the particle on the fluid phase via particle induced stable density stratification, which is only effective in the vertical direction. Here the bulk Richardson number Ri is defined as: 2 0()pf ave fg Ri C U (2 6) in which, f is fluid density and C is the volume averaged concentration over the entire computational domain, i.e. (,,)aveCxyzd C (2 7) Here, is the volume of the domain. In the present formulation, C is normalized by aveC Consistent with Eq. 2 1 in the present formulation the momentum equation does not include frictional or drag forces arising from the inertial effect of particles. Hence, preferential accumulation of particles is not expected to be observed in our simulations, which is appropriate for particles of very small response time. Finally, the particle concentration is calculated by the conservation of particle mass: 2 21 ()[()] Ref sCVCC t Sc Ue (2 8) The term on the right hand side of Eq. 2 8 represents an effective diffusive flux. Here, the Schmidt number is /fSc where is the effective diffusivity of the particles. The governing equation of particle moti on, when applied in the Lagrangian framework, contains no PAGE 29 29 diffusive effect, i.e. Sc However, in the present simulations, a diffusive term is required not only for numerical stability but also to account for subscale particle motion through sub grid interactions. Birman, Martin, & Meiburg (2005) show that the effect of Sc within the range of (0.2, 5.0) to be quite insensitive in moderately energetic gravit y flows. Similarly Bonometti & Balachandar (2008) have investigated the effect of Sc in the range of ) ,[ 1 for gravity currents with similar governing equations and observed Schmidt number effect to be small in high Reynolds number flows. It is observed that in vigorous flows, mixing is dominated by advection and the effect of molecular mixing represented by diffusion to be quite small. Sc is selected to be 0.5 in the present simulations. The velocity boundary conditions at both the top and bottom boundaries are noslip wall boundary conditions. Periodic boundary conditions are imposed at the lateral boundaries of the computational domain (see F igure 2 1). For particle concentration, the following relation is imposed at both top and bottom boundaries of the computational domain: 2 21 ReSCC V yScy (2 9) This boundary condition imposes no net transport of particles across the top and bottom boundaries. Periodic boundary conditions for particle concentration are applied along the streamwise and spanwise directions. These bounda ry conditions in turn imply that the total volume of particles remains constant in the domain throughout the computation. This integral property has been confirmed as a verification test in the present computations. In the numerical description given here, the effective parameters for the particle flow are the settling velocity, sV and the bulk Richardson number, Ri In the fluid momentum equation, the back particle effect is taken into account by the term, 2RiC e The effect of sV is PAGE 30 30 incorporated through advection diffusion equation of concentration, C as follows: the higher the sV the more the particles settle and become concentrated near the bed and modulate the boundary layer turbulence. While sV discerns the effect due to particles of different size, Ri reflects the effect of volu me averaged concentration (see Eq. 2 6 and Eq. 2 7 ). For non cohesive sediment size of particle s and flow energy determine the erosion capacity of the flow from the bed. The interaction between the bed and the carrier flow leads to quasi constant concentration, the equilibrium concentration (see for example Van Ri jn, 1981). For fine sediment, due to the small settling velocity turbulence modulation owing to particle phase (fluid mud) becomes important. According to Winterverp & Van Kesterens (2004) description for steady and tidal flow, a two layered fluid is formed and at the interface of these two layers there emerges a huge decay of turbulence causing a total collapse in turbulence. In this perspective, volume averaged concentration, correspondingly Ri and evaluating the flow turbulence within a range of Ri will help to understand the concept of saturation concentration (Winterverp & Van Kesteren, 2004) for WBL. Therefore, in this study the first part of the proposed simulati ons set 1, covers a range of Ri (see Table 2 1 ), while making a wash load assumpti on at the bottom boundary (see Eq. 2 9 ). The remaining simulations, set 2, are employed to investigate the effect of particle settling velocity, Vs, on the carrier flow turbulence and particle suspension and mixing characteristics. PAGE 31 31 Figure 2 1. Two dimensional descriptive sketch of the problem The velocity profile is the laminar veloci ty profile at t Table 2 1. List of simulations performed in this study The first column shows the set number. The first set of simulations aims to further identify the effect of volume averaged concentration and the second set of simulations aims to identify the effect of p article settl ing velocity. Case 2 shown with asterisk is used as the reference case for evaluation of both particle settling velocity and volume averaged concentration. T he second colum n shows the number of the case, the third column shows the Ri spec ified for that sim ulation, fourth column s hows the settling velocity. Fifth and sixth column shows t he corresponding Schmidt number and Reynolds number Set number Case number Ri V s Sc Re 1 1 0 9 10 4 0.5 1000 2* 1 10 4 9 10 4 0.5 1000 3 3 10 4 9 10 4 0.5 1000 4 6 10 4 9 10 4 0.5 1000 2 5 1 10 4 3 10 4 0.5 1000 6 1 10 4 27 10 4 0.5 1000 7 3 10 4 27 10 4 0.5 1000 PAGE 32 32 CHAPTER 3 NUMERICAL METHODOLOGY A three dimensional pseudo spectral (PS) flow solver is employed to solve particle laden oscillatory turbulent channel flow. Pseudospectral methods (PSM) have comp licated mathematical background. I n C hapter 3 the ma in objective is to describe the main components of the PSM for the current set of two phase flow equations More details on the numerical methodology can be found in Canuto et al. (1987) and Peyret (2002). Section 3.1 is an overview of the domain and grid setup. In Section 3.2, general inform ation on PSM is given. Section 3.3 describes the implementation of PSM in the stated problem of this study. In Section 3.4 the implementation of Helmholtz and the Poisson solvers are described. 3.1 Overview The size of the computational domain is selected to be the 60 60 30 along the streamwise, vertical and spanwise directions, respectively. Thus, the distance from the bottom boundary to the center plane of the channel is 30 times the Stokes layer thickness. The selected domain size is exactly the same as that of Spalart & Baldwin (198 9 ) and their domain is the largest and the most conservative among the different numerical studies. As it will be shown in Chapter 4 the two poin t correl ation functions along the streamwise and spanwise directions suggest that this domain size is sufficiently large to capture the main turbulent characteristics at Re = 1000 in the case of clear fluid. The adequacy of spatial resolution is examined in term s of wall units. Friction velocity in dimensional terms is defined as : w fu (3 1) w here w is the wall shear stress. In non dimensional terms the friction velocity becomes PAGE 33 33 1 Rewu u y (3 2) where w represents an average over the bottom wall. The time history of non dimensional friction velocity and the corresponding centerline velocity are shown in Figure 3 1 for the case s of Ri = 0.0, 104 and 3104. Also marked in Figure 3 1 is the definition of phase ( ) to be adopted in this study. According to this definition, = 0 corresponds to the phase of maximum positive streamwise velocity, corresponds to maximum negative velocity, and /2 correspond to flow reversals. When Ri =0.0 or 104 the flow near the bottom wall remains turbulent and the friction velocity (correspondingly wall shear stress) peaks at about the same time as velocity maximu m. For the case of Ri = 304, as discussed in C hapter 4, the flow near the bottom wall is laminarized and as a result the peak friction velocity leads the velocity maximum by /4 in agreement with the expected behavior of lamina r oscillatory Stokes layer. The non dimensional coordinate in wall units can then be expressed as : ,maxRe u xx (3 3) w here ,maxu is the maximum non dimensional friction velocity at the bottom wall. The r esults to be reported here were performed with a grid of 192193192 (~ 7.1 million) points along the streamwise, wall normal and spanwise directions. For the case of clear fluid ( Ri = 0.0) the + = + = 7.8. Along the wall normal direction, the resolution y+ = y+ = 24.5. For the other three cases where particles attenuate turbulence, the peak friction velocity is lower and as a result the grid resol ution is higher The grid resolution used here is fine enough to capture the entire range of scales, which can be verified from the power spectrum of the velocity fluctuation. The adequacy of the present grid for the clear fluid has been established by Spa lart & Baldwin (1989) PAGE 34 34 and by Ozdemir et al. (2009). The detailed turbulent energy spectra shown for 0 Ri in C hapter 4 and for Ri = 104 and 3104 in C hapter 5 confirm the adequacy of resolution. 3.2 Pseudo S pectral Methods The main advantage of using PSM is that the decay rate of numerical error is exponential rather than polynomial. In the case of finite difference (FD) approach, the decay rate is polynomial as opposed to exponential. Hence much higher spatial resolution is required for FD scheme to achieve the same level of accuracy obtained in PSM. P S representation of a function f at a single point xi is given as: 0()()N illi lfxax (3 4) In Eq. (3 4) ()ifx is the function to be evaluated at point ix where 1 iN and N is the number of discrete points. l is the orthogonal base function and la is the associated coefficient. Function f is the sum of ()lliax which is a high order polynomial sum of all components Though it seems unnecessary to evaluate the local value of a derivative of a function by usi ng all the discrete points in the domain, the exponential decay of error in PSM with the increase in the number of points make it possible to achieve high accuracy with fewer points unlike its FD and finite volume counterparts. On the other hand, PS M are n ot very suitable for irregular geometries or functions with very sharp gradients such as very strong shock waves, because t hey use smooth base functions for the polynomial sum ( Eq. 3 4). In the current study, these restrictions are no t problematic because the geometry is simple and the governing equations conta in diffusive terms which smooth en the flow field. As mentioned in Chapter 2, boundary conditions along the streamwise and spanwise directions are peri odic. Due to periodic boundary conditions, Fourier series are selected as base PAGE 35 35 functions in these directions. Also the use of the Fast Fourier Transform (FFT) algorithm w ill ease the computational cost Comparing the computational costs shows that the cost of FFT for N grid points is on the order of O ( log NN ), whereas wit the use of matrix multiplication for derivative evaluation the computational cost would be O ( 2N ). The boundary conditions at the top and the bottom of the domain are no sl ip and nopenetration Therefore, Chebyshev polynomials are selected as base functions which are well suited to non periodic boundary conditions in bounded domains. ChebyshevGaussLabatto quadrature point s are selected. Chebyshev GaussLabatto quadrat ure points are defined as follows: cos()jjy (3 5) (1)jj (3 6) N (3 7) where 11 jN In other words, these points are the vertical component of the points around a unit circle with equispaced angles of value In this configuration the points are clustered at the top and the bottom of the domain. This gives the advantage of higher resolution near the wall boundaries where boundary layer turbulence plays an important role. T he increa sed resolution near the boundaries allows the capture of the dissipative Kolmogorov scale It should also be noted that selection of Chebyshev G aussLabatto collocation points include the points at 1 y making it simpler for the implementation of a combination of Drichlet and Neumann boundary conditions at the top and bottom (refer to Canuto et al. 1987 and the references therein). To avoid confusion it should also be mentioned that the domain size in the vertical direction can be PAGE 36 36 manipulated with proper transformation. In our simulations it has been transformed with the following relation: *(1) yy 30 (3 8) As stated in Eq. 3 1 the variables to be solved can be represented in 3D as : 11 22 0 2 2 (,,,) (,,,)()x z y xizi xz z xN N N ixiz ijk xyz jj NN juxyzt uteeTy (3 9) In Eq. 3 9 u is the coefficient of the triple expansion and is a function of wave numbers i.e x j z in all three directions of the domain. ()jjTy is the Chebyshev polynomial of order j This representation can b e applied to variables of interest H owever for the sake of brevity only streamwise velocity compo nent is shown here. The c ollocation points can be presented as follows: (1) (,,)(,cos(),) 1 xx z xyziL kL j ijk NNN (3 10) First order derivatives in the periodic x and z directions can be evaluated in spectral space simply as the multiplication of each term with x i and z i respectively. Similarly, the seco nd order derivatives in these direction s are simply the Fourier transform scaled by 2 x and 2 z These are summarized in Eq. 3 11 to Eq 3 14 for the streamwise velocity as : () xu iu x (3 11) () zu iu z (3 12) PAGE 37 37 2 2 2() xu u x (3 13) 2 2 2() zu u z (3 14) The evaluation of the derivatives is performed in the spectral space for the periodic directions. On the other hand, Chebyshev derivatives, first and second order, c an be evaluated in real space by simply multiplying by derivative matrices. The first order derivative matrix elements can be given as follows (Canuto et al. 1987) : (1)(1) () 2sin[()/2]sin[()/2]jl j jl la ajlNjlN D for jl (3 15) (1) 2() 2sin() Dj jly j N for 11 jlN (3 16) 2 (1)21 () 6 DjlN for 0 jl (3 17) 2 (1)21 () 6 DjlN for jlN (3 18) Higher order derivatives are simply: ()(1) 1DDn n p (3 19) where D is the Chebyshev derivative matrix and superscript (n) stands for the order of derivative. At this point it should be mentioned that the coefficients of D (n) are dependant upon the location of the collocation points. The main reason behind this is the dependence of the Chebyshev transform coefficients i.e. al in Eq. 3 1 on the location of the collocation points, i.e. PAGE 38 38 y*. This is also the main reason that makes it possible to evaluate the derivative in the real space (for detailed derivation refer to Peyret 2002) 3.3 Solution of Governing Equations The velocity, s tated in Chapter 2, is solved in four steps using a three stage integ ration scheme given thro ugh Eq 3 20 to Eq 3 28: () mHt (1) (1) ()1()mmAcnlmH U (3 20) ()(1) 2()mmcnlm UU *()(1)()() ()[()()] mmmmHcdmDDRiC UU (3 21) 2() ()1 2()mmp cdm U (3 22) ()() ()mmcdm UU () mp (3 23) (1)(3) ()nt UU (3 24) In Eq. 3 20 to 3 24, () A stands for the advective terms () D () () and 2() are diffusion, gradient, divergence, and Laplace operators The double cap on the top of these operators and the other terms shows that they are evaluated in spectral space in x and z direction () mU and () mH are intermediate variables evaluated in the spectral space T he coefficients are as follows: (1,2,3) m (3 25) 5153 1()(0,,) 9128 cnlm (3 26) 1158 2()(,,) 31615 cnlm (3 27) 5 ()(,,) 6248 ttt cdm (3 28) PAGE 39 39 1 cnl is the coefficient of the nonlinear terms at the first step, and 2 cnl is the coefficient of the non linear terms in the second step. cd stands for the coefficients of the diffusive terms. As can be noted f rom Eq. 3 20 to Eq. 3 24, this sc heme has low storage properties while keeping the accuracy at third order, i.e O ( 3t ) T he only variables that require storage are U and H This is advantageous for DNS simulations when the number of grid points increase s to satisfy the desired grid resolution. Another observation is that the diffusive terms are treated implicitly with a Crank Nicolson integration scheme. Therefore, this scheme can be called a mixed implicit/explicit time discretization (Canuto et al. 1987). All these calculations are performed in Fourier space except in the fi rst step where non linear advection terms are evaluated first in the real space to avoid complicated convolution sums. After evaluating these terms in real space they are transformed back to the Fourier space. In this step, nonlinear d ifferential operato rs introduce aliasing errors which can be described as the cascade of ener gy beyond the highest wave number bounded by the number of collocation points N. This problem is remedied by using the 3/2 dealiasing rule along the periodic directions and the implementation of Arakawa method (Arakawa 1966). The 3/2 rule is nothing but extrapolation of a va riable of size N to 3 N/2 in Fourier space by setting the higher frequencies to zero. Keeping the number of discrete points as 3 N/2, the variables ar e transformed into real space and non linear terms are calculated and back transformed to Fourier space. These terms of size 3 N/2 are then interpolated and sized back to N. Hence, the da ta leakage to higher wave numbers is minimized. The reason of using fa ctor 3/2 is derived in detail in Canuto et al. (1987) and Peyret (2002). At this point it is necessary to mention that /2 rule only eliminates the aliasing errors in periodic directions. The aliasing errors introduced by the Chebyshev PAGE 40 40 derivatives in the wall normal direction are eliminated by the use of Arakawa method (Arakawa 1966). In the second step intermediate velocity, () mU i s evaluated using advective, () mH and diffusive terms, (1)()mDU and ()()mDU This evaluation is done through the solution of the Helmholtz equation: () ()(1) ()()2()mmmcdmD cnlm UUU *() (1)() ()()mmmHcdmDRiC U (3 29) This intermediate velocity is then corrected at step three and four i.e, Eq. 3 22 and Eq. 3 23. At the third step the pressure field, () mp is evaluated in the spectral space by solving th e Poisson equation in Eq 3 22. Then the intermediate flow velocity is corr ected at the fourth step as in Eq. 3 23. As part of t he three stage Runge Kutta scheme this process is repeated three times to complete one full time step. Periodic boundary conditions are automatically satisfied in streamwise and spanwise directions due to FFT implementation in these directions. No slip and no penetration wall boundary conditions are implemented at the top and bottom wall and are implemented by proper modifications at the first and the last row s of the matrix operators of the Helmholtz solver. Although there is no formal boundary condition s pecified for pressure, the Neumann boundary condition is specified at the walls. The reason behind this is that intermediate velocity is evaluated without the pressure gradient term. To keep the slip velocities at the wall minimal, pressure gradients shoul d be set to zero at the top and bottom. To elaborate, the final velocity at the wall at the end of the fourth step can be given as: ()() () 1 1 0,60 [2()]mmm ycdmp eU=eU (3 30) ()() () 3 3 0,60 [2()]mmm ycdmp eU=eU (3 31) PAGE 41 41 w here 1e and 3e stand for the unit vectors directed along streamwise and spanwise direction s respectively. Tangential pressure gradient at time step m, () mp is unknown. Therefore, approximate value is used. The approximation can be done by using a Taylor series expansion of third order: () (1) (2) 2 11 1() () [1] []() (1) (1) mm mcdm cdm pp p Ot cdm cdm ee e (3 32) Note that selection of third order accuracy is in compliance with the third order accuracy achieved in time integration scheme. The solution procedure for the momentum equation has been given in Eq. 3 20 to Eq. 3 24. T he advection diffusion equation of concentration is concurrent ly evaluated with the momentum equation, i.e. Eq. 3 20. The equation of particle conservation is similarly evaluated as the momentum equation as in the following steps: () mHt (1) (1) ()1()mmACcnlmH (3 33) ()(1) 2()mmCCcnlm () (1)() ()[()()] mmmHcdmDCDC (3 34) (1)(3) ()nCtC (3 35) It is noted that only the Helmholtz equation should be solved in Eq. 3 34. The coefficients are again the same as shown in Eq. 3 26 to Eq. 3 28. It is critical that Eq. 3 30 is solved before Eq. 3 21 and the resultant C value is used in the evaluation of Eq. 3 21. 3.4 Solution of Helmholtz and Poisson Equations In Section 3.3 the general methodology of the solution p rocedure for the fluid and particle phases are given. Two major linear algebraic operations involved in solving the governin g equations are the FFT and the matrix operations. As they are widely used in advanced PAGE 42 42 comput ing both these operations are well documented in literature. However, it is worth discussing the solution procedure of Helmholtz/Poisson equations. To get the pre ssure, velocity and concentration a matrix inverse operation should be performed. Eq. 321 suggests the coefficient matrix is dependant on wave numbers, x and z The implication of this is that the inverse operation has to be performed for each combination of wave numbers and stored, which consequently increase s the computational cost. Simply diagonalizing the (2)D matrix eliminates the necessity to store the inverse of the coefficient matrix and yields an efficient solver. Eq 3 21 can also be written as follows: (2)() 22()()Re () () Dm mm xz xRHS cdm UU (3 36) Here, (2)D is represented as 1TT where T is the modal matrix of (2)D and is the diagonal matrix made up of the eigen values of (2)D Then Eq. 3 36 becomes: 1() 22()()Re () ()m mm xz xTT RHS cdm UU (3 37) Eq. 3 37 is then manipulated into the form: () 1()1 ()mm xTTRHS I U (3 38) where 22Re ()xzcdm (3 39) With the implementation of Eq. 3 38, one time calculation of T its inverse, 1T and is sufficient. This reduces the computational cost. For illustration purposes this is only shown for Eq. 3 21 though it can be extended to Eq. (3 22) and Eq. (3 3 4 ). PAGE 43 43 Figure 3 1. Time series of friction velocity for (A) Ri = 0.0, (B) Ri = 4110 and (C) Ri = 4310 PAGE 44 44 CHAPTER 4 CLEAR FLUID SIMULATI ONS: PASSIVE CASE The clear fluid simul ation results presented in C hapter 4 reveal the sediment transport mechanisms in the simplest scenario where settling sediment particles are completely passive to the carrier fluid. According to Noh & Fernando (1991) when the mass concentration of sed iment particles is on the order of O (1) g/l particles remain passive to the carrier fluid. Also clear fluid simulations serve as a validation of the numerical model and two systematic error types are discretization errors and modeling errors of relevance t o the simulations performed he rein (Ferziger & Peric, 2002). The accuracy of the simulations performed is assessed on the basis of evaluating each of these given types of errors. First, the discretization errors in space are minimal as highly accurate pseu do spectral methods are used for spatial derivatives as mentioned in Chapter 3. The resolution in time is around 1/100000 of the wave period which is sufficiently small to minimize the error due to temporal discretization. One source of modeling error may be the violation of underlying assumptions of the mathematical model used to describe particle laden twophase flow. The two major assumptions are: dilute particle concentration and the particles have small response time to fluid motion. The second assumpt ion is evaluated in Chapter 2 and it is concluded that particle response time is within acceptable range (Shotorban & Balachandar, 2006). On the other hand, the first assumption is dependant on the final concentration results, which is discussed in Chapter 5. Insufficient domain size and resolution might constitute another source of error. In direct numerical simulations (DNS) all the scales of turbulence should be captured. This necessitates that the grid resolution is fine enough to capture the smallest K olmogorov scales while the domain size is wide enough to incorporate the largest eddy observed for the specified Reynolds number. These criteria dictate the grid characteristics as well as the size of the domain. PAGE 45 45 One of t he purpose s of C hapter 3 is to prove the validity of the simulations against any systematic modeling errors due to insufficient domain size or the resolution to simulate the flow turbulence. In Section 3.1, the sufficiency of domain size is shown by twopoint correlation functions for each fluctuating velocity and a power spectrum analysis for the turbulent energy. In Section 3.2 mean and turbulent quantities are validated with the results of Spalart & Baldwin (1989) and the particle phase behavior is evaluated in the absence of ba ck particle effect which at least represents very dilute fine particle laden flow. These results are then used as reference to evaluate the particle effect in Chapter 5. 4.1 Grid Resolution and Domain Size Based on the descr iption of the problem given in C hapter 2, the statistically averaged flow is one dimensional and the turbulence is homogeneous in spanwise and streamwise directions. Integral length scale, defined by two point correlation function, is a common practice to evaluate the largest eddy size in the domain. Two point correlation fu nctions are defined in Eq. 4 1) and Eq. 4 2 for streamwise and spanwise directions, respectively: '''' 1 '' ,(,)(,) ()(,) (,)(,) eiiii uux iiuxrtuxt Rtfrt uxtuxt (4 1) '''' 3 '' ,(,)(,) ()(,) (,)(,) eiiii uuz iiuzrtuzt Rtgrt uztuzt (4 2) In Eq. 4 1 and Eq. 4 2 r represents the distance from a single p oint in the domain at instance t The dir ection vectors 1e and 3e in both Eq. 4 1 and Eq. 4 2 are directed towards positive x and z directions where turbulence is homogeneously distributed. The integral length scales in x xx L and z zzL defined in Eq. 4 3 and Eq. 4 4 respectively, is an indicator for the largest eddy size in the domain. PAGE 46 46 0(,)dxxLfrtr (4 3) 0(,)dzzLgrtr (4 4) Therefore, the adequacy of domain size can be tested by the two point correlation function in streamwise and spanwise direction shown in Figure 4 1 an d Figure 4 2 respectively. The results for all three components of velocity at the six differ ent phases are shown in Figure 4 1 and 4 2 The two point correlation functions are computed at two different elevations from the bed of y+ = 0.2 and y+ = 78. He re, y+ is the distance i n wall normal units defined in C hapter 3 The sudden decay of the two point correlation functions suggests that the domain lengths in the streamwise and spanwise directions are sufficient to contain all the significant turbulent flo w features. Although it contains the smallest scale turbulence information, twopoint correlation function, in real space, does not exhibit a clear description of all turbulent scales. In the spectral space, on the other hand, better information is obtaine d. The Fourier transform of the two point correlation function is equal to the power spectrum of the fluctuating velocities in all three directions (for detailed derivation, refer to Pope 2000) for a homogeneous turbulent field. This can be represented as in discretized spectral spa ce as in Eq. 4 5 and Eq. 4 6 : ''''* ,()iiii uuxuuRt (4 5) ''''* ,()iiii uuyuuRt (4 6) In Eq. 4 5 and Eq. 4 6 cap on top of ui stands for its Fourier transform The term ui represents the complex conjugate of ui and stands for the wavelength. Each wavelength can be calculated in the streamwise and spanw ise direction as in Eq. 4 7 and Eq. 4 8 respectively: PAGE 47 4 7 ,2x ix xNLi (4 7) ,2z iz zN Li (4 8) The representation of the two point correlation function in spectral space gives a better picture of al l turbulent length scales (Pope 2000). The power spectrum of fluctuating velocities reflects the energy from largest eddies to inertial subrange to Kolmogorov scales in a homogeneous turbulent field. Therefore, current simulations should reflect each of these features in streamwis e and spanwise directions where turbulence is homogeneous. The energy spectra of the three fluctuating velocity components at different phases for elevations y+ = 0.2 and y+ = 78 are shown in Figures 4 3 and 4 4 respectively. All these three energy ranges of homogeneous turbulence can clearly be observed. At the phases when the flow reaches its maximum, the inertial subrange is emphasized compared to other phases. However, the inertial range is very subtle at the flow reversal. This is due to the tendency of decay in turbulence at the intermittently turbulent regime which has also been reported by Hino et al. (1976). The decay rate at the Kolmogorov scale is a measure to evaluate the resolution of the domain and demonstrate its validity. It is observed that all the spectra show 103 to 106 magnitude decay in energy at different phases, indicating that the turbulent flow features are well resolved with the present grid resolution. 4.2 Comparison of Mean and Turbulent Quantities The numerical flow solver is validated with prior DNS re sults of Spalart & Baldwin (1989) for an oscillatory boundary layer of Re = 1000 without particles. Among all the available three dimensional simulations of oscillatory boundary layer in the lit erature, Spalart & Baldwin (1989 ) u sed the most conservative domain (the largest domain size) and the highly accurate pseudo PAGE 48 48 spectral flow solver. There are also experimental measurements in the literature for Reynolds numbers cl ose to 1000 (e.g., Jensen et al. 1989; Chen et al., 2007). In both these studies, velocity is measured by acoustic Doppler velocimet er (ADV). However, Jensen et al. (1989) discuss in their work that ADV measurements very close the bottom wall should be used with caution as reflection from the smooth wall may contamin ate the observed turbulent quantities. Figures 4 5 to 4 8 show the comparison between the present simulation results and those of Spalart & Baldwin (1989) for statistically averaged streamwise velocity, streamwise, spanwise and wall normal component of turbulent fluctuati ons at six different phases of / Excellent agreements are observed. According to results shown in Figure 4 6, streamwise turbulent intensity is strong near the end o f the acceleration phase ( = ) and the beginning of the deceleration phase ( /6). However, such turbulence decays rather quickly and near flow reversal ( /3), turbulence become small. These features are qualitatively similar with that observed by Hino et al. (1983) for sligh tly smaller Reynolds number (Re = 876) in a wind tunnel. Figure 4 9 presents a more detailed view of the threedimensional turbulent flow structures at the six difference phases during a cycle. Plotted in F igure 4 9 are vortex stru ctures identified with the iso surface of swirling strength ( ci; see Zhou et al., 1999; Chakraborty, Balachandar, & Adrian 2005). The swirling strength is defined as the imaginary portion of the complex eigenvalues of the local velocity gradient tensor. If all three eigenvalues are real then swirling strength is zero. Swirling strength has been tested to extract vortex structures in turbulent flows. The vortices are bounded within a close distance from the bottom wall at the flow peak. Then they start to swell towards the midchannel at the flow reversal and subsequently lose their strength. Consistent PAGE 49 49 with the turbulent intensity shown in Figure 4 6, the ci value is much smaller near flow reversal, suggesting low level of turbulence. Figure 4 10 presents iso surfaces of particle concentration at C/ aveC =0 .65, where aveC is the volume averaged concentration in the computational domain. The results are shown again at six different phases. It is observed that sediments are suspended at all phases but with different characterist ics. During the peak flow ( ) where turbule nce intensity is high, the particles are suspended intensely with small pocket shaped structures similar to the turbulent vorti cal structures shown in Figure 4 9 Du ring deceleration phase, (from = 0 to = start to lose their intensit y and become less energetic. However, at flow reversal ( /2) more organized spatially patchy bursts can still be clearly seen. At the early acceleration phase ( = /3) spatially patchy bursts further organized into several long streaks. Rapid turbul ence decay during the deceleration phase and weak turbulence during flow reversal shown previously discourage turbulent suspension and promote gravitational settling and result in the observed increase of near bed concentration. Proving the adequacy of the domain size and resolution, in Chapter 5 the simulations for different Ri are going to be presented. PAGE 50 50 Figure 4 1. Two point correlation functions in spanwise direction at y+ = 0.2 and y+ = 78 for phases (A) Figure 4 2. Two point correlation functions in streamwise direction at y+ = 0.2 and y+ = 78 for phases (A) PAGE 51 51 Figure 4 3. Power spectrum of fluctuating velocities in streamwise and spanwise direction for phases at y+ = 0.2 A) C) Figure 4 4 Power spectrum of fluctuating velocities in streamwise and spanwise direction for phases at y+ = 78 ( A) (B) (C) PAGE 52 52 Figure 4 5. Phase and planar average of streamw ise velocity ( U ) compared with that of Spalart & Baldwin (1989) (SB198 9 ). Figure 4 6. Phase and planar average of s treamwise turbulent intensity ( Urms) compared with that of Spalart & Baldwin (1989 ) (SB1989). PAGE 53 53 Figure 4 7 Phase and planar average of vertical turbulent intensity ( Vrms) compared with that of Spalart & Baldwin (1989) (SB198 9 ) Figure 4 8. Phase and planar average of spanwise tur bulent intensity ( Wrms) compared with that of Spalart & Baldwin (1989) (SB1989) PAGE 54 54 Figure 4 9. Vortex structures identified by swirl intensity ( ci)2 (Zhou, et al. 1999) at different phases. Figure 4 10. Is o contours of concentration at C = 1.3 at different phases. PAGE 55 55 CHAPTER 5 TWO WAY COUPLED SIMULATI ONS : EFFECT OF PARTICLE V OLUME CONCENTRATION The pre sent formulation of the problem yields four key parameters: the Reynolds number (Re), the bulk Richardson number ( Ri), the nondimensional settling velocity ( SV ) and the Schmidt number ( Sc ). The turbulent nature of the flow is dict ated by the Reynolds number with the intensity of turbulence increasing with increasing Re. The stabilizing effect on flow turbulence is due to particle induced density stratification. It can be deduced that turbulence modulation is a ffected by two parameters: Ri and settling velocity ( SV ). Settling velocity determines the vertical variation of particle concentration profile, while Ri can be thought of as a proxy for volume averaged particle concentration and to directly affect the carrier flow in the momentum equation. A s a first step, in this study the focus is on an oscillatory boundary layer of Re = 1000 and the nondimensional settlin g velocity is held constant at Vs = 904. Four different values of bulk Richardson number are considered ( Ri = 0, 104, 3104 and 6104) in this study. These four simulations are denoted as set 1. As discussed previously the clear fluid (Ri = 0.0) case is in the turbulent regime On the other hand, the flow is fully laminari zed at Ri = 6104 and thus a range of behavior is encountered in between in this study 5.1 Laminar Solution The role of the above defined parameters can be explored in the laminar solution. The momentum equation ( Eq. 2 3) can be solved under laminar flow condition and in the limit of channel half height being much larger than Stokes layer thickness the nondimensional streamwise velocity in the bottom half of the channel can be expressed as : 22 (,)sin sin Re Reytt Uyt ey (5 1 ) PAGE 56 56 The laminar solution shown above satisfies no slip boundary condition at the bottom boundary ( y = 0) and exponentially approaches the free stream oscillatory flow. From the choice of length scale, the nondimensional thickness of the laminar boundary layer is unity. The lami nar concentration profile is independent of the laminar velocity and can be expressed as : (,)expReWpCytCVScy (5 2 ) where the nondimensional concentration at the wall can be expressed as 160Re1exp60ReWp pCVSc VSc (5 3 ) For a sufficiently deep channel (i.e., 30Re1pVSc ) the wall concentration becomes 60ReWpCVSc Note that the laminar concentration profile is independent of the bulk Richardson number and for the present choice of parameters, the wall concentration is 27 times the volume averaged concentration. The laminar concentration also decays exponentially away from the bottom boundary and the concentration boundary layer thickness, 1RepVSc is 2.22 times the Stokes layer thickness. From its definition shown in Eq. 2 6 the bulk Richardson number can rewritten as : 218 Res aveV Ri C d (5 4 ) Thus for a given oscillatory flow condition and particle prop erties (i.e., fixed settling velocity, Reynolds number and / d ), Richardson number becomes a measure of volume averaged particle concentration, aveC in the domain. It is of interest here to investigate the critical values of particle concentration (or Ri) at which turbulent flow and particle transport dynamics become drastically distinctive. For example, it is the interest to establish the critical concentration where the effect of particle density stratification is negligibly small on the carrier PAGE 57 57 flow. Also, it is of practical interest to evaluate the critical concentration (or Ricrit) that causes flow laminarization and complete suppression of turb ulence and mixing. The results are presented in this thesis followi ng this objective. In Section 5 2, the results are first qualitatively observed in terms of turbulent vortical flow structures and isosurfaces of concentration. Mean flow quantities and tu rbulent intensities are then presented in section 5 .3 to examine the relationship between carrier mean flow, particle concentration and secondorder turbulent statistics. The instability mech anisms are analyzed in section 5 .4 and the flux and gradient Rich ardson nu mbers are presented in section 5.5. In Section 5.6, turbulence modulation is further analyzed through the TKE budget equation. 5.2 Turbulent Structures Simulation results for the four different particle concentrations ( Ri = 0.0, 104, 3104, 6104) present rather different outcomes of particle turbulence interaction, which shall be qualitatively identified and described in terms of isosurface s of swirling strength (Figure 5 1 ) and particle concentration (Figure 5 2 ). Swirling strength ( ci ) is the imaginary part of the complex eigenvalue pair of the local velocity gradient tensor. If all three eigenvalues of the velocity gradient tensor are real, then local swirling strength is zero. The swirling strength is able to discrim inate against regions of shear and extract only coherent vortex structures in turbulent flows (Zhou et al. 1999, Chakraborty, Balachandar & Adrian, 2005). Turbulent structures shown here do not include the case of Ri = 6104 as the turbulence is fully suppressed by the suspended particles. Hence, the flow is completely laminarized and there are no vortex structures. The effect of increasing particle induced density strati fication is observed in Figure 51 where iso surface of constant swirling strength is plotted at six different phases spanning half a wave cycle. For the case of Ri = 0.0 case 1 particles do not affect the carrier phase and thus PAGE 58 58 turbulence corresponds to that of clear fluid discussed by Spalart & Baldwin (1989 ). Let us first discuss this clear fluid limit. At = /6 the far field (in the present case the mid channel) is just past its peak forward velocity and a dense pack of near wall vortical structures can be observed close to the bottom boundary. The vortex structures observed at = when the far field reverse flow reaches its peak, are qualitatively similar. Owing to temporal periodicity vortex struc tures at = 0 can be expected to statistically similar to those observed at = but with reversed orie ntation along the streamwise direction and thus are consist ent with those observed at = /6. As th e free stream decelerates ( = /3) the vortex structures thin out and move away from the bottom boundary. This process of outward migration conti nues durin g flow reversal ( = /2) and even past flow re versal ( = 2 /3). Clear signatures of typical shear layer turbulent structures, such as quasi streamwise vortices and hairp in vortices are observed. At = 5 /6 when the reverse flow is accelerating in streng th, only few streaky vortical structures can be observed. These streaky vortices are close to the bottom boundary and it appears that between = 2 /3 and = 5 /6 the turbulence that was observe d away from the boundary at = 2 /3 decays rapidly and the near wall structures observed at = 5 /6 grow rapidly into mature turbulence when the reverse flow peaks at = The iso contours of concentration for the passive particle case ( case 1 ) shown in Figure 5 2 A depicts a similar picture. The strong turbulent near wall fluc tuations can be seen during ( = ) and immediately fo llowing the peak velocity ( = /6). The intensity of these fluctuations somewhat decrease during the deceleration phase, but persist even slightly after flow reversal. The impact of near wall streaky vortical structures on particle concentration can be clearly observed at = 5 /6. PAGE 59 59 The vortex structures for the case of Ri = 104, case 2, show n in Figure 5 2 B are statistically similar to those for the passive case. Similar observations can be made for the iso surface of particle concentration (Figure 5 2 B) with the exception that the iso surface level for case 2 is chosen to be two times the volume averaged concentration, which is slightly higher than that for the passive case. As will be dis cussed in detail when the mean concentration profiles are presented increasing Richardson number from 0 to 104 has a strong influence on mean particle concentration. Sin ce particle concentration in the bottom half of the channel increases substantially, a correspondingly higher contour level must be chosen to extract similar turbulent structures. For the case of Ri = 3104, case 3, vortex structures show a drastically different behavior th an those observed in the case 1 and case 2 First of all, there is no turbulent vortica l structure observed in Figure 5 1 C at peak flow ( = ) and during the enti re deceleration phase (see = /6 and = /3). This can be explained b y turbulence suppression due to particle induced stable density stratification. Although at flow reversal ( = /2 in Figure 5 1 C) coherent vortex structures are not observed, the effect of incipient instability in terms of small amplitude waves can be obs erved in the iso surface of pa rticle concentration in Figure 5 2 C. These instabilities tend to grow from an earli er near laminar condition ( = /6, /3) and reach peak amplitudes in the earl y acceleration phase around = 2 /3. The instability appears as well organized spanwise vortex rollers formed as a result of shear instability. Signature of secondary instability can be seen as spanwise waviness of the vortex rollers (see = 2 /3 in Figure 5 2 C). While relatively clean, small amplitude, and nearly two dimensi onal waves can observed at = /2, a short time later at = 2 /3 a more complex wave train with spanwise disturbances resembling wave breaking formation is observed. At l ater stages of acceleration ( = 5 /6), both the vortex PAGE 60 60 structure (Figure 5 1 C) and the iso su rface of concentration (Figure 5 2 C) show decay of the instabilities and only a remnant of the vortex structure and its impact on particle concentration can be observed. Interestingly, at the end of acceleration when flow reaches pack reverse flow ( = ) all disturbance seems to have decayed and the flow appears to be in laminar state. Note that the concent ration contour shown in Figure 5 2 C is ten times the volume average concentration, which is much higher than those of previous case s of lower Ri As shall be seen in section 3.3 for case 3, the flow largely laminarizes and as a result mean particle concentration close to the bottom wall significantly increases. In order to extract the effect of instability, which is located close to the bottom wall, it is necessary to choose a larger concentration value for the iso surface of concentration field. The se instability features for case 3, observed in Figure s 5 1C and 5 2 C, resemble the shear instability noted by Strang & Fernando (2001). They investigate the behavior of these interfacial shear instabilities with respect to changes in the bulk Richardson number. While for small bulk Richardson number they observe pure KelvinHelmholtz (KH) billows, with the increase in Ri they observe a com problem it is difficult to separate out the exact characteristic of these structures and the associated instabilities as they occur in a transient manner over only a short duration near flow reversal. This is unlike the case considered in Strang & Fernando (2001) where the instabilities occur in a quasi steady manner and the effect of stratification can be considered almost constant. 5.3 Mean Quantities and Turbulent Intensity In this section, interacti on between carrier flow and particles for case 1 to case 3 are further analyzed in terms of statistically averaged particle concentration and velocity statistics. Statistics are obtained by averaging over horizontal planes and also over time instants of identical phase. PAGE 61 61 These averages will be referred to as phase average and denoted by < >. Note that the st atistics during forward flow ( /2 < < /2) are related to st atistics during reversed flow ( /2 < < 3 /2) as follows: ,,,(,),,,(,)rmsrmsrms rmsrmsrmsCUVWytCUVWyt (5 5 a) ,,,(,),,(,) UUVUCytUUVUCyt (5 5 b) These symmetries are used to further e nhance the statistics. Figures 5 3 and 54 show the averaged concentration and streamwise velocity profiles at six different distinct phases during the wave cycle for the cases of Ri = 0, 104 and 3104. Note that the concentration scale for case 3 is different from that of the other two cases in order to accommodate the rapid increase in concentration towards the bottom wall. The averaged streamwise velocity profiles for case 1 and case 2 are almost identical (compare triangle symbo ls and solid curves in Figure 5 4 ). The corresponding root mean square (RMS) turbulent velocities in the streamwise, wall normal and spanwise directions at the six different p hases are presented in Figure 5 5 Again, it can be observed that the turbulent intensities for case 1 and case 2 are of similar magnitude over much of the channel. In particular, there is little difference in turbulent intensities close to the bed between y = 0 and y = 10. However, closer to the mid channel (from y the overall magnitude of turbulent intensities become small, the turbulence intensity for case 2 ranges only between 10% to 50% of that in case 1 The effect of reduced turbulence above y concentr ation profiles shown in Figure 5 3 For the passive case, the concentration profiles are more or less linear and particles are well mixed over the entire channel, except for a thin layer of large concentration gradient very close to the bottom. On the other hand, for case 1 it can be observed that the particles are entrapped in the lower half of the domain and the concentration PAGE 62 62 profile shows a clear feature of lutocline (sharp negative concentration gradient ). This shoulder shaped concentration profile, called a lutocline, is predicted by Noh & Fernando (1991) for fine particle suspended in oscillating grid flow due to the damping of turbulence kinetic energy by particle induced stable density stratification. It is also observed in the laboratory experiment with similar oscillating grid setup (Huppert, Turner, & Hallworth, 1995). In the present simulation, the formation of a lutocline can be clearly seen for fine particle transport in the oscillatory boundary layer. Closer to the mid channel, turbulence and mixing are suppressed due to the stabilizing effect of the large concentration gradient. Suppressed turbulence further prevents upward migration of particles, which would otherwise occur driven by turbulence In this case, since suppression of turbulence occurs near the mid channel where the overall magnitude of turbulent kinetic energy is already small, turbulence modulation does not make significant difference in the mean velocity profile. M ore field eviden ces of a lutocline for wave driven fluid mud transport are discussed in Section 4. The characteristic of a lutocline in oscillatory boundary layer can be further investigated by closely examining the second quarter of the domain from the bottom boundary ( y = 15~30). Figure 5 6 shows the iso surface of concentration for case 1 and case 2 at a lower value of contour level in order to capture the turbulent structures in the region where turbulence is suppressed. A more quiescent behavior can be observed for ca se 2 compared to the passive case. Only larger wave length undulations at reduced amplitude can be seen for case 2 Thus, as Richardson number is increased from 0.0 to 104, particle concentration is more noticeably affected than turbulent vortices. Examining the iso surface of swirling strength in the second quarter from the bottom wall (not shown here) also suggests that the turbulent structures below y PAGE 63 63 = 23 are statistically similar for both case 1 and case 2 Some reduction in vertical structures can be observed only closer to the mid channel. On the other hand, for case 3 a significant difference in averaged mean velocity profiles (see Figure 5 4 ) is observed as compared to those of case 1 and case 2 The velocity and particle concentration profiles in the bottom half of the channel are very close to the laminar profiles given in Eq. 5 1 and Eq. 5 2 whi ch are also plotted in F igures 5 3 and 5 4 The turbulent intensity in the wall normal and spanwise components is nearly zero at time instants of velocit y maximum (Figure 5 5 = ) and during deceleration phase ( = /6 and /3). The RMS streamwise fluctuation is also quite small close to the bottom boundary and appears to decay during the deceleration phase. However, st arting from flow reversal ( = /2 ), the RMS velocity components begin to grow during the acceleration phase in the region close to the bottom boundary (0 y ppear to reach a peak around = 2 /3 and then start to decay. During the entire period the RMS velocity fluctuations for case 3 in the bottom half of the channel are substa ntially lower than those in case 1 and case This reduction in turbulence is consistent with the vortical turbulent st ructures presented in Figure 51 The travelling wave character of the Stokes solution can be observed both in the laminar and turbulent streamwise velocity profiles. In Figure 5 4 the variation of averaged streamwise velocity along the vertical directio n at all phases shows multiple local ma xima and minima. From Eq. 5 1 it can be seen that for the laminar solution, a local minimum starts at the bottom boundary at = /4 and travels upward away from the bottom wa ll as the phase increases from /4. The l inear migration of the location of this minimum for the laminar solution is given by and can be clearly followed during t he acceleration phase from = /2 to = In fact, in Figure 5 4 the location of this minimum can be tracked further by switching to the lowest local velocity PAGE 64 64 maximum at = /6 (owing to left right reflectional symmetry of the streamwise velocity profile). This qualitative behavior of the Stokes solution can be observed in the turbulent velocity profile as well Figure 5 7 shows the vertical location of the velocity extrema for both the turbulent (case 1 and case 2) and the laminar cases as a function of the phase. In the turbulent case a local minimum again forms at the bottom boundary slightly before /2, but it travels away fr om the bottom boundary at a faster rate. In effect, at Re = 1000 the turbulent Stokes boundary layer is about 4.5 times thicker than the corresponding laminar solution. The travelling wave character of the averaged solution is important, since it implies inflection points in the velocity profile and can contribute to onset of inviscid instability. 5.4 Flow Instability The clear appearance of coherent vortical structures in Figure 5 1 C and their temporal growth over part of the wave cycle, as implied in the iso surf ace of concentration in Figure 5 2 C raise inte resting questions regarding the shear instability. Better understanding of this instability shall shed more light on the origins of fullydeveloped turbulence in oscillatory boundary layers, with and without stratif ication. In this section, the primary concern is the instability of the laminar profile. 5.4.1 Stability of Stokes Layer The instability of the oscillatory Stokes layer and transition to turbulence have not been fully understood even in the absence of particle concentration gradient induced stratification effects. Large discrepancy exists between experimental/computational o bservations and the results of stability analysis (Blennerhassett & Bassom 2008). Initial attempts (Hall 1978) have not been fully successful in identifying unstable modes that could explain experimentally observed disturbances. Recent Floquet stability an alysis of an oscillatory Stokes layer of semi  PAGE 65 65 infinite extent by Blennerhassett & Bassom (2002) has yielded growing stationary and travelling mode disturbances with a critical Reynolds number of Re,crit = 1416. They observed the critical Reynolds number t o somewhat decrease with the presence of an upper boundary, but for the present half channel height of 30 Stokes layer thickness, the effect on critical Reynolds number is minimal. In contrast to the above theoretical result of the stability analysis, the present and prior simulations of the oscillatory Stokes layer show fully developed mature turbulence even at Re = 1000. As pointed out by Blennerhassett & Bassom (2008) this disagreement between theory and experiments/computations remains unexplained. Nev ertheless, here the results from the Floque t analysis of Blennerhassett & Bassom (2002) are used to better understand the unstable vortical structures observed for the case of case 3 From their analysis the critical streamwise wavelength can be estimated to be 15 Although the growth rate of the Floquet disturbance has not been fully mapped out over the Reynolds number and wave number parametric space, the streamwise wavelength of the most unstable disturbance can be estimated to b e similar at other Re. The four spanwise oriented vortex rolls observed in Figure 5 1 C, combined with the periodic streamwise extent of the computational domain, yields a wavelength of 16.77 Of course in the present simulation, since the streamwise length of the computational domain is chosen to be 60 the possible disturbance wavelengths are limited to discrete choices. An alternate view of the instability at case 3 is presented in Figure 5 8 where contours of span averaged spanwise perturbation vorticity is plotted on the xy plane at = 2 /3. A periodic vortex roll up can be observed and vortex structures are cent ered about y 3. From F igure 5 4 it can be observed that this location reasonable coincides with the inflection point of the mean velocity profile. The eigenfunction obtained by Blennerhassett & Bassom (2002) shows the PAGE 66 66 disturbance to be located around y 5 and occur during the accelerat ion phase of the wave cycle. Shown in Figure 5 8 B are the time history of the peak perturbation spanwise vorticity (note both the positive and negative perturbations are plotted). Rapid increase in perturbation vorticity during flow reversal and subsequent slow decay can be seen. Despite the rea sonable qualitative agreement there is a need to exercise caution in comparing the present simulation results with the Floquet stability analysis. First, the d isturbances observed in Figure 5 1 C for case 3 cannot be considered as unstable in the framework of Floquet analysis. This disturbance only grew part of the wave cycle, but decayed over the rest of the wave period. These disturbances were observed in the simulation at every cycl e, however as shown in F igure 5 8 B, where the peak spanwise vorticity was plotted as a function of time, the disturbance evolution was not strictly periodic. Second, although laminarization of the flow at case 3 results in an oscillatory Stokes boundary layer and facilitates comparison wi th the Floquet stability analysis, the presence of the strong stable density stratification has not been accounted in the stability analysis o f Blennerhassett & Bassom (2002). Clearly, the density gradient will contribute to the stability of the flow. 5.4. 2 Effect of Density Stratification The study on the effect of density stratification on the stability of shear layers dates back to the pioneering work of Taylor (1931). Miles (1961) and Howard (1961) established the condition that the gradient Richardson number, defined in terms of the present nondimensional variables as 2()gCU RiRi yy (5 6 ) greater than 0.25 everywhere in the flow to be the sufficient c ondition for stability. Thus, somewhere in the flow would indicate ins tability. Based on the laminar velocity and PAGE 67 67 concentration profiles given in Eq. 5 1 and Eq. 5 2 the gradient Richardson number as a function of distance from the bottom wall can be computed using the maximum value of velocity gradient at that location during the wave cycle as (note that laminar concentration gradient is independent of time) : 222 ,laminarReexp2RegppRiRiHVScVScy (5 7 ) In the present set of simulations, Re0.45sVSc and as a result the minimum value of gradient Richardson num ber occurs at the bottom wall and is linearly depend ent on bulk Richardson number: ,laminar(0)6.075gRiy Ri Thus, for all values of Ri considered in this study gradient Richardson number at the wall is much smaller than 0.25. Since Rig,laminar as define d above exponentially increases with y, only a layer close to the bottom boundary satisfies the condition For case 2 and case 3 the layer of fluid given by 3.88 y and 3.18 y satisfy the condition for instability. The difference in the thickness of the unstable layer, as given by gradient Richardson number, is only 20%, but the resulting change in the cha rac ter of the flow (see Figure 5 1B vs Figure 51 B) is dramatic. Of course, the condition 0.25gRi arises from inviscid analysis and thus its relevance for the stability of a layer close to a solid boundary can be questioned. More importantly, the analysis of Miles (1961) and Howard (1961) are for a ste ady base flow and their applicability for the present oscillatory flow becomes relevant only if the time scale of instability growth is much faster than the period of oscillation. It can be conjectured that at Re=1000 with decreasing stratification effect (decreasing Ri ) the growth rate of instability (seen in F igure 5 1 C) can increase dramatically to trigger onset of turbulence quickly within part of the wave cycle. Thus, nonlinear effects can be expected to set in rapidly, and a Floquet analysis that ass umes a linear behavior for the disturbance over the entire wave cycle may not capture this fast transition to PAGE 68 68 turbulence. Clearly further investigation is needed on the instability mechanisms of an oscillatory boundary layer and the effect of stable strati fication. 5.5 Flux and Gradient Richardson Numbers Particle settling leads to density stratification, which attenuates boundary layer turbulence and the suppressed turbulence further encourages particle settling. According to such paradigm, it is believed that there exists a critical level of particle concentration for a given turbulent flow such that major collapse of turbulence and particle settling may occur. Such behavior in steady flow has been well described by flux and gradient Richardson numbers (e. g., Winterwerp 2001). It should be noted that in the present study, the forcing of the flow varies in time, which means that the competition between these mechanisms also varies over time. Therefore, interrelation among particle settling flux, turbulent suspension flux, and turbulence production are more complex. Nevertheless, the role of particle induced density stratification on turbulence modulation through gradient Richardson number, Rig, and flux Richardson number, RiF are examined. For the present case of nonstationary turbulent flow, Rig and RiF are defined in terms of phase averaged flow quantities as follows: 2 gC y RiRi U y (5 8 ) '' '' FCV RiRi U UV y (5 9 ) where primes i ndicates perturbation from the phase average. Both Rig and RiF are measures of turbulence attenuation due to stable density stratification versus turbulence production. Under the assumption that turbulent stresses and turbulent suspension flux can both be represented by PAGE 69 69 gradient diffusion process (i.e., / CVdCdy and / UVdUdy ) these two Richardson numbers can be shown to be qualitatively equivalent. Although gradient Richardson number Rig relies on the validity of gradient diffusion assumption, it is a commonly used parameter because it is easier to measure in laboratory and field experiments. Both Rig and RiF are presented for the case 2 and case 3 at six different phases durin g the wave cycle (see Figures 5 9 and 5 10). For the case 2 (Figure 5 9A and Figure 5 10A ), it can be observe d that the magnitudes of Rig and RiF are generally small below y = 12~18 (the thin vertical line represents the critical value of 0.25 to be used as a reference). The small values of Richar dson numbers suggest that in the lower portion ( y <12~18) of the particle laden oscillatory boundary layer, flow remains turbulent. The domination of turbulence is in accordance with the turbulent intensities discussed in Section 3.3. On the other hand, be tween y = 12 and y = 18 both Rig and RiF increase considerably, suggesting damping of turbulence due to stable density stratification. This regime also coincides more or less with the location of lutocline where large negative particle concentrati on gradie nt exists (see Figure 5 3 ). Between y = 15 to 24, RiF sometimes takes values that cannot be plotted on loglog scale, which can be seen as gaps in the plots at these elevations (see Figure 5 10). Such out of range values of RiF obtained here are due to negative (or near zero) particle suspension flux < CV >. The ratio of particle suspension flux to the set tling flux RSTS, is shown in Figure 5 11 A This ratio, by definition, is identically zero at the bottom wall. It quickly increases away from the wall and roughly reaches a constant O (1) in the lower portion of the boundary layer duri ng the deceleration phase ( = /6, /3). At = this layer is somewhat narrow and extends from the wall to only about y < 4.5. During decel eration and flow reversal ( = /6, = /3 and = /2) the thickness of this layer, where particle turbulent concentration flux nearly PAGE 70 70 counterbalances the settling flux, increases and extends to about y < 18. Immediately following flow reversal, turbulent flux decreases in the acce leration phase. Once again only at the end of the acceleration phase ( = ) the turbulent flux within the layer close to the wall substantially increases. This cyclical waxing and waning of turbulent flux is consistent with the intensity of turbul ent stru ctures shown in F igure 5 1 The behavior for y > 18 is different. During substantial portion of the wave cycle the turbulent particle flux is quite small and in fact < CV > at times becomes negative. In this region, turbulence effect is considerably weake r and particle motion is dictated by the instantaneous local downward advective motion, especially immediately after flow reversal. As a result the net plane averaged concentration flux in the vertical direction is negative. Clearly such downward advection cannot be described by gradient diffusion assumption. The above observation can be interpreted as turbulence production to be significant only over region 0 < y < 18. Note that from the averaged streamwise velocity profile (Figure 5 4 ) it can be seen that for Ri = 104, y = 18 is the approximate upper edge of the oscillatory boundary layer, above which the velocity remains nearly a constant. For the case 3 it can be seen that the regime of small Rig and RiF exists only very close to the bed ( y = 0 ~ 3). The sharp increases in the Richardson numbers for case 3 are due to local averaged velocity gradients becoming zero rapidly away from the bed. The velocity gradients become zero at where the averaged streamwise velocity reaches maximum and minim um values. As discussed above in Section 3.3, the initiation of a velocity minimum at the wall around = /4 and its migration away from the bottom wall with time can be clearly followed in Figure 5 9 B in terms of the location of the sharp peaks. As seen in the averaged streamwise velocity profiles, such maxima and minima are also observed for case 2 however, owing to stronger PAGE 71 71 turbulent mixing near the bed, the regime of sharp increase of Rig is located farther away from the bed. Above y 3 both Rig and RiF maintain their large magnitudes ( Rig and RiF >> 0.25) for the case 3 Again, this is consistent with the location of large negative concentration gradient observed near the bed (see Figure 5 3). In Figure 5 1 C, shear instability generation is observed during flow reversal ( = / 2), which intensi fies at early acceleration ( = 3 /2) with a subseque nt decay till peak velocity ( = 5 /6, ). This interesting time dependant feature cannot be observed in the gradient Richardson number profiles as they are qualitatively dictated by the laminar profiles at all phases. This is because the fluctuation level is not strong enough to create major change in the statistically averaged quantities. On the other hand, flux Richardson number profiles show a slightly bett er picture of the instantaneous mixing process. Velocity fluctuations appear to have a greater effect as shear instability intermittently enhances perturbations during flow reversal. However, the overall level of velocity fluctuation is relatively low in t his case and hence both Rig and RiF are not very effective parameter to identify shear instability. On the other hand, RSTS, is a good parameter to identify shear insta bility during flow rev ersal as observed from Figure 5 11B. Suspension flux is very small across the channel over much of the wave cycle when compared to the settling flux. At flow reversal and during the acceleration phase, the suspension flux become lar ge near the bed and in fact at = 2 /3 the ratio of particle suspension flux to s ettling flux reaches 0.75 at y = 2.5. The vertical location of this increased suspension flux corresponds well with the location of the vor tex structures seen in Figure 5 8 5.6 TKE Budget Simulation results presented so far identify several critical phenomena, s uch as the formation of lutocline and occurance of instability during flow reversal, as Ri (i.e., particle PAGE 72 72 concentration) increases. It is clear that these phenomena are due to different degrees of particle turbulence interaction in the oscillatory boundar y layer. In this section, the mechanisms causing these critical phenomena can be explained via the budget of turbulent kinetic energy (TKE), denoted as k in this study. The TKE budget equation for the present problem can be written as: '' '' '' 2 '2 ''' '' 1 21 Re 111 2Re2i ii ij ii j jj ii iij ii ij ju uu k uu RiCu tx xx P u uuu uu xxx (5 10) The term on the left hand side is the time derivative of TKE. The six terms on the righthand side are turbulent production, turbulent dissipation, particle induced buoyancy dissipation/production, pressure transport, turbulent advection, and viscous diffusion, respectively. For convenience, pressure transport, turbulent advection, and viscous diffusion (the last three ter ms on the righ t hand side of Eq. 5 10) are collectively presented as total transport in the TKE budget p resentat ion (see Figures 5 12~5 21). For case 1 where particles are passive to carrier flow, the leading contributions are from turbulent production and turbulent dissipation within the oscillatory boundary layer (see Figure 5 12 for y = 0~12). Production and dissipation are more or less in balance with each other except within y < 1 (corresponds to within 50 wall unit), where total transport becomes comparable especially during and after flow reversal ( = /2, 2 /3, 5 /6). The time ra te of change in TKE (i.e., / kt ) is relatively unimportant within the oscillatory boundary layer. The temporal flow characteristics of the oscillatory flow can be observed from turbulence production which correlates with mean velocity gra dient. From flow reversal ( = /2) to dur ing the accelerating phase ( = 2 /3, 5 /6, ), the production term increases significantly (by three order of magnitudes) close to the bottom wall and fades away moving away from the wall. On the other PAGE 73 73 hand, during deceleration ( = /6, /3, /2), turbulent production near the bed almost fade way completely (Hino et al. 1983). Near the top and above the oscillato ry boundary layer (see Figure 5 13A for y = 15~25), a very different TKE budget balance can be observed. Firstly, the intensities of these terms are two to three orders of magnitudes smaller than those within the oscillatory boundary layer. Moreover, turbulent production is of less importance due to lower shear above the oscillatory boundary layer and hence it is the observation that importance of total transport increases which entrains turbulent kinetic energy from the lower oscillatory boundary layer to balance turbulent dissipation. Time rate of change in TKE is also of more importance in this regime. For case 2 the TKE budget is similar to that of case 1 within the oscillato ry boundary layer (see Figure 5 14 for y = 0~12). This is consistent with results shown in Figure 5 4 that the averaged velocity profiles for case 1 and case 2 are almost identical. Particle induced buoyancy dissipation term is negligible in the TKE budget in this regime. On the other hand, the TKE budget becomes very different near the top and above the oscillato ry bounda ry layer (see Figure 5 13B for y = 15~25). Particle induced buoyancy dissipation contributes to about 20% to 50% of the turbulent dissipation. Comparing to Ri = 0.0 (Figure 5 13A), the overall intensities of these terms are smaller. Most noticeably, the tu rbulent production for case 2 is about 40% to several factors smaller than that of case 1 Also between y = 15~25 the production term penetrates towards the midchannel more than the case 2 This reduction in turbulence production explains the observed re duction in turbulent fluctuation (see Figure 5 5 ) and more importantly the formation of lutocline near the top of the oscillatory boundary layer. It is emphasized here that for case 2 because the region where turbulent production is significantly reduced is already of very small magnitude compared to that near the bed, the total turbulence production in the entire PAGE 74 74 domain is almost unchanged. Hence the averaged velocity profiles for case 1 and case 2 are similar. Another observation on the production term i s the change in the sign along the vertical direction for case 1 For example, at the phases = /6 and = 5 /6 between y = 18 to y = 22 negative production is observed. This is different from the production results of case 1 where the production stays p ositive along the vertical direction at all the phases. This means the phase and the planar averaged derivative of the mean flow along vertical direction and the Reynolds stresses become of the same sign be tween y = 18~22 at phases = /6 and = 5 /6. In case 3, a different picture in TKE budget is observed (see Figure 5 15) throughout the entire domain as compared to the other two cases. Firstly, the magnitudes of all terms in the TKE budget are significantly smaller at all the phases due to suppression of turbulence by strong particle induced stable density stratification. However, at the onset of instabi lity during flow reversal ( = /2), turbulent production is observed to occur around y = 1~3, which is consistent with the location of inst abi lity shown in Figure 5 14. Due to such abrupt generation of turbulence, time rate of change in TKE as well as total transport become the more important term (than turbulent dissipation) to balance production. Such feature also triggers the upward partic le burst observed during the flow instability. It should also be noted particle induced buoyancy loss in TKE becomes comparable to viscous dissipation with the onset of instability which falls between = /2 and = 2 /3. Also peak dissipation due to part icle buoyancy moves with the suspended particles in the vertical direction from the point where KH billows are first observed. Finally, because turbulent production is significantly reduced throughout the entire domain, flow tends to laminarize (see averag ed ve locity shown in Figure 5 4 ) and mixing is suppressed. PAGE 75 75 5.7 Flow Regimes and Field Interpretations Numerical simulations presented here qualitatively identify several regimes of fine particle transport and their corresponding dynamics in oscillatory cha nnel flow. These regimes are consistent with prior laboratory and field studies on fine sediment transport in wave boundary layer. Hence, this section provides a summary of flow regimes revealed by our numerical results and their field implications. It is illuminating here to convert, at least qualitatively, the bulk Richardson number Ri to the near bed sediment mass concentration in the wave boundary laye r. As mentioned, the present Re = 1000 gives an amplitude of oscillatory velocity of 0.56 m/s and Stok es boundary layer thickness of 0.0018 m. Using Ri given in Eq. 2 6 and a sediment density p = 2650 kg/m3, the volume averaged concentrations for case 2 and case 3 can be deduced to be aveC = 0.0011 and 0.0033, respectively. The volume averaged concentration is used to normalize the particle concentration in the calculations. Hence, according to the averaged concentratio n profiles presented in Figure 5 3 the near bed sediment mass concentration for Ri = 104 and 3104 can be estimated to be aveC = 10 and 50 g/l, respectively. According to the simulation results, four flow regimes are observed as bulk Richardson number increases from clear fluid condition ( Ri = 0.0 ca se 1 ). These regimes are described in sections 5.7.1 to 5.7.4. 5.7.1 Regime 1 Regime I describes a condition of very dilute sediment concentration ( Ri ~ 0) where the effect of sediment induced density stratification is negligible and sediment is passive t o carrier flow. Sediment concentration profile is more or less well mixed in and above the wave boundary layer (See the concentration profile of case 1 in figure 5 3 ). Simulation results suggest sediment concentration for this regime to be much smaller tha n 10 g/l. In fact, field observations identify PAGE 76 76 such regime as low concentration mud suspension (e.g., Winterwerp & van Kesteren 2004) where near bed concentration is no more than O (1) g/l. 5.7.2 Regime 2 Regime II describes moderate interaction between sediment and carrier flow as sediment concentration near the bed exceeds several grams per liter. Sediment induced stable density stratification attenuates carrier flow turbulence near the top of the wave boundary layer. Hence, the associated turbulent mix ing of sediment is significantly suppressed causing pronounced shoulder shape concentration profile, characterized by a rapid decrease in concentration near the top of the wave boundary layer. On the other hand, because the magnitudes of turbulent produc tion and turbulent kinetic energy around the top of the wave boundary layer are already small comparing to those within the wave boundary layer, the averaged flow velocity and turbulent intensity are only slightly modulated. The shoulder shape concentrat ion profile (See the concentration profile of case 2 in figure 5 3) called a lutocline, exist s not only for fine sediment suspension in an oscillating grid (Huppert et al. 1995), steady current or tidal flow (Trowbridge & Kineke 1994; Winterwerp 2001) but also observed is in the field and laboratory during wave driven fine sediment transport with a near bed concentration more than 10 g/l (e.g., Ross & Mehta 1989; Tray kovski et al. 2000, Traykovski et al. 2007; Lamb, DAsar o, & Parsons 2004). Our simulation results are consistent with these observations on wave induced fluid mud transport. 5.7.3 Regime 3 In regime III, near bed sediment concentration is around several tens to 100 g/l where strong coupling between carrier fl ow and sediment transport is observed. Due to the increasing effect of a stable sediment concentration gradient that directly suppresses turbulent production in the wave boundary layer, the mean flow velocity tends to laminarize and causes a smaller wave PAGE 77 77 b oundary layer thickness. This feature of carrier flow is quite different from that of regime II. Moreover, as turbulent suspension of sediment in the wave boundary layer is shut down during a portion of the wave period, the flow destabilizes during flow re versal and episodic sediment burst s and high turbulence occur via flow instability. Subsequently during the rest of the wave cycle the instability decays and a quiescent laminar flow is observed. From a comparison of averaged velocity and concentration profiles with the laminar counterparts, it is observed that they are nearly identical. It is not completely clear if such occurrence of instability during flow reversal has been observed in the field or laboratory under the flow and sediment conditions considered in this study. Foster, Beach, & Holman (2006) observe similar episodic near bed sediment burst events occur only during flow reversal in a sandy beach surf zone of two meter flow depth. However, both the Stokes Reynolds number (Re~2000) and nondimens ional settling velocity ( Vs = 203) are about 2~3 times larger than the conditions considered in this study. One might argue that the larger Reynolds number contributes to enhanced turbulence and as a result larger particle settling velocity (or larger c oncentration) is required to establish the condition of instability during flow reversal. The effect of increased particle settling velocity at a fixed bulk Richardson number on turbulent suppression and laminarization has been discussed by Cantero et al. (2009) in the context of turbidity currents. Traykovski (personal communication ) observe lutocline wave s during wave driven fluid mud events. They report lutocline waves that have wave lengths much smaller than the surface wave length. However, the near bed sediment concentration that they observe is of few hundred g/l (Traykovski, personal communication) and hence the measured lutocline wave may be mostly PAGE 78 78 controlled by rheological stresses and associated non Newtonian flow effects, which are not incorpora ted in the present study. 5.7.4 Regime 4 In regime IV, where near bed sediment concentration is greater than O(100) g/l, sediment induced stable density stratification is severe enough that it further suppresses the instability observed in regime III and the flow remains laminar at all times. Complete laminarization of turbulent wave boundary layer due to sediment induced stable density stratification has been hypothesized in many models to study wave dissipation over muddy seabed. In such studies the wav e boundary layer is calculated based on a laminar formulation with enhanced effective viscosity or non Newtonian (e.g., Bingham plastic) viscosity that further damps the wave energy (e.g., Dalrymple & Liu 1978; Mei & Liu 1987). Our simulation results discu ssed here support such laminar flow assumption. In a laboratory study of fine sediment transport in an oscillatory tunnel, Lamb, DAsaro, & Parsons (2004) report as wave velocity amplitude increases, wave boundary layer thickness also increases and more se diment is suspended from the bed. However, a critical wave velocity amplitude seems to exist such that when the wave amplitude is further increased, wave boundary layer turbulence eventually collapse and wave boundary layer thickness decreases. This observ ation of carrier flow behavior is consistent with our simulation results shown here. In many coastal modeling applications, the effects of wave boundary layer are parameterized by a friction factor (or by a drag coefficient). It is well known that when the presence of sediment attenuates carrier flow turbulence, the friction factor experienced by overlaying hydrodynamics is reduced (e.g., Glenn & Grant 1987; Style & Glenn 2002). Our simulation results suggest such paradigm for fine sediment may need to be r efined for wave induced fine sediment transport. Sediment induced stable density stratification affects sediment PAGE 79 79 suspension and carrier flow in wave boundary layer in two different levels. For milder concentration, there exists a regime (regime II) of weak turbulence attenuation where sediment mixing is suppressed only at the top of the wave boundary layer to form a lutocline. However, the mean flow field is intact and hence there is no need to revise the friction factor associated with mean flow energy dis sipation. Only at higher concentration, such as regime III and IV, turbulence is significantly attenuated and friction factor must be reduced. In most field studies on wave propagation over m uddy seabed (e.g., Forristal & Reece 1985; Elgar & Raubenheimer 2 0 08), there are mostly reports of increased energy dissipation due to the presence of mud in the wave boundary layer instead of drag reduction. Our simulation results here imply the presence of lower concentration of fluid mud may not give reduction of fri ction factor (e.g., regime II). As sediment concentration further increases such that sediment induced stable density stratification directly reduces turbulent production, the mean flow laminarizes. In the mean time, rheological stress may further become e ffective due to high sediment concentration to attenuate wave energy. Hence, it is likely that the range of sediment concentration in the wave boundary layer that attenuates turbulence and causes reduction of friction is limited. PAGE 80 80 Figure 5 1. Turbulent vortical structures extracted in terms of iso surfaces of swirling strength (ci). for cases (A) Ri = 0, (B) Ri = 104, and (C) Ri = 3 x104. The contour level chosen is ci = (0.1~0.2)*max( c i). Vortical structures are shown at six different phases, from PAGE 81 81 Figure 5 2. Turbulent vortical structures extracted in terms of isosurfaces of particle concentration for case 1, case 2, and case3. Concentration iso contours are shown at six different phases, from Ri = 0 showing C = 1.65, (B) Ri = 104 showing C=2.0, and (C) Ri = 304 showing C = 10.0. PAGE 82 82 Figure 5 3. Mean concentration profiles for Ri = 0.0, 104 and 3104 (Refer to the upper axis for Ri = 304). Six different phases during half a wave cycle are shown and the results for the other half of the wave cycle can be obtained from symmetry. PAGE 83 83 Figure 5 4. Mean streamwise velocity profiles for Ri = 0, 104 and 3104. Six different phases during half a wave cycle are shown and the r esults for the other half of the wave cycle can be obtained from symmetry. PAGE 84 84 Figure 5 5. Turbulent rms velocity fluctuation along (A) streamwise, (B) wall normal and (C) spanwise directions. Six different phases during half a wave cycle are shown for Ri = 0, 104 and 3104. PAGE 85 85 Figure 5 6. Iso surfaces of concentration at six different phases from 1 and case 2. (A) Ri = 0 showing C=1.3, (B) Ri = 104 showing C = 1.3. Note only the middle third of the bot tom half of the channel from y = 15 to y = 30 is shown. PAGE 86 86 Figure 5 7. Time variation of extremum velocity for Ri = 104 and Ri = 3 x 104. Note that there exists phase slots where 2 extremum points observed for both cases which can also be seen planar and phase averaged velocity profiles. Figure 5 8. Spanwise vorticity contours and time series of maximum and minimum perturbation vortcity for case 3. (A) Contours of span averaged spanwise vorticity at (B) time history of the maximum vorticity o bserved in the aforementioned region. PAGE 87 87 Figure 5 9. Plot of gradient Richardson number at six different phases for (A) Ri = 104 and (B) Ri = 304 in log scale. PAGE 88 88 Figure 5 10. Plot of flux Richardson number at six different phases for (A) Ri = 104 and (B) Ri = 3104 in log scale. PAGE 89 89 Figure 5 11. The ratio of turbulent particle suspension flux over settling flux for (A) Ri = 104 and (B) Ri = 304. PAGE 90 90 Figure 5 12. Profiles of TKE budget terms for Ri = 0.0: Plotted are the turbulent production (Prod), viscous dissipation (Vis. Diss.), total kinetic energy transport (Total Transp.), which is the sum of diffusion, turbulent advection, and pressure transport terms, and time rate of change in TKE ( k / t ). PAGE 91 91 Figure 5 13. Profiles of TKE budget terms for Ri = 0.0 and Ri = 104 at y = 15~25 : Plotted are the turbulent production (Prod), viscous dissipation (Vis. Diss.), total kinetic energy transport (Total Transp.), which is the sum of diffusion, turbulent advection, and pressure transport terms, and time rate of change in TKE ( k / t ). PAGE 92 92 Figure 5 14. Profiles of TKE budget terms for Ri = 104: Plotted are the turbulent production (Prod), viscous dissipation (Vis. Diss.), total kinetic energy transport (Total Transp.), which is the sum of diffusion, turbulent advection, and pressure transport terms, dissipation due to particles (Partcle Diss.) and time rate of change in TKE ( k / t ). PAGE 93 93 Figure 5 15. Profiles of TKE budget terms for Ri = 3 104: Plotted are th e turbulent production (Prod), viscous dissipation (Vis. Diss.), total kinetic energy transport (Total Transp.), which is the sum of diffusion, turbulent advection, and pressure transport terms, dissipation due to particles (Particle Diss.) and time rate o f change in TKE ( ( k / t ). PAGE 94 94 CHAPTER 6 TWO WAY COUPLED SIMULATIONS: EFFECT OF PARTICLE SETTLING VELOCITY P article settling velocity is one of the key parameters that affect particle distribution and turbulence particle i nteraction. Therefore, in Chapter 6 the focus is on the particle settling velocity. It should again be emphasized that settling velocity, in the formulation adopted (see Eq. 2 8), ca n be taken as a proxy for particle size. As noted in Chapter 5, the equil ibrium state in the vertical concentration profile is achieved when the settling and turbulent diffusion fluxes are in balance (see Eq. 2 8). Settling flux is settling velocity times the local volumetric concentration of particle phase. Therefore, the vari ation in particle settling velocity has a direct impact on particle transport characteristics by modulat ing the concentration profile. With increased settling velocity, particle phase concentration migrates towards the bottom boundary. This shift in concen tration profile is likely to result in more intense turbulence modulation at the bottom boundary layer and effect the location of the lutocline as well. In contrast the reduction in settling velocity would result in migration of the concentration profile i n the opposite direction of gravity and result in a milder concentration gr adient. This may discourage laminarizatio n Again, a reduced settli ng velocity would likely impact the location of the lut ocline and turbulent suspension in regime 2 as described in Chapter 5. All these described scenarios, based on the variation in settling vel ocity, provide the motivation for numerical simulations with variable particle settl ing velocity to be presented in this chapter In addition, the simula tion results for diffe rent Ri presented in Chapter 5, also suggest that similar particle transport and flow regimes are likely to exist for varying particle settling velocity Together with the simulations presented in Chapter 5, the simulations presented in C hapter 6 bound a rectangular region in parametric space (see Figure 6 1) in which a realistic combination of flow and particle parameters in a moderately en ergetic flow field is studied. In Figure 6 1, S et 1 represents the simulations PAGE 95 95 conducted with variable Ri and set 2 represen ts the simulations within a range of Vs whose results are presented in this c hapter Exactly following the example given in Chapter 5 with the same flow characteristics, Vs = 3 104 Vs = 27 104 l these particle sizes fall within the category of silt without flocculation. Therefore, in Chapter 6 the effect of fine particles of variable size on the carrier flow turbulence and its tr ans port characteristics is investigated via different particle settling velocity. The flow character istics and particle transport are qualitatively evaluated and are presented in Section 6 1 via isocontours of concentration and vortex structures. In secti on 6 2, ensembleaveraged concentration and velocity profiles are shown. Also in Section 6 2 turbulent characteristics of the flow field is evaluated with r.m.s. of all three components of velo city fluctuations Similar to the analyses in Section 5 3, in S ection 6 3 suspension and mixing characteristics are evaluated via gradient and flux Richardson numbers, i.e. Rig and RiF respectively, in addition to the ratio of turbulent to settling flux. The turbulence level and possible modulation sources of turbulence are presented in Section 6 4 via TKE budget. 6.1 Qualitative Observations: Turbulent Structures Based on the simulation resul ts presented in Chapter 5 for different Ri, similar flow and particle transport regimes are qualitatively observed for varying particle settling velocity. The simulation cases which fall into the same flow regime categories can be paired as: i) case 2 and case 5, ii) case 3 and case 6, iii) case 4 and c ase 7. While temporal and spatial differences are observed in turbulent flow fi eld cases in the same pairs are of similar characteristics The quantita tive results presented in the next three sections conform to qualitative observations presented in this section. Therefore, from this point forward, discussions on the results of a specific case are linked to the results of its pair in order to identify the distinctions and PAGE 96 96 similarities between them Turbulence activity is evaluated qualitatively through vortex structures, namely concentration iso contours and iso contours of s wirl strength, i.e. ci. The defini tion of swirl strength is given Section 5.1 and therefore is not repeated herein. As laminar flow field is achieved in cas e 7, the results are not shown in this section As shown in Figure 6 2A, concentration iso contours for case 5 ( Vs = 3 104 and Ri = 1 104) at C = 0.4 show similar patterns through the wave cycle as in that of case 2 ( Vs = 9 104 and Ri = 1 104). Note that to extract the flow characteristics the selected iso contour value, C = 0.4, in Figure 6 2A, is significantly smaller than the one in c ase 2, where C = 2.0 wa s selected. This implies particle concentration at the bed is lower in case 5 and higher particle suspension is present in the water column. Closely packed, energetic small perturbations grow into larger clusters of spikes at flow reversal. Then these spikes evolve into long streaks at ear ly acceleration and then deform again into small perturbations at t he end of the acceleration This is very similar to structures observed at early dece leration. Sim ilarly, in Figure 6 3A, dense vortex structures are uplifted in vertical direction at their intensity at H orizontally oriented vortices are observ ed at the bottom wall at phase which are likely to be responsible for the long streaks observed in the concent ration iso contours. Again, at very dense vortices emerge due to increase in far field velocity and hence higher turbulence. Concentration iso contours for case 6 ( Vs = 2 7 104 and Ri = 1 104) are shown in Figure 6 2B. As it is shown in Figure 6 4 where time series of maximum and minimum perturbation vorticity are presented, shear instability at flow reversal, which is the signature of the third regime described in Chapter 5, the observed vorticity magnitudes at flow reversal during the second half are higher than the ones observed in the first half for case 6. Therefore, to emphasize maximum variation in turbulence and particle suspension as a result of shear PAGE 97 97 instability, the second half of the wave cycle is selected for comparison between case 3 and case 6. As a result, the phase in c ase 6 corresponds to phase ase 3. Similar to case 3, quiescent flow field is disturb ed due to shear instability in the form of w aves of small amplitude in concentration isocontours. Then these waves grow into more complex 2D vortex structures at early acceleration followed by temporal decay of these structures. It should be noted that the orientation of the structures are in the o pposite direction due to selection of the second half of the wave cycle where the flow is in the o pposite direction to the first half of the cycle. Si milar observations are made in the iso contours of swirl strength, i.e. ci. Quasi 2D r ollers are first o bserved at flow reversal, i.e. at into more complex 3D vortices at = which can be observed in terms of large hairpin vortices generated at are observed. Though somewhat similar, there exist significant differences between case 3 and case 6. The particle concentration iso contour level selected for case 6 is C = 4.0 whereas for case 3 it is selected as C = 10.0. This is due to sharp inc rease in concentration profile very close to the bottom wall for larger particle settling velocity which is discussed in more detail in section 6.2. This results in a lower particle concentration at the location where shea r instability takes place, i.e. the upper limb of the theoretical laminar velocity profile. As seen in Figure 6 2B, small er wave amplitudes observed in the iso concentration at ompared to ones at case 3. Then they rapidly grow into m ore complex 2D structures at turbulent fluctuations at contours of C = 4.0. When compared to case 3 at case 6 are larger in size lea d ing to higher level s of turbulence. Although in Figure 6 3B turbulent structures are not observed, still there exist low level fluctuations. This suggests higher level of turbulence as PAGE 98 98 opposed to the quiescent flow observed at oncentration level for case 6 at the instability location, mentioned earlier in this section, combined with smaller Ri cause smaller particle buoyancy effect i.e. 2RiC e (See Eq. 2 3 ). Therefore, turbulence is less suppressed causing higher turbulence generation and smaller temporal decay in turbulence. This issue is further evaluated quantitatively in section 6.2. 6.2 Mean and Turbulent Quantities The effect of particle induced density stratification on turbulence modulation and the e ffect of flow turbulence on particle transport can be monitored through ensemble averaged velocity and concentration profiles as well as r.m.s. of velocity fluctuations. The details of how to obtain ensemble averaged quantities is previously described in s ection 5 2 and therefore is not repeated in this section. The computed ensemble averaged particle concentration, velocity, and r.m.s. of velocity fluctuations are given in Figure s 6 5, 6 6, and 6 7 respectively. The mean velocity profiles for case 5, sho wn in Figure 6 4 are very close to that of case 1 and case 2. The oscillatory boundary layer, again, extends from y = 0 to y = 20 and incorporates large turbulence mixing which can also be seen from the r.m.s. velocities (figure 6 6) The main difference between case 1 and case 2, as stated in Section 5 3, is the ensemble averaged concentration profiles. Similar observations are found between case 1 and case 5. Particle phase concentration is mostly entrapped within the lower half of the domain showing a s harp concentration gradient, i.e. lutocline, between the regions of high concentration and low concentrati on. Though case 2 and case 5 have many similarities, case 5 has the following distinctions from case 2: i) concentration profile below lutocline shows a more well mixed characteristic whereas in case 2 a strong er conc entration gradient is observed below the lutocline ii) the location of lutocline for case 5 is further away from the bottom extending from y PAGE 99 99 e extends from y = 15 to y = 25, iii) the lutocline has a milder concentration gradient than the one in case 2. All these observations can be attributed to the reduced settling velocity in case 5. The rate of decrease in positive y direction shows the effe ct of partic le settling velocity and is already present in case 1 (see Figure 5 3) where particles are passively settling without influencing the carrier flow turbulence. As the settling velocity becomes larger, the concentration gradient gets more pronoun ced. For lower particle settling velocity, the variation in particle concentration is small and particle phase below the lutocline, y = 5 to y = 20, is more well mixed. The mass conservation equation (see Eq. 2 8) adopted in the present formulation clearly states that the equilibrium particle concentration is achieved when settling flux is counterbalanced with particle phase diffusion, turbulent suspension and time rate of change in concentration profile. Therefore, reduced particle settling flux, due to re duced settling velocity, causes a shift in lutocline location in positive y direction to achieve the equilibrium state. Again based on the mass conservation equation the reduced settling flux is balanced with lower diffusion which causes milder concentra tion gradient at the lutocline. The r.m.s. of velocity fluctuations, on the other hand, do not show significant differences between case 2 and case 5. Shown in Figure 6 6, the r.m.s. velocities are nearly identical throughout the entire lower half of the domain. However, it should be kept in mind that turbulence level is very small at the far stream for case 1 and even smaller for case 2 and 5 where decay in turbulence is present due to particle induced density stratification. For case 6, ensemble averaged concentration and velocity profiles are very close to theoretical laminar solution (see Eq. 5 1 and Eq. 5 2 ). Although for case 6 the presence of turbulence is observed through r.m.s. velocities, its intensity does not suffice to make a significant change in the ensemble averaged velocity and concentration profiles. Similar PAGE 100 100 observations are also made for case 3. Similar to case 3, in case 6 the r.m.s. velocities increase after the onset of shear instability and consequen t turbulence generation which are qu alitatively observed in section 6.1. As also observed in section 6.1, the onset of instability takes place at = and vertical direction at the same phase. The difference between maximum Urms values at and ent from the ones observed at and Vrm s and Wrms are nearly zero before the occurrence of shear instability Further increase in the rms velocity at ase 3) and ase 6) is clearly in coherence with the qualitative observations made in section 6.1. While maximum Vrms and Wrms in case 6 are slightly smaller than the ones observed in case 3 in accelerating phases, Urms values especially at is nearly two times of that observed for c ase 3 at in section 6.1 through iso contours of C = 4.0 (compare Figure s 6 2B and 5 2C). The existence of large bumpy structures for case 6 as opposed to smal l vertical variations in iso contour of C = 10 in case 3 suggests a higher turbulence level. This is verified by higher r.m.s. velocity levels observed in c ase 6 compared to that of Case 3. The observations made heretofore are critical in understanding mixing and deposition of particle phase and is further analyzed in section 6.3. 6.3 Mixing and Suspension Mixing of the particle phase and i ts transport due to the variation of particle settling vel ocity in water column are investigated through Rig, RiF, and ratio of particle suspension flux to particle settling flux ( RSTS). The definitions of these parameter s are given in Section 5 5. Figure 6 8 presents the Rig and RiF profiles within the lower half of the domain for case 5 ( Figure 6 8A) and case 6 ( Figure 6 8B). The first distinct observation for case 5 from case 2 is the close PAGE 101 101 agreement of Rig and RiF. Within the oscillatory boundary layer, where turbulence dominates, both Rig and RiF are smaller than 0.25 (see the thin line in Figure 6 8) and nearly an exact match is observed between these two numbers. Close to the velocity peak at y = 20 these two numbers show singular like behav ior and attain high values that cannot be shown within the range of the plot. From the location of the peak velocity to the lo cation of the lutocline, i.e. between y = 20 and y = 25, close agreement persists. Just above the lutocline the values of Rig and RiF start to fluctuate and the agreement between these two numbers diminishes. The same reason as in case 2 is valid here for case 5, that is, large time dependant fluctuations dictate the transport of suspended particles within the lutocl ine in the absence of high turbulent fluctuations (See Figure s 5 9 and 5 10).In addition to Rig and RiF profiles Figure 6 9A shows the RSTS profiles at six different phases spanning in the second half of the wave cycle for case 5. Closer examination of Figure 6 9A shows that from the bottom wall, where there is no turbulent suspension flux, RSTS increases from 0 to 1 up to y unity until y this ratio remains unity between y at all the phases shown. In other words, settling flux is counterbalanced by turbulent suspension. Between y y = 30 large time depende nt fluctuati ons are present as described above. In case 5 the described fluctuations start just above the point y = 17, whereas in case 2 the location where the fluctuations start is approximately at y = 13. It is also noted that above the locations mentioned RSTS v alues vary within the range of [ 2, 4 ] for case 5, and [ 2, 2 ] for case 2. As shown in Figure 6 5A the location of the lutocline in case 5 is above the one observed in case 2. Also as previously mentioned in section 5.2 and section 5.5, the fluctuatio ns are due to time dependant local advection within the lutocline. Therefore, the location s of the start ing point of these fluctuations are closely related to the location of the lut ocline. The observed start ing point of these fluctuations in case 5 above the one PAGE 102 102 observed in case 2 is in agree ment with their corresponding lutocline location The difference in range of RSTS, however, is a ffected by the settling velocity. Since the settling velocity is smaller in case 5, the denominator of RSTS becomes smalle r and a higher variation is expected at the lutocline compared to that of case 2. Laminar velocity profil es at different wave phase show that laminar Stokes boundary layer thickness is nearly one fourth of its turbulent counterpart. Therefore, for case 6 t he velocity gradient term that appears in the denominator of both Rig and RiF leads to singular like behavior outside the oscillatory boundary layer. This is also observed in case 3 shown in Section 5 5. In case 6, until the point of peak velocity, both Ri chardson numbers are smaller than 0.25 which is taken as the reference for instabi lity initiation. Within a narrow region bounded by velocity peak (y = 2~4 in Figure 68B) and the region where far stream velocity is attained (y = 4~6 in Figure 6 8B), both Rig and RiF are close to 0.25 with very close agreement at phases and on of the shear instability at the values of both Rig and RiF become larger and become even greater at he deviation betw een Rig and RiF i s observed at until coherent with the observations for ensemble averaged me an and turbulent quantities discussed in Section 6 2 which effects the mixing and suspen sion processes represented through Rig and RiF. Therefore, it should be mentioned that both Rig and RiF are not good parameters to obtain the overall mixing characteristics as they depend on ensemble averaged velocity that lead to singular like behavior through most of the domain. RSTS profiles at diffe rent phases are further presented for case 6. Although quantitative difference is immanent, similar qualitative results are observed for phase suspension flux at y y y similar to that PAGE 103 103 of case 3. The peak value of RSTS is 0.8 for case 3 and 2.0 for case 6 su ggesting a higher suspension during the accelerating phase However, it should be reminded that the peak value of RSTS is affected by the local concentration and the settling veloc ity of the particles which may explain this significant difference. The same reason can be used to explain the RSTS profile observed at ved at y the aforementioned phase which is different from that observed in case 3. Here the ensembleaveraged concentration profile has sharp variation within the region bounded by y = 0 and y = 4. Above y = 4 the concentration becomes nearly zero and therefore causes the settling flux to suspension flux t o be comparable. As a result, a second peak is present at y = 10 in RSTS profile at which are also different from case 3. As noted in sections 6 1 and 6 2, the turbulence level after the flow reversal is higher compared to that of case 3 and the remnant of turbulence due to shear in stability is still observed at = ons in RSTS profiles. Therefore, complementary TKE budget analysis is necessary to evaluate the turbulence generating and dissipating mechanisms. 6.4 TKE Budget As previously me ntioned in section 6.1, turbulence modulation mechanisms can be explained throu gh TKE budget. The TKE budget equation is given in Eq. 5 10 and the physical meaning of each term is explained in section 5.6. Similarly, all the transport terms, including turbulent convection, viscous di ffusion, and pressure transport, are presented as a total transport. In section 5.6, the difference between case 1 and case 2 is observed within a vertical band where lutocline is observed. The main difference between case 1 and case 2 comes from the reduced production and viscous dissipation terms in case 2 compared to those in case 1 while the buoyancy dissipation also becomes comparable to production and viscous dissipation at the PAGE 104 104 lutocline. However, as observed from the ensemble averaged concentration profi les, (see Figure 6 5A ), due to reduction in Vs, the lutocline location shifts in positive y direction. Therefore, in addition to analyses of budget terms in the whole water column (see Figure 6 10), these term s should also be analyzed in the lutocline region o f case 2 and case 5. T wo separate figures a re presented to compare the TKE budget terms. The comparison s are shown between y = 15 to y = 25 ( Figure 6 11) and y = 25 to y = 30 ( Figure 6 12) which approximately correspond to the region where lutocline is observed for case 2 and case 5, respectively. The overall TKE budget terms in the boundary layer are in agreement with case 1 and case 2 (see Figure 6 10). The production and the viscous dissipation terms are the major contributing terms in the TKE budge t equation in the oscillatory boundary layer. These terms lose their intensity at the deceleration until the flow reversal and increase at the accelerating phase. S imilar to case 2, time rate of change in TKE ( ) for case 5 is signified at the instances very close to the flow reversal for case 5. If the budget terms of case 2 and case 5 are compared between y = 15 and y = 25 slightly different picture emerges. First of all, magnitudes of production, dissipation, and the bu oyancy decay terms are different. The magnitude of the buoyancy d ecay term for case 5 is 30%~ 50 % of that observed in case 2. The production and the viscous dissipati on terms for case 5 are 20 ~ 40% larger than those observed in case 2. All these observati ons show that within y = 15 to y = 25 in case 5, as opposed to ca se 2, turbulence is still large with smaller magnitude of buoyancy decay due to vertical shift in lutocline which is observed through ensemble averaged concentration profiles (see Figure 6 5A ). If moved to the lutocline location of case 5, while the magnitudes of production and viscous dissipation terms of both case 2 and case 5 are very close and small, the buoyancy dissipation for c ase 5 is at least 2 times larger than that of c ase 2. Also i t should be noted that the magnitude of buoyancy dissipation term for case 2 between y = 15 to y = 25, and PAGE 105 105 the same term for case 5 between y = 25 to y = 30 are comparable. Therefore, it can be concluded that at the lutocline the buoyancy dissipation term is intensified due to larger suspension. As stated in section 6.2, the TKE budget terms for case 6 is going to be presented for the second half of the wave cycle due to more intense shear instability between re very small due to suppressed turbulence in ca se 6. If compared to case 3 at of magnitude smaller Maximum production at of case 3 is O (106) whereas the order of production term is O (107) at case 6. On the other hand, the comparison between phase the order of magnitude of the TKE budget terms are comparable and all of them are O(107). At the flow reversal, i.e. ase 3 and = water column is quite similar, the magnitude of the production and dissipation ter ms as well as the buoyancy dissipation for case 6 are comparatively 40 % smaller than those in case 3. The wave amplitudes in C = 4.0 isocontours at the flow reversal of case 6, given in Figure 6 2B, are way smaller than the ones in case 3. This means the temporal growth rate of perturbations for case 6 is smaller than that of case 3. The difference in temporal growth rate of per turbations is observed only at the early stage of instability which might be the main cause of the difference in TKE budget terms. Nevertheless, at nd become comparable to each other again and the variation within the water column is again v ery similar. At compared to case 3 except the production In both cases there exist two peaks in the production. In case 3, the first peak from the bot tom wall has a larger value than the second peak. At this point, it should be reminded that the second peak from the bottom wall is within the vertical band where shear instability is initiated. In other words, the value of the production term at the second peak from PAGE 106 106 the wall is a measure of the generated turbulence due to shear instability. The value of the second peak at larger than that of case 3 at concluded that on average the turbulence intensity of case 6 due to shear instability is larger than the one observed in case 3. This has also been observed from the isocontour of concentration for case 6 and case 3 and the r.m.s. velocity profiles where higher turbulence level is observed after shear ins tability at the same phase. Another observation is that, in case 6 the peak point of the buoyancy dissipation term moves with the second peak of the turbulence production from the bottom wall which is quite similar to the observations made for the sa me ter m in case 3. Finally, at all budget terms except total transport term are larger compared to those at urbulence generated after the shear instability is very high and at the location of the instability buoyancy term in t he momentum equation ( Eq. 2 3) is smaller compared to case 2 due to lower local particle concentration and smaller Ri Therefore, the generated turbulence is not sufficiently absorbed by the buoyancy term in Eq. 2 3 With very close shear gradient values o bserved from Figure 6 3, it can be concluded that higher turbulence production is observed due to higher Reynolds stress which is the remnant of the shear instability. PAGE 107 107 Figure 6 1. Representation of simulation parameters in parametric space. All these simulations are conducted with constant Reynolds number Re = 1000. Closed circles indicate the simulations to identify the volume averaged concentration effect on the flow turbulence and particle trans port (Set 1) and open circles represent the simulations to identify the particle settling velocity effect on the carrier flow (Set 2). Simulation with Ri = 1 104 and Vs = 9 104 is used as the reference case for evaluation of both particle settling velocity and volume averaged concentration effects PAGE 108 108 Figure 6 2. Iso contours of concentration for case 5 and case 6. ( A) c ase 5 ( Vs = 3 104 and Ri = 1 104 ) at C = 0.4 and B) case 6 ( Vs = 27 104 and Ri = 1 104) at C = 4.0. The phase information is given in each frame. PAGE 109 109 Figure 6 3 Iso c ontours of swirl strength for case 5 and case 6. (A) c ase 5 ( Vs = 3 104 and Ri = 1 104 ) at ci 2 = (0.1~0.2) ( ci)2 max and (B) case 6 ( Vs = 27 104 and Ri = 1 104) at ci 2 = (0.01~0.02) ( ci)2 max. The phase information is given in each frame PAGE 110 110 Figure 6 4 max min) for case 6. Maximum perturbation vorticity is always positive and minimum perturbation vorticity values are always negative. Therefore, the same line legend is used for maximum and minimum vorticity. Also included in the figure is the free stream velocity of the wave with the dashed line. PAGE 111 111 Figure 6 5. Ensemble averaged velocity profiles for case 2, case 3 case 5 and case 6 The velocity profiles are sampled at six different phases in the first half of the wave cycle. Also included is the negative of the ensemble averaged velocity profile s of at six different phases spanning in the second half of the wave cycle. PAGE 112 112 Figure 6 6 En semble ave raged concentration profiles of case 2, case 3, case 5, and case 6. (A) C ase 2 ( Vs = 9 104 and Ri = 1 104 ) and case 5 ( Vs = 3 104 and Ri = 1 104 ), (B) case 3 ( Vs = 3 104 and Ri = 3 104 ) and case 6 ( Vs = 27 104 and Ri = 1 104 ) at six different phases. The phases for case 2, 3, and 5 span in the first half of the wave cycle whereas the phases of case 6 span in the second half of the wave cycle. PAGE 113 113 Figure 6 7 r.m.s of velocity fluctuation profiles for case 2 case 3 case 5 and case 6 at six different phases. The phases for case 2, 3, and 5 span in the first half of the wave cycle whereas the phases of case 6 span in the second half of the wave cycle. PAGE 114 114 Figure 6 8 Flux Richardson number and gradient Richardson number profiles for ( A) case 5 ( Vs = 3 104 and Ri = 1 104 ) and (B) case 6 ( Vs = 27 104 and Ri = 1 104 ). The phase information is included at each frame. PAGE 115 115 Figure 6 9 RSTS profiles for (A) case 5 ( Vs = 3 104 and Ri = 1 104 ) and (B) case 6 ( Vs = 27 104 and Ri = 1 104 ) at six different phases. The phase information is included at each frame. PAGE 116 116 Figure 6 10. TKE budget terms of case 5 ( Vs = 3 104 and Ri = 1 104 ). Plotted are the turbulent production (Prod), viscous dissipation (Vis. Diss.), total kinetic energy transport (Total Transp.), which is the sum of diffusion, turbulent advection, and pressure transport terms, and time rate of change in TKE ( ). The phase information is included at each frame. PAGE 117 117 Figure 6 11. Profiles of TKE budget terms for case 2 ( Vs = 9 104 and Ri = 1 104 ) and case 5 ( Vs = 3 104 and Ri = 1 104 )at y = 15~25 : Plotted are the turbulent production (Prod), viscous dissi pation (Vis. Diss.), total kinetic energy transport (Total Transp.), which is the sum of diffusion, turbulent advection, and pressure transport terms, and time rate of change in TKE ( ). The phase information is included at each frame. PAGE 118 118 Figure 6 12. Profiles of TKE budget terms for case 2 ( Vs = 9 104 and Ri = 1 104 ) and case 5 ( Vs = 3 104 and Ri = 1 104 )at y = 25~30 : Plotted are the turbulent production (Prod), viscous dissipation (Vis. Diss.), total kinetic energy transport (Tota l Transp.), which is the sum of diffusion, turbulent advection, and pressure transport terms, and time rate of change in TKE (t ). The phase information is included at each frame. PAGE 119 119 Figure 6 13. TKE budget terms of case 6. Plotted are the turbulent production (Prod), viscous dissipation (Vis. Diss.), total kinetic energy transport (Total Transp.), which is the sum of diffusion, turbulent advection, and pressure transport terms, and time rate of change in TKE ( ). The phase information is included at each frame. PAGE 120 120 CHAPTER 7 CONCLUSIONS Highly accurate numerical simulations are carried out for fine particle transport in oscillatory boundary layer. Two way coupled simulations adopting an e quilibrium Eulerian method is used to simplify particle phase formulation that is appropriate for small particle response time or St Thorough analyses are conducted to assess the domain size sufficiency and the computational resolution adequacy As a result of these analyses it is concluded that f low turbulence is fully resolved from the largest to the smallest scale at a Stokes Reynol ds number of Re = 1000. In an effort to identify the effect of particle phase concentration and particle settling velocity (or particle size ), eight simulations are conducted and the results are presented. In this study higher order particle inertia term s a re neglected in the particle velocity and the resulting turbulence modulation is only due to particle induced density stratification. The first four simulations, which are collectively denoted as set 1, are aimed at evaluating the volume averaged partic le concentration effect, quantified by the bulk Richardson number. In all the simulations of set 1, particle settling velocity is kept constant whi ch can be considered as varying volume averaged concentration for fixed particle s ettling velocity In the s econd set of the simulations, denoted as set 2, the effect o f particle settling velocity on the characteristics of particle transport and the c arrier flow turbulence is investigated. Simulation results reveal four different regimes of carrier fluid flow an d particle transport characteristics in the spectrum of both bulk Richardson number and non dimensional particle settling velocity. Regarding the effects of volume averaged concentration, the first set of simulations identifies several distinct regimes of particle laden oscillatory boundary layer according to different magnitudes of bulk Richardson number (or sediment concentration). The corresponding modes of transport and their relationship to general the fluid mud transport paradigm are PAGE 121 121 illustrated in Figure 7 1. R egime 1 is associated to the limited supply of fine sediment and hence sediment is dilute and well mixed. In this case, only dilute suspension regime can be identified. Regime 2 is associated with classic fluid mud transport paradigm shown in Figure 7 1, where the formation of a lutocline separates the dilute suspension from the more concentrated mobile fluid mud transport. Regime 3 and regime 4 are associated with highly turbulence attenuated wave boundary layer (laminarization) where the thic kness of mobile fluid mud layer is significantly reduced. In reality, it is expected that sediment particles start to deposit and consolidation of cohesive sediment bed in the concentrated aggregate network A t Ri = 104 in set 1 simulations which corresp onds to a near bed sediment concentration of 10 g/l for the Reynolds number and settling velocity selected, sediment induced stable density stratification attenuates flow turbulence and reduces mixing of sediment near the top of the turbulent wave boundary layer, which lead s to the formation of a sharp concentration gradient, i.e., lutocline. This test case also gives a qualitative estimate on the minimum concentration for which particle s cannot be considered as passive to the carrier flow turbulence. On the other hand, turbulence near the bed is not affected by the sediment and mean flow velocity is almost identical to clear fluid condition. At Ri = 304 (near bed concentration about 50 g/l), particle induced stable density stratification is strong enough to attenuate turbulence in the entire wave boundary layer. Flow tends to be laminarized for two third of the wave cycle and mean flow velocity and concentration profiles become similar to laminar solution. However, shear instability, which can be seen cle arly both in coherent vortical structure and instantaneous particle concentration, occurs during flow reversal This instability trigger s large fluctuation s in velocity and concentration and particle suspension that last about one third of the wave period. At Ri = 304 or greater, turbulence is attenuated across the entire oscillatory boundary layer PAGE 122 122 and hence the oscillatory boundary layer thickness and friction are reduced due to particle induced stable density stratification. Finally at Ri = 604 (nea r bed concentration more than 100 g/l), oscillatory boundary layer flow become completely laminarized at all times due to intense sediment induced stable density stratification. In the spectrum of settling velo city, four distinct regimes are observed. The change in local concentration certainly affects the carrier flow turbulence. The variation of settling velocity within regime 2 indicates that a decrease in settling velocity moves the location of lutocline in a counter gravitational direction and also cau ses milder concentration gradient which defines the lutocline. In addition to the change in the lutocline, it is observed that production and the viscous dissipation in the TKE budget terms ceases at the lutocline while buoyancy decay increases at the same level. With further increase in particle settling velocity, it is again observed that the laminarized flow field is followed by shear instability at the flow reversal i.e., the distinct signature of regime 3. The instability takes place at the upper limb of the theoretical laminar velocity profile and hence the order of magnitude of the term, 2e RiC at that location determines the intensity of turbulence and subsequen t suspension. In case 6, the order of magnitude is lower compared to case 3, which also falls in to regime 3. Therefore, generated turbulence level is higher compared to case 3 and the extent of the particle suspension moves up to y = 15 as opposed to y = 5 which is observed in case 3. Although there are small distinctions are present in the spectrum of Ri and Vs, the de scribed four regimes still remain Though more number of simulations are necessary to draw the borders of these regimes, to build a road ma p for future studies these regimes are presented in 2D parametric space of Ri and Vs including a rough estimate of regions where these regimes exist for given flow conditions, i.e. Re =1000 (see Figure 6 2) The borders lying between each PAGE 123 123 regime should ideally be determined by hydrodynamic stability analysis which is not within the scope of the thesis. However, still a thorough investigation is necessary for theoretical instability analysis. Even the turbulent clear fluid flow in oscillatory chan nel is not fully understood and the discrepancies between the theoretical and experimental studies are not yet identified. It should also be reminded that Schmidt number, Sc used in this study is 0.5 to mimic the random subgrid scale particle motion. How ever, the value of Sc is still an open research question and also the borders of these regimes may be Schmidt number dependant. Therefore, although the present numerical study forms a useful plot to begin to understand the state of the muddy seabed the re gions of each regime in 2D parametric map should be improved by theoretical, field and laboratory experiments. Apart from, Sc the existence of these regimes is also dependant on the intensity of carrier flow turbulence. The competition between the turbule nce generation and the suppression capacity of the particles within a given Ri, Vs, and Sc determines the regimes described. Hence, Re is a dictating parameter that determines the extent of the regimes described. Therefore, the variation of regimes withi n the spectrum of Re is an other important step that will widen the understanding on the limits of the se regimes and their variability. A set of simulations is under progress as a next step to the current study and the illustration of the set of simulations is given in Figure 7 3. In the current study, the particles are assumed to be of small St Therefore, higher order inertial terms are neglected With the increase in the particle size, the response of the carrier flow turbulence, particle suspension and deposition shall not be predicted in desired accuracy by only considering the particle induced stable density stratification. Therefore, as a next step to the PAGE 124 124 current study, simulations with higher order particle inertia terms in the particle velocity sha ll improve the prediction for larger particles. The importance of regime 2 comes from the high suspension of particle s In the presence of mild downslope gravitational force, it is suspected by the author that most wave supported gravity driven mud flows t ake place in this regime which should be investigated as a future study. In addition background conditions in lutocline fo rmation and the prediction of its location are point s of interest not only to predict the fate of sediment particles but also practi cal applications such as backscatter of the sound waves due to suspended sediment near the bed This is very critical to detect the submarine threats via sonar and t herefore, special attention should be diverted to the background conditions of the formatio n of the lutocline and its location from the bed. Another important phenomenon is the existence of inter particle stresses which is not considered in the current study. In most fine sediment transport studies, when sediment concentration is greater than O( 100) g/l, rheological stresses caused by interactions among the floc aggregates and interstitial water need to be considered. Rheological stresses give rise to enhanced effective viscosity that can further increase energy dissipation and friction factor. I t is the authors belief that inclusion of these nonNewtonian flow features shall form an important step forward to evaluate the effect of rheology in an effort to understand the fate of terrestrial sediment in the coastal environment and hydrodynamic dis sipations PAGE 125 125 Figure 7 1. Illustrative sketch that shows the of fluid mud transport profile Figure 7 2 The location of the simulations in two dimensional parametric space of Vs and Ri. The dashed line shows the estimate of the borders that defines each of these regimes. Filled circles show the first set of simulations and the diamonds show the second set of simulations. PAGE 126 126 Figure 7 3 A complete set of simulations These simulations en velope a cubic region in three dimensional parametric space PAGE 127 127 LIST OF REFERENCES AKHAVAN, R., KAMM, R. D., & SHAPIRO, A. H. 1991a An investigation of transition to turbulence in bounded oscillatory Stokes flows. Part 1. Experiments. J. Fluid Mech 225, 395. AKHAVAN, R., KAMM, R. D., & SHAPIRO, A. H.1991b An investigation of transition to turbulence in bounded oscillatory Stokes flo ws. Part 2. Numerical simulations. J. Fluid Mech 225, 423. ALISHAHI, M. R. & KRONE, R. B.,1964 Suspension of cohesive sediment by windgenerated waves. Technical report HEL 2 9, Hydraulic Engineering Laboratory, University of California, Berkeley, Califor nia, August 24. ARAKAWA, A., 1966 Computational design for longterm numerical integration of the equations of fluid motion: Two dimensional incompressible flow. part 1. J. Comput. Physics 1 119143. BAGCHI, P. & BALACHANDAR, S. 2003 Effect of turbulenc e on the drag and lift of a particle. Phys. Fluids 15, 3496. BAGCHI, P. & BALACHANDAR, S. 2004 Response of the wake of an isolated particle to an isotropic turbulent flow. J. Fluid Mech 518,95. BALACHANDAR, S. & EATON, J.K. 2010 Turb ulent dispersed multiphase flow. Ann. Rev. Fluid Mech BARENBLATT, G. F., 1953 On the motion of suspended particles in a turbulent stream. Prikladnaja Mat. Mekh., Engl. Trans. 17, 261. BATES, C. C., 1953 Rational theory of delta formation. American Association of Petrol ium Geologists Bulletin 27, 21192216. BEST, J.L., 2005 The fluid dynamics of river dunes: a review and some future research directions. J. Geophys. Re s. Earth Surface, 110, F04S02, doi:10.1029/2004JF000218. BLENNERHASSET, P. J. & BOSSOM, A. P., 2006 The linear stability of highfrequency oscillatory flow in a channel. J. Fluid Mech 556, 1 25. BONOMETTI, T. & BALACHANDAR, S. 2008 Effect of Schmidt number on the structure and propagation of density currents Theor. Comput. Fluid Dyn. 22, 341 361 DOI 10.1007/s0016200800852 BURTON, T M & EATON J.K. 2005 Fully resolved simulations of particle turbulence interaction. J. Fluid Mech 545,67. CANTERO, M. I., BALACHANDAR, S. & M. GARCIA 2008 An Eulerian Eulerian model for gravity currents dri ven by inertial particles, Int. J. Multiphase Flow 34, 484501. PAGE 128 128 CANTERO, M.I., S. BALACHANDAR, CANTELLI, A., PIRMEZ, C., & PARKER, G. 2009. Turbidity current with a roof: Direct numerical simulation of self stratified turbulent channel flow driven by suspended sediment J Geophys. Res 114, C03008. CANUTO, C., HUSSAINI, M. Y., QUARTERONI, A., & ZANG, T. A., 1987 Spectral Methods in Fluid Dynam ics Spr inger Verlag CHAKRABORTY, P., BALACHANDAR, S. & ADRIAN, R.J. 2005 On the relationships between local vortex identification schemes. J. Fluid Mech 535:189214 CHEN, D., CHEN, C., TANG, F., STANSBY, P., & LI, M. 2007 Boundary layer structure of oscillatory open channel shallow flows over smooth and rough beds COSTAMAGNA, P., VITTORI, G., & BLONDEAUX, P. 2003 Coherent structures in oscillatory boundary layers. J. Fluid Mech 474, 1. DALRYMPLE, R. A., & LIU, P. L., 1978 Waves over soft mud beds: a two layer fluid mud model. J. Phys. Oceanogr 8 11211131. EINSTEIN, H. A. & CHIEN, N. 1955 Effects of heavy sediment concentration near the bed on velocity and sediment distribution, MRD Sediment Ser. 8, Univ. of Calif., Berkeley. ELGAR, S. & RAUBENHEIMER, B., 2008 Wave dissipation by muddy seafloors, Geophys. Res. Lett. 35 (L07611), doi:10.1029/2008GL033245 ELGOBASHI, S., 1991 Particle laden turbulent flows: direct simulation and closure models. Appl. Sci. Res. 52, 301. ELGOBASHI, S. & TRUESDELL, G.C. 1992, Direct simulation of particle dispersion in a decaying isotropic turbulence, J. Fluid Mech. 242, 655. ELGOBASHI, S. & TRUESDELL G.C. 1993. On the 2 way interaction between homogeneous turbulence and dispersed solid pa rticles 1.turbulence modification, Phys. Fluids A, 5, 1790. FERRY, J. & BALACHANDAR, S., 2001 A fast Eulerian method for disperse two phase flow, Int. J. Multiphase Flow 27 1199. FERZIGER, J. H., & P, M. 2002. Computational methods for fluid dynamics Springer Verlag. FORRISTALL, G. Z. & REECE, A. M. 1985 Measurements of wave attenuation due to a soft bottom: The SWAMP experiment, J. Geophys. Res. 90, 3367 3380 FOSTER, D. L., BEACH, R. A., & HOLMAN, R.A., 2006 Turbulence observations of the nearshore wave bottom boundary layer. J. Geophy. Res. 111, 1 11. PAGE 129 129 GELFENBAUM, G. & SMITH, J. D. 1986 Experimental evaluation of a generalized suspended sediment transport theory, in Shelf Sands and Sand Stones, edited by R. J. Knight and J. R. McLean, 133144, Can. Soc. of Pet. Eng., Calgary, Alberta GLENN, S. M., & GRANT, W. D. 1987 A suspended sediment stratification correction for combined wave and current flows. J. Geophys. Res. 92(C8), 8244 8264. GOLDSMI TH, S. T., CAREY, A. E., LYONS, W. B., KAO, S. J., LEE, Y. T., & CHEN, J. 2008 Extreme storm events, landscape denudation, and carbon sequestration: Typhoon Mindulle, Choshui River, Taiwan. Geology 36(6), 483486. HALL, P. 1978 The linear stability of flat Stokes layers. Proc. R. Soc. Lond. A 359, 151166. HARRIS C.K., TRAYKOVSKI, P.A., & GEYER, W. R. 2005 Flood dispersal and deposition by near bed gravitational sediment flows and oceanographic transport: a numerical m odeling study of the Eel River shelf, northern California. J Geophys. Res 110(C9), (1 16). HILL P.S., T.G., MILLIGAN, & W.R. GEYER. 2000 Controls on effective settling velocity in the Eel River flood plume Continental Shelf Research, 20: 20952111. HINO, M. 1963 Turbulen t flow with suspended particles. J. Hydraul. Div. Am. Soc. C iv. Eng 89(HY4), 161185. HINO, M., SAWAMOTO, M. & TAKASU, S. 1976 Experiments on transition to turbule nce in an oscillatory pipe flow. J. Fluid Mech 75, 193. HINO M., KASHIW AYANAGI, M., NAKAYAMA, A. & HARA, T., 1983 Experiments on the turbulence statistics and the structure of a reciprocating oscillatory f low J. Fluid Mech., 131, 363 400. HOWARD, L. N., 1961 Note on a paper by John W. Miles. J. Fluid Mech 10 509. HSU, T. J. & HANES, D.M., 2004 The effects of wave shape on sheet flow sediment transport. J. Geophys. Res. 109(C5), C05025, doi:10.1029/2003JC002075. HSU, T J, TRAYKOVSKI, P. A., & KINEKE, G. C., 2007 On modeling boundary layer and gravity driven fluid mud tr ansport, J. Geophys. Res. 112, C04011, doi:10.1029/2006JC003719. HSU, T J., OZDEMIR, C. E., & TRAYKOVSKI, P. A. 2009 High resolution numerical modeling of wave supported sedim ent gravity driven mudflows J. Geophys. Res ., 114, C05014, doi:10.1029/2008JC005006. HUPPERT, H. E., TURNER, J. S., & HALLWORTH, M. A. 1995 Sedimentation and entrainment in dense layers of suspended particles stirred by an oscillating grid. J. Fluid Mech ., 289, 263 293. PAGE 130 130 JACKSON, J. R., 1973 A model study of the e ffects of small amplitude waves on the resuspension of fine grain cohesive sediments. M. S. Thesis, University of New Hempshire, Durham, New Hempshire, 53. JENSEN, B. L., SUMER, B. M., & FREDS J., 1989 Turbulent oscillatory boundary layers at high Reyno lds numbers. J. Fluid Mech 206, 265. LAMB, M. P., DASARO, E. & PARSONS, J. D. 2004 Turbulent structure of highdensity suspensions formed under waves. J. Geophys. Res. 109, C12026. LARSON, M. & KRAUS, N. C., 1995 Prediction of cross shore sediment transport at differe nt spatial and temporal scales, Marine Geology 126 111127. MAA, J.P. Y. & MEHTA, A.J. 1987 Mud erosion by waves: A laboratory study Cont. Shelf Res 7(11/12), 12691284. MAXEY, M. R. 1987 The gravitational settling of aerosol partic les in homogeneous turbulence and random flow fields, J. Fluid Mech 174, 441. MCLAUGHLIN, J.B., 1989 Aerosol particle deposition in numerically simulated channel flow, Phys. Fluids 1 12111224. MEHTA, A. J. 1987 On estuarine cohesive sediment suspension behavior J Geophys. Res 94(C10):1430314314 MEI, C. C., & LIU, K. F., 1987 A Bingham plastic model for a muddy seabed under long waves J. Geophys. Res 94(C13) 14,58114,594. MILES, J. W. 1961 On the stability of heterogeneous shear flows. J. Fluid Mech 10, 496. MULDER, T & SYVITSKI, J. P. M. 1995 Turbidity currents generated at river mouths during exceptional discharges to the world oceans, J. Geol ., 103, 285 299. NOH, Y. & FERNANDO, J. S., 1991 Dispersion of suspended particles in turbulent flow, Physics of Fluids A, 3 (7), 17301740. OZDEMIR, C. E., HSU, T J, & BALACHANDAR, S. 2010 Simulation of fine sediment transport in oscillatory boundary layer, Journal of Hydro environment Research, 4, doi: 10.1016/j.jher.2009.10.013. PEYRET, R. 2002. Spectral methods for incompressible viscous flow Applied mathematical sciences 148, Springer Verlag. POPE, S. 2000 Turbulent flows Cambridge Univ Press. ROGERS, C. B. & EATON, J.K. 1991 The effect of small pa rticles on fluid turbulence in a flat plate, turbulent boundary layer in air. Phys. Fluids A 3 928 37. PAGE 131 131 ROSS, M. A.& MEHTA, A. J 1989 On the mechanics of lutoclines and fluid mud, J. Coast. Res., Spec. Issue 5 51 61. SALON, S., ARMENIO, V., & CRISE, A., 2007 A numerical investigation of the Stokes boundary layer in the turbulent regime, J. Fluid Mech 570, 253. SARPKAYA, T. 1993 Coherent structures in oscillatory boundary layers. J. Fluid Mech. 253, 105. SCULLY, M.E., FRIEDRICHS, C.T., & WRIGHT, L.D., 200 3 Numerical modeling of gravity driven sediment transport an deposition on an energetic continental shelf: Eel Riv, Northern California. J. Geophys. Res 108 (C4), 17 1 1714). SHOTORBAN, B. & BALACH A NDAR, S. 2006. Particle concentration in homogeneous shear turbulence simulated via Lagrangian and equilibrium Eulerian approaches, Phys. of Fluids 18, 065105. SPALART, P. R. & BALDWIN, B. S. 1989 Direct simulation of a turbulent oscillating boundary layer. Tu rbulent Shear Flows 6 Springer. SQUIRES, K. D. & EATON, J. K. 1991 Preferential concentration of particles by turbulence, Phys. Fluids A 3 1169. STERNBERG, R. W., CACCHOINE, D. A., PAULSON, B., KINEKE, G. C. & DRAKE D. E. 1996 Observations of sediment transport on the Amazon subaqueous delta. Continental Shelf Research 16, No. 5/6, 697715. STYLES, R. & GLENN, S. M., 2000 Modeling stratified wave and current bottom boundary layers in the continental shelf, J. Geophys. Res. 105 2411924139. TANAKA, T., & EATON, J .K. 2008 Classification of turbulence modification by dispersed spheres using a novel dimensionless number. Phys. Rev. Lett 101, 114502. TAYLOR, 1931 G.I., Effect of variation in density on the stability of superposed streams of fluid. In: Proc. R. Soc. London, Series A 132 (1931), p. 499. TRAYKOVSKI, P., GEYER, W. R., IRISH, J. D., & LYNCH, J. F. 2000 The role of wave induced fluid mud flows for cross shelf transport on the Eel river continental shelf. Cont. Shelf. Res ., 20, 21132140. TRAYKOVSKI, P., WIEBERG, P.L., & GEYER, W. 2007 Observations and modeling of wave supported sediment gravity flows on the Po prodelta and comparison to prior observations from the Eel shelf Cont. Shelf. Res., 27, 375399. TROWBRIDGE, J. H., & KINEKE, G. C., 1994 Structure and dynamics of fluid mud on the Amazon continental shelf, J. Geophys. Res., 99(C1), 865 874. TURNER, J. S. 1973 Buoyancy Effects in Fluids, Cambridge Univ. Press,New York. PAGE 132 132 United Nations Atlas of the Oceans, http:/ /www.thew2o.net/atlas Contents Copyright 2009. VANONI, V. A. 1946 Transportation of suspended sediment by water ASCE Trans ., 111, 67 133. VANRIJN, L. C., 1981 Entrainment of fine sediment particles; development of concentration profiles in a steady, u niform flow without initial sediment load. Rep. No M1531, Part II, Delft Hydraulic Lab., Delft, The Netherlands. VITTORI, G., 2003 Sediment suspension due to waves. J. Geophys. Res 108 (C6), 3173. VITTORI, G ., & VERZICCO R., 1998 Direct simulation of transition in an oscillatory boundary layer. J. Fluid Mech 371:207 232. VON KERCZEK, C. & DAVIS, S. H.1974 Linear stability theory of oscillatory Stokes layers. J. Fluid Mech. 62, 753773. WINTERWERP, J. C., 2001 Stratification effec ts by cohesive and noncohesive sediment. J. Geophys. Res 106 (C10), 22,55922,574. WINTERWERP, J. C., 2002 On the flocculation and settling velocity of estuarine mud. Continental Shelf Res, 22:1339 1360 WINTERWERP, J. C., & VAN KESTEREN, W. G. M., 2004 Introductio n to the physics of cohesive sediment in the marine environment, Elsevier WRIGHT, L. D., & C.A. NITTROUER, 1995 Dispersal of river sediments in coastal seas: six contrasting cases, Estuaries, 18, 494508. ZENG, L ., BALACHANDAR, S., FISCHER, P ., & NAJJAR, F .M. 2008 Interactions of a stationary finite sized particle with wall turbulence. J. Fluid Mech 594:271 305 ZENG, L., BALACHANDAR, & NAJJAR, F M. 2009 Wake response of a stationary finite sized particle in a turbulent channel flow. Int. J Multiphase Flow Under consideration ZHOU, J., R., ADRIAN, J. & BALACHANDAR, S & KENDALL, T. M. 1999 Mechanisms for generating coherent packets of hairpin vortices in channel flow. J. Fluid Mech 387, 353 396. ZHOU, D. & NI, J. R. 1995 Effects of dynamic interaction on sediment laden turbulent flows, J. Geophys. Res. 1 00(C1), 981 996. PAGE 133 BIOGRAPHICAL SKETCH Celalettin Emre Ozdemir was born as the only son of Cumhur Ozdemir, who was a social worker and Zeliha Ozdemir, who is currently a retiree from being high school teacher, and the brother of Dr. Nilufer Aslihan Ozdemir who is right now pursuing an academic career as post doctoral research associate in University Catholique Louvein. He graduated from Middle East Technical University wi th B.S. degree in civil engineering in 2001 and got his M.S. degree from the same department in 2003. After two years of educational experience at IIHR Hydroscience and Engineering, he found the right place to make quality research He would like to pursue a career in research He believes high quality research comes from creative thinking supported by deep background, discipline, experience, and a working environment free of discrimination based on prejudices, corruption, and tyranny. 