Course Hero - We put you ahead of the curve!

You have requested the below document.

Coursehero >> Texas >> Texas >> BOUGIEJL 81498

by Copyright Jonathan Lee Bougie 2004 The Dissertation Committee for Jonathan Lee Bougie certifies that this is the approved version of the following dissertation: Continuum Simulations of Fluidized Granular Materials Committee: Jack B. Swift, Supervisor Graham F. Carey Richard Fitzpatrick William D. McCormick Harry L. Swinney Continuum Simulations of Fluidized Granular Materials by Jonathan Lee Bougie, B.A. DISSERTATION Presented to the Faculty of the Graduate School of The University of Texas at Austin in Partial Fulfillment of the Requirements for the Degree of DOCTOR OF PHILOSOPHY THE UNIVERSITY OF TEXAS AT AUSTIN August 2004 This dissertation is dedicated to my wife Julie. Acknowledgments This work would not have been possible if it weren't for the help, guidance, and friendship of many people. My advisor, Jack Swift, has been an invaluable source of guidance, patience, and support. He has been a great person to work with and learn from during my years in Austin. I would also like to thank Harry Swinney for all that he has contributed to this research and to my education and development as a scientist. In addition, Bill McCormick has provided great insight and discussion which cannot be overlooked. Many thanks to all those at the Center for Nonlinear Dynamics, past and present, who have been instrumental through their work, their insights, and their questions. These people have made the CNLD a stimulating and rewarding place to work. People whose discussion has been vital to my graduate career include Sung Joon Moon, Erin Rericha, Dan Goldman, Jennifer Kreft, Andy Lee, Matthias Schro ter, Mark Shattuck, and many more. Everyone at e the CNLD deserves mention both for creating a space for exciting research as well as for all the friendships over the years. Outside of the CNLD, my family has been a constant source of support. Thanks go out especially to my parents, Michael and Bonnie, and to my brother Devin. My wife Julie has always been there for me and deserves special thanks. v Finally, I want to thank all my friends, colleagues, and comrades. There are way too many to list here, but thanks go out to you all. This work was supported by the Engineering Research Program of the Office of Basic Energy Sciences of the Department of Energy (Grant DE-FG0393ER14312) and a grant from the Texas Advanced Research Program (Grant ARP 003658-0055-2001). vi Continuum Simulations of Fluidized Granular Materials Publication No. Jonathan Lee Bougie, Ph.D. The University of Texas at Austin, 2004 Supervisor: Jack B. Swift A successful hydrodynamic theory of granular media could allow scientists and engineers to exploit the powerful techniques of fluid dynamics to describe granular phenomena. We use computer simulations to test a set of continuum hydrodynamic equations for granular media which were proposed by Jenkins and Richman nearly twenty years ago [59]. We use these continuum simulations as well as molecular dynamics (MD) simulations to investigate phenomena in vertically shaken layers of grains. When a layer of grains is shaken vertically by a plate with maximum acceleration greater than the acceleration of gravity, the layer will be thrown off the plate with each cycle of the plate. Continuum and MD simulations show that normal shocks form in the layer upon contact with the plate later in the cycle. We show that increasing the coefficient of restitution of the particles increases the speed of the shock in the layer, and that the limit of perfectly elastic particles is not singular. In addition, deeper layers of particles exhibit denser vii packing fractions near the plate and higher shock speeds than shallow layers. Pressure gradients produced by these shocks play an important role in the dynamics of standing wave patterns formed in oscillated granular layers. Both continuum and molecular dynamics simulations produce standing waves with wavelengths which agree with previous experiments. Continuum simulations reproduce stripe patterns found in MD simulations of frictionless particles, but do not reproduce square or hexagonal patterns found in experiments and MD simulations of frictional particles. Finally, we show that fluctuations present in molecular dynamics simulations are not captured by our current continuum model. By fit to Swift-Hohenberg theory, we find these fluctuations in MD simulations to be several orders of magnitude larger than fluctuations found in ordinary fluids. Differences between patterns in continuum and MD simulations near onset are found to be consistent with the presence of fluctuations in MD simulations without friction. viii Table of Contents Acknowledgments Abstract List of Tables List of Figures Chapter 1. Introduction 1.1 Modeling of Granular Flows . . . . . . . 1.1.1 Molecular Dynamics Simulation . . 1.1.2 Continuum Equations . . . . . . . 1.2 Vertically Oscillated Granular Layers . . 1.3 Patterns in Vertically Oscillated Granular 1.4 The Role of Friction in Shaker Patterns . 1.5 Shocks in Granular Materials . . . . . . . 1.6 Outline . . . . . . . . . . . . . . . . . . . Chapter 2. v vii xiii xiv 1 3 3 5 8 9 12 15 16 . . . . . . . . . . . . . . . . Layers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Previous Solutions to Granular Continuum Models 19 23 23 27 28 28 29 30 31 31 Chapter 3. Simulation Methods 3.1 A Continuum Model of Granular Materials . . 3.2 A Continuum Simulation of Granular Shakers 3.2.1 Continuum Fields . . . . . . . . . . . . 3.2.2 User input . . . . . . . . . . . . . . . . 3.2.3 Initial conditions . . . . . . . . . . . . . 3.2.4 Algorithm . . . . . . . . . . . . . . . . . 3.2.5 Boundary Conditions . . . . . . . . . . 3.2.6 Spatial Derivatives . . . . . . . . . . . . ix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.7 Size of Spatial Grid . . . . . . . . . . . . . . . . 3.2.8 Forward timestepping . . . . . . . . . . . . . . . 3.2.9 Artificial Dissipation in Low Density Regions . . 3.3 A Molecular Dynamics Simulation of Granular Shakers . . . . . . . . . . . . . . . . 33 39 41 44 47 47 49 50 50 51 52 52 54 54 54 56 57 57 59 64 67 71 73 73 75 75 75 76 77 77 78 Chapter 4. Shocks in Shaken Granular Layers 4.1 Background on Shock Waves . . . . . . . . . . . . . . . . . . . 4.2 Properties of the Shaken Layer . . . . . . . . . . . . . . . . . . 4.3 Simulation Parameters . . . . . . . . . . . . . . . . . . . . . . 4.3.1 Molecular Dynamics Simulations . . . . . . . . . . . . . 4.3.2 Continuum Simulations . . . . . . . . . . . . . . . . . . 4.4 Shock Formation and Propagation . . . . . . . . . . . . . . . . 4.4.1 Tracking the Layer, the Shock, and the Shock Width . . 4.4.2 Shock Profile During Formation and Propagation . . . . 4.4.2.1 f t = 0: The layer falls towards the plate . . . . 4.4.2.2 f t = 0.15: Shock forms as the layer strikes plate 4.4.2.3 f t = 0.22: Shock propagates through the layer . 4.4.2.4 f t = 0.40: The layer begins to leave the plate . 4.4.3 Shock Dynamics . . . . . . . . . . . . . . . . . . . . . . 4.5 Mach Number and Shock Width . . . . . . . . . . . . . . . . . 4.6 Effect of Inelasticity . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Shocks in Layers of Different Depths . . . . . . . . . . . . . . . 4.8 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 5. Patterns in Shaken Granular Layers 5.1 Patterns Near Onset in Shaken Granular Layers . . . . . . . . 5.2 Properties of the shaken layer . . . . . . . . . . . . . . . . . . 5.3 Simulation Parameters . . . . . . . . . . . . . . . . . . . . . . 5.3.1 Molecular Dynamics Simulations . . . . . . . . . . . . . 5.3.2 Continuum Simulations . . . . . . . . . . . . . . . . . . 5.3.3 Visualizing Patterns . . . . . . . . . . . . . . . . . . . . 5.4 Standing Wave Stripe Patterns in Simulations of Granular Layers 5.5 Dispersion Relations in Continuum, MD, and Experiment . . . x 5.6 Dynamics of Standing Wave Patterns . . . . . . . . . . . . . . 5.6.1 Wave Patterns in MD Simulations . . . . . . . . . . . . 5.6.1.1 f t = 0: The layer collides with the plate . . . . 5.6.1.2 f t = 0.25: The layer is compressed on the plate 5.6.1.3 f t = 0.5: The layer begins to leave the plate . . 5.6.1.4 f t = 0.75: The layer begins to leave the plate . 5.6.1.5 1.0 f t < 2.0: The cycle repeats . . . . . . . . 5.6.2 Wave Patterns in Continuum Simulations . . . . . . . . 5.6.2.1 f t = 0: The layer collides with the plate . . . . 5.6.2.2 f t = 0.25: The layer collides with the plate . . . 5.6.2.3 f t = 0.5: The layer begins to leave the plate . . 5.6.2.4 f t = 0.75: The layer is leaving the plate . . . . 5.6.2.5 1.0 f t < 2.0: The cycle repeats . . . . . . . . 5.7 Shocks and Patterns . . . . . . . . . . . . . . . . . . . . . . . . 5.7.1 Shock Formation and Propagation . . . . . . . . . . . . 5.7.1.1 f t = 0.875: The layer is off the plate . . . . . . 5.7.1.2 f t = 1: The layer collides with the plate . . . . 5.7.1.3 f t = 1.125: The shock propagates through the layer . . . . . . . . . . . . . . . . . . . . . . . . 5.7.1.4 f t = 1.25: The shock leaves the layer . . . . . . 5.8 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 6. 81 83 83 85 85 86 86 86 87 91 91 92 93 93 94 96 96 97 101 101 Onset of Standing Wave Patterns and the Role of Fluctuations 103 6.1 Onset of Patterns in Shaken Granular Layers . . . . . . . . . . 103 6.2 Simulation Parameters . . . . . . . . . . . . . . . . . . . . . . 104 6.2.1 Molecular Dynamics Simulations . . . . . . . . . . . . . 104 6.2.2 Continuum Simulations . . . . . . . . . . . . . . . . . . 105 6.3 Onset of Patterns for Increasing . . . . . . . . . . . . . . . . 105 6.3.1 Layers Above and Below Onset of Patterns . . . . . . . 105 6.3.2 Growth and Ordering of Patterns With Increasing . . 106 6.3.3 Characterizing patterns . . . . . . . . . . . . . . . . . . 108 6.3.3.1 Onset of patterns in continuum simulations . . 109 xi 6.4 6.5 6.6 6.7 6.3.3.2 Onset of patterns in molecular dynamics simulations . . . . . . . . . . . . . . . . . . . . . . . . 110 Fluctuating Hydrodynamics . . . . . . . . . . . . . . . . . . . 113 6.4.1 Swift-Hohenberg Simulations . . . . . . . . . . . . . . . 114 6.4.2 Comparing Swift-Hohenberg and molecular dynamics simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 The Effect of Noise on Onset . . . . . . . . . . . . . . . . . . . 118 The Effect of Changing Frequency on Noise Strength . . . . . 119 6.6.1 Growth of Patterns and Fluctuations for Varying Frequency119 6.6.2 Noise Strength Variation with Frequency . . . . . . . . 122 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 129 129 130 133 135 136 136 138 139 140 Chapter 7. Future Work: Modifying the Continuum Model 7.1 Granular Fluctuating Hydrodynamics . . . . . . . . . . . . . 7.1.1 Further study of noise in the granular shaker . . . . . 7.1.2 Investigating noise in other systems . . . . . . . . . . 7.2 Adding Friction to the Continuum Model . . . . . . . . . . . 7.2.1 Continuum Models with Friction . . . . . . . . . . . . 7.2.1.1 Including Rotations in Continuum Model . . . 7.2.1.2 Additional approximations . . . . . . . . . . . 7.2.1.3 A possible method for modeling friction . . . 7.2.2 Proposed procedure for incorporating friction . . . . . . . . . . . . . . Chapter 8. Conclusions and Future Work 142 8.1 Discussion of major results . . . . . . . . . . . . . . . . . . . . 143 Bibliography Vita 150 166 xii List of Tables 6.1 6.2 Coefficient of RB noise strengths FRB at the densest region of the cell for f = 0.25 (FRB (0.25, f t)), and for f = 0.4174 (FRB (0.4174, f t) are tabulated for various times f t in the cycle (MD simulation). Also shown is the wavelength in units of particle diameter ((f )/) for each frequency. . . . . . . . . The ratio of the RB noise strength for f2 = 0.4174 to the RB noise strength for f1 = 0.25 at the densest region of the cell q0 (f2 )FRB (f2 , f t)/q0 (f1 )FRB (f1 , f t) is tabulated for various times f t in the cell. Shown separately are the ratio of nondimensional wavenumbers q0 (f2 )/q0 (f1 ), and of the RB de pendence on n, u, and T : FRB (f2 , f t)/FRB (f1 , f t). . . . . . . 124 125 xiii List of Figures 1.1 1.2 1.3 Electromagnetic shaker system (a) for visualizing patterns in oscillated granular layers. The container is illuminated by light incident at low angles from the side (arrows in (b)). The light reflects off the particles and is visualized by a CCD camera, yielding bright spots for peaks and dark spots for valleys of the pattern. From [113]. . . . . . . . . . . . . . . . . . . . . . . . Stripe patterns (a), square patterns (b), and hexagonal patterns (c) formed in vertically oscillated granular layers. All patterns oscillate at f /2, where f is the frequency of oscillation of the plate; in the hexagon pattern, the right and left sides of the cell are out of phase with respect to each other (peaks in the center of the hexagon on the right, peaks at the border on the right). The shaking parameters , f , and H are given by (a) 3.3, 67 Hz, 7, (b)2.9, 25 Hz, 4, and (c) 4.0, 67 Hz, 7, respectively. In (a) and (c), the diameter of the container is 770 and in (b), the container is square with sides 1100, where the particles used are bronze spheres of diameter = 0.165 mm. (a) and (c) are from [76], (b) is from [39], and the entire figure was composed in [109]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Phase diagram for a vertically oscillated granular layer of bronze spheres of diameter = 0.165 mm, where the layer depth H = 5, and the container diameter is 770. The diagram indicates the granular patterns observed as a function of the dimensionless accelerational amplitude and dimensionless frequency of the plate oscillation f . The transition from a flat layer to squares is hysteretic; solid lines indicate the transition for increasing , dotted lines indicate the transition for decreasing . Similar patterns are observed for a variety of particle types and layer depths. From [79]. . . . . . . . . . . . . . . . . . . . . . 10 11 13 xiv 1.4 1.5 Defect creation and "melting" of a square pattern in a shaken granular layer for (a) experiment and (b) MD simulation for = 0.17 bronze spheres oscillated at = 2.9 and f = 32 Hz. In the experiment, graphite has been added to reduce friction, and at t = 0, a frequency modulation is added to experiment at fL = 2 Hz. In MD simulation, no frequency modulation is added, but at time t = 0, friction was turned off in the simulations. The numbers represent the number of oscillations since the chagne at t = 0. The insets in each picture show the power spectrum (white is low power, black is high power) corresponding to each image. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Horizontal component of the velocity field (vx ) of a flow of stainless steel spheres with diameter = 1.2mm falling under gravity onto a wedge as determined by (a) experiment, (b) MD simulations, and (c) continuum equations to Navier-Stokes order. Each picture shows a region of size 130 104. A shock separates an undisturbed region in which the horizontal component of flow is nearly zero, from a compressed region in which streamlines (solid lines with arrows) are diverted around the obstacle. From [93]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . Center of mass height zcm (x, y, t = 0) as a function of horizontal location x, y for a grid with a number of horizontal gridpoints (a) Nx = Ny = 42, (b) Nx = Ny = 28, (c) Nx = Ny = 21, (d) Nx = Ny = 14, (e) Nx = Ny = 11, and (f) Nx = Ny = 6. These images show the initial conditions for the convergence test, based on linear interpolation between gridpoints from the oscillatory state of the Nx = Ny = 42 simulation (a). . . . . . Center of mass height zcm (x, y, t = 10/f ) as a function of horizontal location x, y for a grid with a number of horizontal gridpoints (a) Nx = Ny = 42, (b) Nx = Ny = 28, (c) Nx = Ny = 21, (d) Nx = Ny = 14, (e) Nx = Ny = 11, and (f) Nx = Ny = 6. These images are snapshots of the layer taken 10 cycles of the plate after the images in Fig. 3.1. . . . . . . . . . . . . . . . . Percent difference from the reference simulation of Nref = 1764 as a function of changing number of gridpoints N = Nx Ny , shown on a linear scale (a) and a log-log scale (b). . . . . 14 16 3.1 35 3.2 36 3.3 38 xv 3.4 3.5 Functional form of the strength of artificial dissipation f () as a function of volume fraction , a) for the range 0 max , and b) for the range 0 0.1max , where max = 0.65 is the random close-packed volume fraction in 3D. The artificial viscosity is only nonnegligible for volume fractions which are very small compared to the random close-packed volume fraction. Thus, this artificial viscosity has a negligible effect on the shocks and patterns found in the system, while playing an important role in ensuring numerical stability in extremely low density regions. . . . . . . . . . . . . . . . . . . . . . . . . . . 43 Standing wave patterns from experiments using lead spheres and from molecular dynamics simulations for different nondimensional frequencies f and accelerational amplitude relative to gravity ; showing (a) square patterns (b) stripes, (c) and (d) alternating phases of hexagons, (e) a flat layer, (f) squares, (g) stripes, and (h) hexagons. Patterns (a)-(e) oscillate at f /2, while patterns at higher (f)-(h) oscillate at f /4. From Ref. [7]. 46 Dimensionless temperature T /(g) and volume fraction as functions of dimensionless height z/ at four times f t in the cycle. For each time, a picture from MD simulation is shown in the left column, with particles color-coded according to temperature: high T in red, low T in blue, and the bottom plate of the container shaded solid gray. The right column shows horizontally averaged (blue, labeled on the bottom graph) and T /(g) (red, labeled on the top graph) as functions of z/ (ordinate) for the same four times. The plate is shown as a horizontal black solid line, results from MD simulation are dots, and continuum results are solid lines. At times when a shock is present in the graph, the shock region (see Sec. 4.4.1) is shaded in gray. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Location of the shock (solid lines for continuum, triangles for MD) and the center of mass of the layer (dashed lines for continuum, circles for MD) as a function of time f t during one cycle of the plate (thick solid line) for particles with (a) e = 0.99, (b) e = 0.90, and (c) e = 0.85, starting from the same initial conditions at f t = -0.33. The plot is shaded according to the volume fraction from the continuum simulation, so that high volume fraction is dark and low volume fraction is light. The "top" and the "bottom" of the layer from MD (when the volume fraction drops to less than 4% of random close-packed) are shown as +'s and 's, respectively. . . . . . . . . . . . . . . . 4.1 55 4.2 58 xvi 4.3 4.4 4.5 4.6 4.7 4.8 Mach number (M a) (dashed line for continuum and triangles for MD) and rescaled volume fraction (10) (solid line for continuum and circles for MD) as functions of z/ when f t = 0.22. At this time the shock is fully developed and is propagating through the layer (see Fig. 4.1). . . . . . . . . . . . . . . . . Mach number (M a) (dashed line for continuum and triangles for MD) and estimated mean free path (/) (solid line for continuum and circles for MD) as functions of z/ when f t = 0.22. At this time the shock is fully developed and is propagating through the layer (see Fig. 4.1). . . . . . . . . . . . . . . . . The maximum Mach number (M a) (dashed line for continuum and triangles for MD) and the ratio of shock width d to mean free path (solid line for continuum and circles for MD) as functions of time f t. The shock begins forming at f t = 0.12 and leaves the layer at f t = 0.35. . . . . . . . . . . . . . . . . Average dimensionless shock speed vshock / g in the reference frame of the plate as a function of coefficient of restitution e. vshock is calculated as the average speed of the shock from when the shock is formed until it leaves the layer. . . . . . . . . . . Location of the shock ( for constant e = 0.9, for velocitydependent e = e (vn )) and the center of mass of the layer ( for constant e = 0.9, + for velocity-dependent restitution e = e (vn )) as a function of time f t during one cycle of the plate (thick solid line) in the MD simulation. . . . . . . . . . . . . (a) The height of the center of mass of the layer, and (b) the location of the shock as a function of time f t for various layer depths H. The trajectory of the center of mass (a) is shown for a complete cycle in the oscillatory state from the MD simulation with H = 4.5 ( ), H = 9 ( ), and H = 13.5 ( ). The center of mass location from continuum simulations agrees with MD to within a root mean square difference of 3% over one cycle in each of the three cases. The shock location (b) is shown for the fraction of the cycle in which the shock is in the layer, from the oscillatory state of H = 4.5 ( for MD, - - for continuum simulation), H = 9 ( for MD, -- for continuum), and H = 13.5 ( for MD, - - - for continuum ). . . . . . An overhead view of a layer of grains, showing the center of mass height zcm as a function of horizontal position (x, y) in a cell with horizontal dimensions Lx Ly = 126 126, from (a) MD simulations, and (b) continuum simulations. Peaks of the layer corresponding to large center of mass height zcm are shown in white; valleys corresponding to low zcm are shown in black. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii 60 61 62 66 68 70 5.1 79 5.2 5.3 5.4 Dispersion relations for stripes which form perpendicular to the long dimension of cells with horizontal dimensions 168 10. Data for continuum simulations are shown as triangles, MD simulations as circles, and points where continuum and MD simulations yield the same wavelength are shown as squares. In both continuum and MD simulations, the dominant wavelength of the final oscillatory state fits very well to the dispersion relation found in experiments = 1.0 + 1.1f -1.32 0.03 (solid line) [113]. Error bars in both cases are calculated exclusively from discretization due to periodic boundary conditions in a finite size box. . . . . . . . . . . . . . . . . . . . . . . . . . . Views of MD simulations from the side of the cell, looking down the length of stripes oriented parallel to the y-axis at various times f t during two cycles of the plate. All particle locations within the range (0 y 10) are projected into the x-z plane, regardless of their y-coordinate, and particle positions are represented by open blue circles. Particles that appear to overlap are at different y-coordinates (in front of or behind each other). The plate height is represented as a thick solid red line. The full cell is 42 42 in horizontal extent, so that the entire span of the cell is shown in the x-direction. The cell extends to a height of z = 200, but the figures show only z 20 since no particles fly higher than z = 20 during these two cycles. . . Volume fraction from continuum simulations at a slice y = 21 parallel to the x - z plane at various times f t during two cycles of the plate. Empty space ( = 0) is white, random close-packed volume fraction ( = 0.65) is black, and the color increases from white to black through shades of blue as volume fraction increases. The plate height is represented as a thick solid red line. The full cell is 42 42 in horizontal extent; the entire span of the cell is shown in the x-direction. The cell extends to a height of 80 above the lower plate, but the figures show only z 21, since the density is quite low above this height. No particles fly higher than z = 20 during this time in the MD simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 84 88 xviii 5.5 5.6 5.7 5.8 Dimensionless horizontal velocity in the x- direction ux / g from continuum simulations, through a slice y = 21 parallel to the x-z plane for various times f t during two cycles of the plate. Flow to the right ux > 0 is represented by red colors, flow to the left (ux < 0) is represented by blues, and the darkness of the color increases from white (ux = 0) to black as the magnitude of ux increases. The plate height is represented as a thick solid black line. The full cell is 42 42 in horizontal extent; the entire span of the cell is shown in the x-direction. The cell extends to a height of 80 above the lower plate, but the figures show only z 21. . . . . . . . . . . . . . . . . . . . . . . . Dimensionless Vertical velocity uz / g from continuum simulations, through a slice y = 21 parallel to the x - z plane for various times f t during two cycles of the plate. Upwards flow uz > 0 is represented by red colors, downwards flow (uz < 0) is represented by blues, and the darkness of the color increases from white (uz = 0) to black as the magnitude of uz increases. Velocity is shown in a stationary reference frame in which the plate is oscillating sinusoidally between z = 0 and z = 2A. The plate height is represented as a thick solid black line. The full cell is 42 42 in horizontal extent; the entire span of the cell is shown in the x-direction. The cell extends to a height of 80 above the lower plate, but the figures show only z 21. . . Dimensionless granular temperature T / (g) from continuum simulations, through a slice y = 21 parallel to the x - z plane for various times f t during two cycles of the plate. Large T is white, T = 0 is black, and the color increases from black to white through shades of red as horizontal velocity increases. The plate height is represented as a thick solid green line. The full cell is 42 42 in horizontal extent; the entire span of the cell is shown in the x-direction. The cell extends to a height of 80 above the lower plate, but the figures show only z 21, since the density is quite low above this height. No particles fly higher than z = 20 during this time in the MD simulations. Volume fraction (left column) and granular temperature T (right column) from continuum simulations at a slice y = 21 parallel to the x - z plane at times 0.875 f t 1.25 when the shock is forming and propagating through the layer. The colorbar below each column represents the scale for all pictures in that column. The plate is a thick red line in the left column and a thick green line in the right column. . . . . . . . . . . . 89 90 95 98 xix 5.9 Vertical velocity uz (left column) and granular pressure P (right column) from continuum simulations at a slice y = 21 parallel to the x - z plane at times 0.875 f t 1.25 when the shock is forming and propagating through the layer. The colorbar below each column represents the scale for all pictures in that column. The plate is a thick black line in the left column and a thick green line in the right column. . . . . . . . . . . . . . . . . . 5.10 Horizontal velocity ux (left column) and the horizontal derivative of granular pressure x P = P (right column) from conx tinuum simulations at a slice y = 21 parallel to the x - z plane at times 0.875 f t 1.25 when the shock is forming and propagating through the layer. The colorbar below each column represents the scale for all pictures in that column. The plate is a thick black line in both columns. . . . . . . . . . . An overhead view of the layer of grains, showing the center of mass height zcm (x, y) of the layer as a function of location in the box, for (a) MD simulations with a plate acceleration with respect to gravity = 1.9, (b) MD simulations with = 2.2, (c) continuum simulations with = 1.9, and (d) continuum simulations with = 2.2. Peaks corresponding to large zcm are shown in white, while valleys corresponding to small zcm are shown in black. The grayscale for all four images use the same units as shown by the colorbars on the right of the figure. . . The mean square deviation 2 of the local center of mass height from the average center of mass height of the entire layer as a function of changing accelerational amplitude for MD (triangles) and continuum (circles) simulations. The vertical dotted line represents the onset of stripe patterns in the continuum simulations (c = 1.955). For this figure, 2 was calculated using bins of size xbin = ybin = 2 for both continuum and MD simulations. . . . . . . . . . . . . . . . . . . . . . . 99 100 6.1 107 6.2 111 xx 6.3 Comparison of MD simulations (triangles) to the Swift Hohenberg model (solid lines) for (a) 2 , and (b) Pmax as a function of control parameter (bottom axis) for SH, and (top axis) for MD. This plot is obtained by using the parameters for SH simulations which give the best fit between MD and SH simulations. This fit yields a noise strength F = (1.2 0.2) 10-2 and a delayed onset of long range order M F = 0.094 correspondc ing to LR = 2.15 0.01 in MD (the vertical dotted line in c the figure). The global ordering jumps sharply at M F = 0.10 c (corresponding to LR = 2.15 in MD), representing a transic tion to stripe patterns, while 2 increases smoothly through that transition. This fit predicts a mean field onset value of M F = 1.965 0.007, corresponding to M F = 0 (the vertical c c dashed line in the figure.) . . . . . . . . . . . . . . . . . . . . The squared residual R between and as a function of the noise strength F used in SH simulations. The best least squares fit is given by F = (1.2 0.2) 10-2 . . . . . . . . . . 2 2 2 M D 2 SH 112 6.4 117 6.5 Comparison of growth of normalized by the mean center of mass height of the layer 2 2 / zcm 2 = (zcm - zcm )2 / zcm 2 for MD simulations with f = 0.25 (squares) and f = 0.4174 (triangles). Both frequencies use a layer depth H = 5.4, with a cell of size 2 2 in the horizontal direction, where the dominant wavelength = 21 for f = 0.4174 and = 42 for f = 0.25. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 xxi Chapter 1 Introduction Systems composed of many interacting macroscopic particles are collectively known as "granular materials" or "granular media." Sand in a dune, rice at the supermarket, pills at a pharmaceutical factory, and rocks in an avalanche are all examples of granular materials. Granular systems can demonstrate huge ranges in size and behavior: from an hourglass to planetary rings; from a static pile of pebbles to an avalanche down a mountainside. Despite the ubiquity of grains in nature and industry, there are many basic questions about the collective behavior of granular materials which are still unanswered. The interactions between individual grains may be well characterized based on the properties of the particles (friction, coefficient of restitution), the type of interaction (static contact or dynamic collision), and whether the particles are in a vacuum (dry) or interact with an interstitial fluid (wet). However, general governing equations for the collective behavior of thousands of grains are not known, and at least until recently the applicability of continuum equations for granular fluids was under debate [25, 36, 60]. An important characteristic of granular materials is that collisions between grains are inelastic. Thus a system of grains which is initially in motion 1 will gradually dissipate energy through collisions until the particles eventually come to rest. Grains may be driven away from a static equilibrium by an input of energy, frequently through forcing from a boundary. If sufficiently driven, granular materials may be excited into a regime in which there are no long-lasting contacts between particles, but in which particles interact with each other mainly through inelastic collisions of short duration compared to the time scales of the macroscopic flow. This regime is known as the regime of rapid granular flow [17]. These regimes may be considered in analogy to states of matter, in which granular particles behave as a solid, liquid, or gas. Thus a static pile in which there are long-lasting interactions between particles with a static configuration displays solid-like properties, while we call granular materials in the rapid flow regime "granular fluids." Granular materials have been studied extensively in static systems as well as in the rapid flow regime (see reviews and books in [4, 17, 26, 50, 52, 85, 109]). In addition, studies have investigated the transition as a granular fluid cools through inelastic collisions [38], and as static piles are fluidized through driving [81]. In this study, we focus on grains in the fluidized regime. In analogy to molecular fluids, we characterize granular fluid flows by the number density n of particles, the net velocity u due to collective motion of particles, and the "granular temperature" T , where 3 T 2 is the average kinetic energy due to random particle motion. Finally, granular materials may be affected by interactions between the particles and an interstitial fluid. Granular flows have been studied in the 2 presence of air between particles [15, 84, 122], in beds of grains fluidized by liquid or gas flow [28, 35], and in evacuated containers [77, 113]. The effect of interstitial gas or liquid is significant in many cases; however, for particles of diameter larger than 1 mm, air friction is often negligible if the velocity is not too large [93]. We focus in this work on modeling dry granular materials (grains in a vacuum) in the rapid flow regime. We use a hydrodynamic-type model to test whether such an approach can accurately model dry grains in the rapid flow regime. 1.1 1.1.1 Modeling of Granular Flows Molecular Dynamics Simulation Molecular dynamics (MD) simulations are commonly used to simulate granular materials [26, 47]. Molecular dynamics simulations are a type of numerical simulation which models a collection of particles (or molecules) by tracking the trajectories of each particle between collisions, and then calculating the outcome of collisions between particles. MD simulations have been carried out for granular materials as well as for simulations of fluids and solids (see reviews in Refs. [2, 19, 90].) MD simulations directly model the grains at the particle level, integrating Newton's equations for particles between collisions, and simulating particle-particle and particle-boundary collisions. Thus these simulations may be carried out for systems in which the governing equations for collective be3 havior are unknown, without additional approximations aside from modeling the behavior of individual particles. These simulations become computationally more time consuming as the number of particles increases, however, as each particle in the system must be tracked. Grains may be modeled as either hard spheres or soft spheres, depending on the type of simulation. For a review of these types of simulations, please see Refs. [26, 47]. Soft sphere models allow particles to penetrate each other, with an interaction force which depends on penetration depth. Such a simulation is particularly good for solid-like regimes, because the soft-sphere model allows finite interaction time between particle. It also lends itself well to the addition of interactions between particles and an interstitial fluid or other long-range forces. Hard sphere models approximate particles as hard spheres with finite diameter , and collisions between particles are modeled as binary, instantaneous collisions when the particle centers are separated by exactly (or the average of the diameters of the two colliding particles in polydisperse systems). For event driven hard sphere simulations, the position and velocity of each particle is specified at an initial time t0 . From its initial position and velocity, the trajectory of each particle is known to be simply a parabola as it falls under gravity until it collides with another particle or a wall. At the time of the first collision, the program calculates the new positions and velocities of the colliding particles, and uses their post-collisional velocities as their new initial conditions. In this way, event driven simulations only need need to make new 4 calculations at times when collisions occur. This method tends to be much faster than the soft sphere method, in which interactions are present at all times. However, the hard sphere model allows no penetration or deformation of particles, and cannot simulate solid regimes in which the finite duration of collisions is important. In either hard- or soft-sphere models, interactions between particles and between particles and walls may include inelasticity, frictional effects, or other (e.g. electrostatic) effects. In this work, we use event-driven hardsphere MD simulations for monodisperse, inelastic, frictionless particles in which electrostatic and other effects are negligible. This model is described in more detail in Chapter 3. 1.1.2 Continuum Equations In molecular fluids, the number of particles involved in a system frequently makes direct molecular dynamics simulation impractical. The most common approach to ordinary fluids is to view the fluid as a continuum, and to use continuum models to describe the time evolution of continuum fields such as density , mean velocity u, and temperature T . These variables are generally described for ordinary fluids by the well-known Navier-Stokes (NS) equations. In addition to allowing for the investigation of systems of many particles, the approach of using continuum models has allowed for the development of sophisticated theoretical tools for understanding properties of fluid flow. 5 Both analytical and numerical techniques such as linear stability theory, amplitude equations, and large numbers of computational fluid dynamics (CFD) techniques are predicated upon a continuum model for fluid dynamics. Unlike ordinary fluids, in which the NS equations are well-established for most fluid flows, a similar set of equations is not established for granular materials. In fact, the situation for granular materials is in many ways opposite to that of molecular fluids: in ordinary fluids, the Navier-Stokes equations were originally proposed as conservation laws based on fluid as a continuum [65]; it was only later that Boltzmann proposed kinetic theory for a rarefied gas composed of individual modecules [16]. In granular theory the "microscopic" approach of individual particles is well understood, while a set of established continuum equations is missing. There are several barriers to establishing a continuum theory of granular materials [110,125]. The first is a challenge to the basic continuum hypothesis. Continuum theory requires small variations in space over long distances so that averages over small fluid elements can replace the behavior of collections of individual particles. This generally requires that microscopic length scales (mean free path of a particle) be much smaller than the smallest macroscopic length scale of the flow field. However, in this dissertation alone, patterns are studied with wavelengths in the tens of particle diameters (Chapter 5 and Chapter 6), in layers 5 to 10 particle diameters deep (Chapters 4, 5, and 6). In these cases, the mean free path of a particle is on the order of a particle diameter or less (Chapter 4). For granular materials, significant separation 6 of scales may only be present for very low density and nearly elastic particles [37, 110, 125]. In addition, the inelasticity of particles means that any granular system in a rapid flow regime is far from equilibrium. Assumptions based on thermodynamics and equilibrium statistical mechanics which hold for molecular fluids near equilibrium may break down for granular systems [17]. Derivations of continuum equations from kinetic theory generally assume a velocity distribution function which is nearly Gaussian in form. However, granular experiments and MD simulations have displayed velocity distributions that are dependent on experimental geometry and details of forcing [64,83,99,116,117]. Finally, for many granular systems, it is not known what effects are negligible and which ones need to be included in the continuum model. For example, the continuum equations used in this paper [59] neglect frictional effects. However, since this study began, some experiments [39] and simulations [80] have shown that friction plays a significant role in pattern formation. Despite these difficulties, many continuum models have been proposed for granular materials. A common approach is to use a kinetic theory approach similar to that used to derive continuum equations for elastic gases in order to derive continuum balance equations for mass, momentum, and energy from a microscopic model of particle interactions. The variety of approximations used in such derivations is remarkable. Such equations have been derived to Euler, Navier-Stokes, and Burnet orders for fluids of inelastic particles with and without friction; in the limit of dilute or nearly elastic gases; and with 7 variations such as modifications to Fourier's heat law [1,12,42,45,56,57,59,71, 73,89,103]. Additional study has been dedicated to developing an appropriate set of boundary conditions for solving these continuum equations for granular flows [13, 53, 58, 105, 106]. With so many possibilities, there is a small amount of data comparing various models directly to experiment or to MD simulations for experimentally realizable systems [45, 54, 93]. In addition, attempts to use such continuum equations to predict properties of granular media have mainly been restricted to idealized systems that are difficult to produce in experiment [43, 45, 54, 61] and to steady state or asymptotic time limits [43, 45, 54, 93]. The viability of these continuum equations for experimentally realizable, three-dimensional (3D), time-dependent granular systems has scarcely been examined. In this work, we compare three-dimensional simulations of continuum equations derived to Navier-Stokes order for monodisperse, frictionless, nearly elastic particles to MD simulations which have been previously validated by comparison to experimentally observed standing wave patterns for varying control parameters with frictional particles [7,79]. We test the ability of such a continuum model to reproduce behavior seen in experiment and MD simulations, and identify potential improvements to such a model by comparison to MD simulation. 1.2 Vertically Oscillated Granular Layers Dry, rapid granular flows are usually driven by forcing from the bound- aries, frequently by vibration, either vertically or horizontally [14, 24, 32, 63, 8 94, 111]. Vertically oscillated granular layers represent an important testbed for granular materials, and have been the subject of extensive investigation in recent years [7, 18, 29, 33, 101, 112, 114]. The system consists of a layer of grains on an impenetrable plate that oscillates sinusoidally in the direction of gravity with amplitude A and frequency f . The dimensionless peak acceleration is = 4 2 f 2 A/g, where g is the acceleration due to gravity. A sample experimental setup for studying shaken layers is shown in Fig. 1.1. This simple oscillating system has been found to yield standing wave pattern formation [75], convection [62], clustering [31], steady-state flow fields far from the plate [11], and shocks [44]. The shaken layer thus provides a useful system for testing continuum models of granular layers. We model vertically shaken layers in our simulations to investigate shocks and pattern formation, and compare results from MD and continuum simulations. 1.3 Patterns in Vertically Oscillated Granular Layers A layer of grains on a plate oscillating sinusoidally in the direction of gravity with frequency f and amplitude A will leave the plate at some time in the cycle if the maximum acceleration of the plate is greater than that of gravity ( > 1). The layer dilates above the plate, then collides with the plate later in the cycle, compressing the layer on the plate. Above a critical value of acceleration, standing wave patterns spontaneously form in the layer. These standing waves produce horizontal variation in the layer with peaks 9 (b) Figure 1.1: Electromagnetic shaker system (a) for visualizing patterns in oscillated granular layers. The container is illuminated by light incident at low angles from the side (arrows in (b)). The light reflects off the particles and is visualized by a CCD camera, yielding bright spots for peaks and dark spots for valleys of the pattern. From [113]. 10 (a) (b) (c) Figure 1.2: Stripe patterns (a), square patterns (b), and hexagonal patterns (c) formed in vertically oscillated granular layers. All patterns oscillate at f /2, where f is the frequency of oscillation of the plate; in the hexagon pattern, the right and left sides of the cell are out of phase with respect to each other (peaks in the center of the hexagon on the right, peaks at the border on the right). The shaking parameters , f , and H are given by (a) 3.3, 67 Hz, 7, (b)2.9, 25 Hz, 4, and (c) 4.0, 67 Hz, 7, respectively. In (a) and (c), the diameter of the container is 770 and in (b), the container is square with sides 1100, where the particles used are bronze spheres of diameter = 0.165 mm. (a) and (c) are from [76], (b) is from [39], and the entire figure was composed in [109]. corresponding to deep areas of the layer and valleys corresponding to shallow areas of the layer (Fig. 1.2). These standing wave patterns are subharmonic with respect to the plate, with a frequency of f /2, where f is the frequency of the oscillation of the plate [75]. Various subharmonic standing wave patterns, including stripe, square, and hexagonal patterns, have been found experimentally, depending on the nondimensional frequency f = f H/g and the nondimensional accelera- tional amplitude [75]. As is increased just above the onset of patterns, the layer forms stripe or square patterns, depending on the value of nondi- 11 mensional frequency f . For f below a critical value, squares are formed, while for f above this value, stripes form in the container (Fig. 1.3). The value of the stripe-square transition varies with different , but occurs near f 1/3 [113]. The behavior near onset is repeatable for a wide range of layer depths and of particle sizes, including for polydisperse systems with particle size ranges of up to approximately 30% and irregular particles such as rice grains [114]. As is increased further, additional bifurcations to different pattern regimes are found, including hexagonal patterns, f /2 flat layers and f /4 patterns in which the layer collides with the plate only every other cycle, and disordered "labyrinthian" patterns (Fig. 1.3). 1.4 The Role of Friction in Shaker Patterns Frictional properties of grains play an important role in pattern forma- tion in vibrated granular layers. Examination of the center of mass of a peak in a square pattern shows harmonic motion about an equilibrium position of the peak, with a natural frequency known as the lattice frequency fL . Experiments in which the plate frequency f was modulated by the lattice frequency fL , show that for sufficient modulation amplitude, defects form which destroy the long range order of square patterns [39]. Adding graphite to the system acts as a lubricant which reduces the friction between particles; the addition of graphite increases the disorder introduced by the modulation frequency, reducing the system to a state of disordered peaks (Fig. 1.4). 12 9 8 7 unexplored randomly moving labyrinths (disordered) phase bubbles stripes hexagons (f/4) squares (f/4) stripes (f/4) 6 5 4 3 2 0.2 0.4 squares (f/2) flat (f/2) hexagons (f/2) stripes (f/2) flat 0.6 0.8 1 f* Figure 1.3: Phase diagram for a vertically oscillated granular layer of bronze spheres of diameter = 0.165 mm, where the layer depth H = 5, and the container diameter is 770. The diagram indicates the granular patterns observed as a function of the dimensionless accelerational amplitude and dimensionless frequency of the plate oscillation f . The transition from a flat layer to squares is hysteretic; solid lines indicate the transition for increasing , dotted lines indicate the transition for decreasing . Similar patterns are observed for a variety of particle types and layer depths. From [79]. 13 Figure 1.4: Defect creation and "melting" of a square pattern in a shaken granular layer for (a) experiment and (b) MD simulation for = 0.17 bronze spheres oscillated at = 2.9 and f = 32 Hz. In the experiment, graphite has been added to reduce friction, and at t = 0, a frequency modulation is added to experiment at fL = 2 Hz. In MD simulation, no frequency modulation is added, but at time t = 0, friction was turned off in the simulations. The numbers represent the number of oscillations since the chagne at t = 0. The insets in each picture show the power spectrum (white is low power, black is high power) corresponding to each image. A similar effect may be observed in MD simulations where the friction between particles is reduced. If the friction is reduced to zero, the square patterns "melt" and become disordered peaks even without frequency modulation (Fig. 1.4). A square lattice is identifiable by four equally spaced peaks in the power spectrum; a disordered pattern of a given wavelength is represented by a ring in fourier space. In Chapters 5 and 6, we investigate f /2 patterns near the onset of 14 standing wave patterns, using MD and continuum simulations. We focus on the region near onset where f /2 stripe and square patterns are observed in shaken layers of frictional particles, although our continuum simulations model frictionless particles. Most of our investigation of patterns is in the regime in which stripe patterns are observed in both frictional and frictionless systems, although we also examine the onset of disordered peaks for frequencies in which squares form for frictional particles. 1.5 Shocks in Granular Materials The effective speed of sound (the speed of propagation of small am- plitude pressure fluctuations) is typically on the order of several cm/s for granular particles [109]. This is four orders of magnitude smaller than the 330 m/s speed of sound in air. Thus a distinguishing feature of granular materials is that granular flows reach supersonic speeds under common laboratory conditions [41, 45, 93]. Hence shock formation, which is achieved only under extreme conditions in ordinary gases, is commonplace in granular materials. Shock propagation has been extensively studied in rarefied gases [21, 124], but the similarities and differences between shocks in ordinary fluids and in granular media is the subject of ongoing research in theory [5, 43], simulation [10, 86, 93] and experiment [44, 93, 95, 96]. A recent laboratory and molecular dynamics study [93] examined time-independent behavior of oblique shocks in granular flow between two closely spaced plates, and the results were compared with 2D simulations of Navier-Stokes order continuum equations 15 Figure 1.5: Horizontal component of the velocity field (vx ) of a flow of stainless steel spheres with diameter = 1.2mm falling under gravity onto a wedge as determined by (a) experiment, (b) MD simulations, and (c) continuum equations to Navier-Stokes order. Each picture shows a region of size 130 104. A shock separates an undisturbed region in which the horizontal component of flow is nearly zero, from a compressed region in which streamlines (solid lines with arrows) are diverted around the obstacle. From [93]. (Fig. 1.5). Experiments [44] and MD simulations [3, 86] suggest that shock waves form in shaken granular layers confined to two-dimensions as the layer contacts the plate. In Chapter 4, we study shock formation and propagation in simulations of three-dimensional shaken granular layers. 1.6 Outline This dissertation aims to investigate the potential for hydrodynamic- type continuum equations to accurately model rapid granular flows. We discuss previous progress towards solving continuum models of granular materials in Chapter 2. To explore the potential of granular hydrodynamics, we seek to estab- 16 lish a correspondence between a continuum model of granular materials and a MD simulation which has been previously compared to experiments. We compare the results from these simulations in modeling a common experiment in which frictional MD simulations have previously been verified: that of the vertically oscillated granular layer. The simulation models and the numerical methods used to integrate these models are described in Chapter 3. In addition to exploring the strengths and weaknesses of continuum theory, this work aims to use these simulations to investigate granular phenomena in vertically shaken layers. The two main phenomena we investigate are shocks and pattern formation in shaken layers. In Chapter 4, we use our simulations to study shock propagation in vertically shaken granular materials. The results from our simulations are used to investigate the role of inelasticity in shock formation and propagation, and to study the dynamics present upon collision of a granular layer with a solid plate. We investigate the formation of patterns in continuum simulations, and compare these patterns to those found in MD simulations in Chapter 5. In addition, we use the results from continuum simulations to explore the dynamics of these patterns and the link between shock propagation and standing wave patterns in shaken layers. Aspects of standing wave patterns are also investigated in Chapter 6, where we use MD simulations to probe the effect of finite particle number on 17 pattern formation. Results from MD simulations are compared to the SwiftHohenberg model equation with the presence of noise. Continuum simulations are used here as a control to compare noisy patterns to those found in a model which views the system as a true continuum. Chapter 7 returns to the question of the validity of the continuum theory, and proposes procedures to modify the continuum equations based on results from the previous sections. This chapter discusses the importance of frictional effects and effects of fluctuations which are not present in the current continuum model. Finally, Chapter 8 summarizes the major results from this dissertation, as well as indicating possible directions for future work. 18 Chapter 2 Previous Solutions to Granular Continuum Models Despite the variety of proposed continuum equations (cf Section 1.1.2), to this date there are few instances in which granular hydrodynamic equations have been solved for experimentally realizable systems. This dissertation represents an attempt to expand this body of work and to test the applicability of hydrodynamic equations for modeling rapid granular flows. Simplified continuum equations, such as amplitude equation models, have been shown to reproduce important aspects of specific granular systems, including pattern formation in oscillated granular layers. For example, a phenomenological model characterizing the layer thickness based on local conservation of mass produces waves similar to the square patterns and solitary wave oscillons observed experimentally [18]. A similar amplitude equation model with a different forcing mechanism yields stripes, squares, and oscillons [112]. Recent experiments [40] have shown that noisy fluctuations in the Swift-Hohenberg equation resemble fluctuations near the onset of patterns in granular materials; in Chapter 6 we use this same model equation to investigate noise in MD simulations. Finally, a continuum shallow water-type theory 19 has produced patterns with a dispersion relation matching those from experiment [27]; a recent experiment demonstrated that dragging a cylinder through a vibrofluidized granular bed produces a Mach cone which matches the Mach relation angle predicted from shallow water theory [46]. These simplified models, as well as other model equations [98, 100, 104, 118], provide important evidence that continuum models can capture some significant granular phenomena. However, these models all have the disadvantage of drastically simplifying the system by considering only certain specific aspects of the flow (e.g. layer thickness or horizontal motion) and neglecting other important phenomena, (e.g. the dynamics of collision of the layer with the plate). Thus, these models cannot provide a general set of governing equations for rapid granular flow, but are only able to act as models for specific granular experiments. General hydrodynamic equations, such as those used in this dissertation, have been solved for a few simplified cases. Attempts to find analytical solutions to granular continuum equations have generally involved simplifications difficult to reproduce experimentally. For instance, cases of shear flow [54] have been analyzed assuming constant granular temperature (difficult to produce in an inelastic gas forced through the boundaries). An early attempt to describe granular flow using hydrodynamic equations [45] found analytic solutions for several cases, including a "steady-state system with no flow," "steady-state Couette flow with no gravity," and "steady-state Couette flow in a gravitational field." While this paper represented a step forward in 20 understanding the role of inelasticity in hydrodynamic equations, the limitations of this work were best expressed by the author: "This paper has been restricted principally to a consideration of analytic solutions to certain static and steady-state problems. This has usually meant, especially in the case of non-zero flow fields, that it was possible to solve only for elastic grains in a gravitational field, or for inelastic grains without gravity. Since all laboratory experiments must involve inelastic grains and since most experiments must also be performed in a gravitational field, it is not generally possible to compare the present calculations with the result of flow experiments, even if the uncertainties over the proper boundary conditions are resolved." Shocks created by a piston moving into a semi-infinite uniform medium have been studied by Euler-like hydrodynamic equations in one-dimension [41]. These equations are shown to yield an analytic result for the shock solution in the asymptotic time limit. In addition, numerical integration of the time evolution of these equations exhibit solutions which evolve to match this asymptotic solution. Despite the fact that this analysis requires reducing the system to a one-dimensional problem with a semi-infinite uniform medium, recent MD simulations of a shock forming between two wedges seem to agree with predictions of the shock speed from the asymptotic solution [92]. This problem holds much interest and merits further investigation. Recent studies of granular flow past a wedge have begun to establish a connection between experiments and continuum simulations [93]. In these experiments, molecular dynamics simulations were able to quantitatively re21 produce shock properties observed experimentally. Continuum simulations were carried out in two-dimensions to find the steady- state behavior of the shock. These simulations exhibited shocks which were qualitatively similar to the shocks found in experiment (cf Fig. 1.5). The continuum simulations quantitatively reproduced results from MD simulations using frictionless particles, but differed quantitatively from experiment and MD simulations with friction To my knowledge, the work in this dissertation represents the first simulations of fully three-dimensional, time-dependent flow using a complete set of granular hydrodynamics equations. These simulations are used to model rapid granular flow in a granular shaker, a system which has been extensively studied experimentally. For the majority of the paper, these simulations are compared to MD simulations using frictionless particles; the addition of friction is the only change necessary to reproduce experimental results using these MD simulations [7] (see Chapter 3 for a more thorough discussion of the MD simulations). In addition, wavelengths observed in continuum simulations of patterns in shaken granular materials yield good agreement when compared directly to a dispersion relation from experiment (see Chapter 5). 22 Chapter 3 Simulation Methods To evaluate the potential for a continuum hydrodynamics-type approach to model granular media, we have developed a C language based code to simulate shaken granular layers based on hydrodynamic equations to NavierStokes order. We use these simulations to investigate properties of shocks and pattern formation in shaken granular layers. We compare these simulations to results from a molecular dynamics code which has previously reproduced patterns found experimentally. MD simulations for this work were executed by Sung Joon Moon and Jennifer Kreft. The continuum code and the molecular dynamics code are the main tools we use in this work to study granular media. In Chapter 6, we interpret results from MD and continuum using simulations of the Swift-Hohenberg model equation to investigate the role of noise on the onset of patterns in this system (we will discuss the Swift-Hohenberg simulations in Chapter 6.) 3.1 A Continuum Model of Granular Materials We numerically solve the continuum equations to Navier-Stokes order proposed by Jenkins and Richman [59]. These equations were derived for a 23 dense gas composed of frictionless (smooth), inelastic hard spheres of uniform diameter by applying Grad's 13-moment method to the inelastic EnskogBoltzmann equation. These equations are derived to full Navier-Stokes order, and include formulas for calculating the transport coefficients. These equations do not include friction between particles, however. The effect of inter-particle friction will be discussed in Chapter 7. The Jenkins-Richman (JR) model yields hydrodynamic equations for number density n (or equivalently, volume fraction = n 3 ), 6 velocity u, and granular temperature T. Here 3 T is the average kinetic energy per unit 2 mass due to random particle motion: T = 1 3 |v - v |2 , where v is the velocity of a particle, and the brackets represent an average over all of the particles in a small fluid element at an instant of time. Throughout this dissertation, we use units such that the mass of a particle m = 1. Previous MD simulations have shown that the mean square deviation may vary for different components of velocity, and that the value of the vertical component of temperature: Tz = |vz - vz |2 is more than double the value of the 1 2 horizontal components of temperature: Tx = |vx - vx |2 + |vy - vy |2 at some locations in a vertically oscillated granular layer. In this work, however, we are able to acheive good comparisons to MD simulations (see Chapter 4) using only a scalar model of granular temperature: T = 1 (2Tx + Tz ). 3 Ordinary hydrodynamic equations express conservation laws for mass, momentum, and energy. For granular media, mass and momentum are conserved, but energy is dissipated through inelastic collisions, as characterized 24 by the normal coefficient of restitution e. The Jenkins-Richman (JR) model express conservation of mass and momentum, and the energy equation is modified by the addition of a temperature loss term : n + (nu) = 0, t (3.1) n u + u u t z = P - ng^, (3.2) 3 n 2 T + u T t = - q + P : E - , (3.3) where the components of the symmetrized velocity gradient tensor E are given by: Eij = 1 (j ui + i uj ) . The components of the stress tensor P are given by 2 the constitutive relation: 2 Pij = -p + ( - )Ekk ij + 2 Eij , 3 and the heat flux is calculated from Fourier's law: q = -T. (3.5) (3.4) In the granular shaker system under consideration, the volume fraction of particles reaches the same order of magnitude as the random close-packed volume fraction in some regions of the cell (cf Chapter 4), leading to correlations between particle positions. Therefore, we use a model that does not approximate the particles as uncorrelated point particles, but rather models 25 the system as a dense gas, where finite particle size and resulting correlations between particle positions are accounted for by the radial distribution function at contact. To calculate the pressure, we use the equation of state and radial distribution function at contact proposed by Goldshtein et al. [43] to include both dense gas and inelastic effects: p = nT [1 + 2(1 + e)G()] , where (3.6) G() = g0 (), and the radial distribution function at contact 4 3 max (3.7) g0 () = 1 - max -1 , (3.8) where max = 0.65 is the 3D random close-packed volume fraction. The form of equations (3.1-3.3) are identical to those for a compressible gas of elastic particles, aside from the energy loss term in Eq. 3.3: 12 nT 3/2 = (1 - e2 ) G(). (3.9) Finally, the JR equations are completed by formulas for the values of 26 the transport coefficients. The bulk viscosity is given by 8 = nT 1/2 G(), 3 the shear viscosity by = 4 5 1 nT 1/2 +1+ 6 16 G() 5 1+ 12 G() , (3.11) (3.10) and the thermal conductivity by 5 1 6 15 nT 1/2 +1+ = 16 24 G() 5 1+ 32 9 G() . (3.12) For an oscillating granular layer, incompressible or isothermal approximations are unrealistic. In shaken layers studied here, n and T vary by orders of magnitude in both space and time during each cycle of oscillation (cf Chapter 4.) Thus the viscosities and thermal conductivity are very far from constant and must be calculated at each timestep at each location in the cell. 3.2 A Continuum Simulation of Granular Shakers To simulate shocks and pattern formation in shaken granular layers, we have developed a C language code to solve the JR continuum equations. This code has been specifically developed to solve the JR equations in a threedimensional rectangular box. It could be modified to simulate general geometries and other granular systems by changing boundary conditions and the grid on which it is solved. 27 3.2.1 Continuum Fields The code is designed to integrate the continuum equations (3.1-3.3) forward in time, starting from specified initial conditions. This integration yields number density n, velocity u, and granular temperature T as functions of space and time. To integrate these equations, we use finite differencing on a fixed grid in 3D. The fields n, u, and T are treated as 3D arrays of size Nx Ny Nz in the two horizontal directions x and y, and the vertical direction z, respectively. In addition, the pressure, radial distribution function at contact, and the transport coefficients depend on the value of local fields n and T ; so P , , , G, , and are also treated as fields in three-dimensions, and specified by three dimensional arrays: double ***n,***T,***u,***v,***w; double ***P,***mu,***K,***G,***lam,***gam; The equations are solved on a fixed grid in 3D, so the number of gridpoints in each direction Nx , Ny , and Nz are defined by the user, but remain fixed after that point. 3.2.2 User input The user provides the main parameters for the system to be modeled. For each simulation, the shaking parameters , f , the coefficient of restitution e, and the layer depth H are provided by the user, making it possible to sim- 28 ulate a variety of particle types, depths, and shaking parameters. In addition, the user may vary the number of gridpoints in each direction Nx , Ny , and Nz , and the corresponding physical linear dimensions of the box Lx , Ly , and Lz respectively. The size of the gridpoints in each direction x = Lx /(Nx - 1), y = Ly /(Ny - 1), and z = Lz /(Nz - 1) are chosen for stability and for the resolution required to visualize shocks and other phenomena, as described later in this chapter. 3.2.3 Initial conditions ^ ^ ^ The five main fields n, u = ux x +uy y +uz z , and T are initialized for the first timestep in the simulation. The user may choose to start from an artificial flat layer or to read in the initial conditions from a file. For the artificial flat layer, densities are assigned to create a flat layer above the plate with random fluctuations of user selectable strength added to the density at each gridpoint to create horizontal variation (we use approximately 10 % variation for initial flat states throughout this work.) The horizontal components of velocity are set identically to zero, and the layer is given a uniform temperature and a small, uniform downward velocity (toward the plate). In addition to an artificial flat state, the initial conditions may also be read from a file. For example, if the user wishes to restart a simulation from the end of a previous run, he or she may read in a file saved from the end of that previous simulation. Once the fields are initialized, the user may specify any of the simulation parameters, including the shaker parameters and the 29 coefficient of restitution. This procedure can be used to set up a steady state and then change the acceleration or frequency of the plate, or to investigate the effect of inelasticity on a system starting from identical initial conditions, as in the shock study in Chapter 4. Once the main fields have been initialized, the initial values of pressure, radial distribution function at contact, and transport coefficients are calculated from equations (3.6-3.12). 3.2.4 Algorithm The equations are integrated using forward-time, centered-space (FTCS) differencing on the regular grid defined above, using variable timesteps. At the beginning of each timestep, the values of n, u, and T for that timestep have been calculated from the previous timestep. The basic algorithm followed by the program is enumerated below and detailed in following subsections: 1. Enforce the boundary conditions (cf Sec. 3.2.5). 2. Calculate P , , , G, , and using eqs. (3.6-3.12). 3. Calculate spatial derivatives from the fields (cf Sec. 3.2.6). 4. Calculate time derivatives using eqs. (3.1-3.3). 5. Calculate the length of the adaptive timestep, and update the main fields for the next timestep (cf Sec. 3.2.8). 6. Add artificial viscosity to the low density regions (cf Sec. 3.2.9). 30 7. Repeat the process starting with enforcing the boundary conditions. The most significant aspects of this algorithm are detailed below. 3.2.5 Boundary Conditions The continuum equations are solved in the reference frame of the container, with the container oscillations accounted for by a sinusoidal external forcing term in the momentum equation, Eq. (3.2). We use impenetrable boundary conditions at the plates: uz = 0 in the frame of the plate. Boundary conditions for horizontal velocities and temperature are set to match those found in the MD simulation. For most of the oscillation, their vertical derivatives at the plates were found to be negligible in the MD simulation, so for simplicity the boundary conditions used for the continuum simulation were: ux /z = uy /z = 0, and T /z = 0, although these derivatives are not identically zero in the MD simulation for the entire oscillation cycle of the plate. For all fields, we use periodic boundary conditions in the horizontal (x and y) directions to eliminate sidewall effects. 3.2.6 Spatial Derivatives We use centered-space derivatives which are accurate to second order in the size of the grid spacing in each direction. Thus for an arbitrary field (x, t), 31 x x0 ,t (x0 + x, t) - (x0 - x, t) + x3 , 2x (y0 + y, t) - (y0 - y, t) + y 3 , 2y (z0 + z, t) - (z0 - z, t) + z 3 . 2z (3.13) y y0 ,t (3.14) z z0 ,t (3.15) This method of taking spatial derivatives is accurate to second-order in grid size x, y, and z. Since we use periodic boundary conditions at the side walls, we are able to implement centered-space derivatives at the side walls as well as in the bulk. At the top and bottom boundaries of the container, however, we do not have periodic boundary conditions, so we use one-sided differencing in the vertical direction: z z0 ,t - (z0 + 2z, t) + 4 (z0 + x, t) - 3 (z0 , t) + z 3 2z (3.16) at the bottom plate, and z z0 ,t (z0 - 2z, t) - 4 (z0 - x, t) + 3 (z0 , t) + z 3 2z (3.17) at the top plate. These one-sided differences retain second-order accuracy in z at the boundaries. 32 3.2.7 Size of Spatial Grid The grid size is selected according to the spatial resolution required by convergence of the quantities of interest, as well as by the resolution in which we wish to visualize the flow (e.g. shock tracking in Chapter 4.) To test the convergence of the solution as a function of changing grid size in the horizontal direction, we test how the solution changes as the size of the grid varies. For our testbed, we use a cell with layer depth at rest of H = 5.4 and size Lx = Ly = 42 in the horizontal directions. Experiments have shown the depth of the layer at rest H = 0.9Np , Lx Ly where Np is the total number of particles in the cell [6]. In continuum simulations we set the total integrated number density in the cell: Np n (x, t) d3 x = Lx Ly H/0.9. We oscillate the cell at = 2.2 and f = 0.4174. This is the same as cell size, depth, and frequency used throughout the majority of Chapter 5 and Chapter 6 to test patterns in oscillated granular layers. For our test, we use an overhead view of the cell, and calculate the height of the center of mass zcm (x, y, t) at each horizontal location (x, y) in the cell at various times t (see Chapter 5 for more discussion of this measurement). We change the horizontal grid size and test how zcm (x, y, t) changes. Since no analytic solution is available, we use a uniform grid with grid spacing x = y = 1 in the horizontal directions, and z = 0.2 in the vertical direction as a reference grid. The simulation was run for 80 cycles until an oscillatory state pattern amplitude was reached (see Chapter 5 for 33 more details of pattern behavior in this system). We then set up square grids of varying size in the horizontal directions: x = y = 1.5, 2, 3, 3.8, and 7, while keeping the vertical grid size constant. These choices correspond to a number of gridpoints in the two horizontal directions Nx = Ny = 28, 21, 14, 11, and 6, respectively. Each new grid was initialized at time t = 0 from the oscillatory state of the reference grid by linear interpolation between gridpoints. Figure 3.1 shows the initial conditions set up for the simulations on each grid. As the grid size is increased, the pictures look visibly coarser, although the basic shape of the waves is maintained except in Fig. 3.1(g), where the pattern is grossly undersampled with only 6 gridpoints across 2 wavelengths in the x- direction. From the initial state, each simulation was run on its own grid for 11 cycles of the plate using the new grid spacing in each case. Pictures from these simulations are shown in Fig. 3.2 exactly 10 cycles of the plate (t = 10/f ) after the initial conditions shown in Fig. 3.1. As can be seen from the picture, for Nx = Ny 11, the simulations remain quite similar. However, by Nx = Ny = 11 (Fig. 3.2(f)), the peaks and valleys look quite washed out, and by Nx = Ny = 6 (Fig. 3.2(g)), the layer looks very flat. It is worth noting that even in the case of Nx = Ny = 6, the layer is thrown off the plate and collides with the plate later in the cycle with a center of mass trajectory similar to the other grid sizes when averaged over the entire cell. However, this grid spacing is not sufficient to capture the peaks and valleys which give the horizontal 34 42 (a) (b) y zcm 9 8 42 (c) (d) 7 6 5 42 (e) (f) 0 0 42 x 42 Figure 3.1: Center of mass height zcm (x, y, t = 0) as a function of horizontal location x, y for a grid with a number of horizontal gridpoints (a) Nx = Ny = 42, (b) Nx = Ny = 28, (c) Nx = Ny = 21, (d) Nx = Ny = 14, (e) Nx = Ny = 11, and (f) Nx = Ny = 6. These images show the initial conditions for the convergence test, based on linear interpolation between gridpoints from the oscillatory state of the Nx = Ny = 42 simulation (a). 35 42 (a) (b) y zcm 9 8 42 (c) (d) 7 6 5 42 (e) (f) 0 0 42 x 42 Figure 3.2: Center of mass height zcm (x, y, t = 10/f ) as a function of horizontal location x, y for a grid with a number of horizontal gridpoints (a) Nx = Ny = 42, (b) Nx = Ny = 28, (c) Nx = Ny = 21, (d) Nx = Ny = 14, (e) Nx = Ny = 11, and (f) Nx = Ny = 6. These images are snapshots of the layer taken 10 cycles of the plate after the images in Fig. 3.1. 36 dependence of zcm . After 10 cycles have passed, the simulations are then compared at 8 equally spaced times over the next cycle of the plate. At each time t0 , zcm (x, y, t0 ) from each different simulations are compared to zcm(ref ) (x, y, t0 ) from the simulations conducted on the reference grid. The simulations are interpolated back onto the same grid; then the root mean square difference between the solution from the new grid and the solution from the old grid is defined as: R (t0 ) = zcm (x, y, t0 ) - zcm(ref ) (x, y, t0 ) 2 where the brack- ets indicate an average over all horizontal locations in the cell at time t0 . This error is then expressed as a percentage of the mean height in the cell at that time: R (t0 ) = 100R (t0 ) / zcm (x, y, t0 ) . Finally, R is averaged over eight equally spaced times in the cycle to give an average percentage deviation from the reference grid, . Figure 3.3 shows as a function of N = Nx Ny for each grid used in the test. Increasing the number of gridpoints on the new grid in each case decreases the percentage error from the reference grid. The grid with N = 21 21 , corresponding to grid spacings of x = y = 2 shows less than a 2.5% root mean square difference from the reference grid. This error was judged to be sufficient for the accuracy needed for this study, so we use horizontal grid spacing 0 < x = y 2 throughout this dissertation. Subject to this constraint, the grid spacing was selected in each case based on considerations of simulation time and resolution required for imaging. A similar test showed that vertical grid spacing z yielded suffi37 20 (a) 15 10 5 0 0 10 2 200 400 600 N 800 (b) 10 1 10 1 10 0 10 2 N 10 3 Figure 3.3: Percent difference from the reference simulation of Nref = 1764 as a function of changing number of gridpoints N = Nx Ny , shown on a linear scale (a) and a log-log scale (b). 38 cient accuracy for pattern simulations. However, much of the work in Chapter 4 and Chapter 5 is based on tracking shock trajectories which propagate quickly in the vertical direction. Thus, we use a finer mesh in the vertical direction, z = 0.2, throughout this dissertation, in order to track shock the vertical location of the shock with high precision. 3.2.8 Forward timestepping We use a fully explicit forward-time stepping scheme which is accurate to first order in the timestep, using adaptive timestep t: n (x, t0 + t) n (x, t0 ) + n t + t2 , x,t0 (3.18) u (x, t0 + t) u (x, t0 ) + u t + t2 , x,t0 (3.19) T (x, t0 + t) T (x, t0 ) + where n , u , t x,t0 t x,t0 T t + t2 , x,t0 (3.20) and T t x,t0 are calculated using (3.1-3.3) and the spa- tial derivatives calculated using our centered-space scheme. This method is accurate to first order in timestep t. Forward-time, centered-space (FTCS) differencing is known to be unstable for advective terms [87]. However, for a simple advective equation, it can stabilized by the presence of a diffusion term D (real or numerical). For 39 the one-dimensional advection-diffusion equation, u 2u u = -v + D 2, t x x FTCS differencing is stable for D = condition: t x/ |v| [87]. Since we are simulating our grains as a viscous fluid, the FTCS scheme may be stable provided we have sufficient viscosity and sufficiently small timesteps. Due to the nonlinearity of the coupled equations, stability criteria may not be obtained analytically as in a simple advection-diffusion equation. To ensure stability, we use adaptive timesteps which obey conditions similar in form to the Courant condition: dt C |xmin | , dt C |xmin | , dt C |xmin | . |ux | |uy | |uz | |x|2 2t (3.21) if the timestep satisfies the Courant The Courant-type condition assures that if velocity becomes too large, the timestep will become small to compensate. The JR equations contain terms besides simple advection terms. As part of our stability criteria, we regulate the system so that the layer does not undergo extremely large changes in a single timestep. Using dimensional analysis, we introduce the following condition to ensure that if diffusion and pressure terms become large, the timestep will become small to prevent extreme changes: min min min dt ns |x8 | , dt ns |x8 | , dt ns |x8 | , 2 2 2 40 dt (1-e2 )G()T , at all locations in the cell where |xmin |2 is the smallest gridspacing between x, y, and z; the coefficients C = 1.5, s = 0.16, and = 0.3 are shown empirically to stabilize the code (the code runs for hundreds of cycles without diverging or becoming numerically unstable), and increasing or decreasing any of these criteria by a factor of two results in negligible changes in the results. Thus, this code is considered to be stable and converged with respect to timestep size when these criteria are used. It is possible that with further testing, even larger timesteps could be used to decrease computation time. Roughly ten cycles of the plate can be simulated in a 48 hour period using a single 2.4 GHz Athlon desktop computer with 1 gigabyte of RAM, using the mesh size, timesteps, and parameters discussed in this chapter. 3.2.9 Artificial Dissipation in Low Density Regions The transport coefficients calculated from the Jenkins-Richman equations are sufficient to produce stable, repeatable results. However, in low density regions, terms of order 1/G() in Eqs. (3.11) and (3.12) become large. These regions have very low densities, and have negligible impact on the solution in the bulk of the layer, which is where we focus our study. For densities near close-packed, on the other hand, terms of order 1/G() are negligible in Eqs. (3.11) and (3.12). Since these terms are only significant in low density regions which do not significantly affect the bulk of the layer, these terms were neglected and artificial dissipation [97] was introduced g 41 in the low-density region above and below the layer. Artificial dissipation (AD) is added to the velocity and temperature fields through an additional dissipation term with coefficient DAD : ux t = DAD 2 ux , AD (3.22) uy t = DAD 2 uy , AD (3.23) uz t = DAD 2 uz , AD (3.24) T t = DAD 2 T. AD (3.25) The strength of the artificial dissipation is calculated as a function of density, so that DAD = where f () = e-200( ) . 6 2x2 f () , 3t (3.26) (3.27) In the bulk of the layer, where the shock forms and propagates (cf Chapter 4) and where the standing wave patterns are observed (cf Chapter 5), the volume fraction is high, and the artificial dissipation has negligible effect on the solution (Fig. 3.4). 42 1 f(n) 0.8 a) b) 0.6 0.4 0.2 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0 0.01 0.02 0.03 0.04 0.05 0.06 Figure 3.4: Functional form of the strength of artificial dissipation f () as a function of volume fraction , a) for the range 0 max , and b) for the range 0 0.1max , where max = 0.65 is the random close-packed volume fraction in 3D. The artificial viscosity is only nonnegligible for volume fractions which are very small compared to the random close-packed volume fraction. Thus, this artificial viscosity has a negligible effect on the shocks and patterns found in the system, while playing an important role in ensuring numerical stability in extremely low density regions. 43 3.3 A Molecular Dynamics Simulation of Granular Shakers We compare the results from our continuum simulations to those from an inelastic hard sphere molecular dynamics simulation that was validated in previous studies [7, 79] by comparison to experimentally observed standing wave patterns for varying control parameters. This MD simulation models particles as inelastic hard spheres, and follows the trajectories of each particle in the container. Collisions between particles are modeled as instantaneous binary collisions which conserve linear momentum, but dissipate energy through inelasticity. As in the continuum model (Sec. 3.1), the dissipation is assumed to be characterized by the normal coefficient of restitution e: e = -vn /vn , (3.28) where vn = (v1 -v2 ) (r1 -r2 )/|r1 -r2 |, and v's and r's are velocity and position vectors for a pair of colliding particles before the collision (no superscript) and after the collision (with a superscript ). Particles are monodisperse, and |r1 - r2 | and |r - r | are identically . 1 2 In previous studies, the normal coefficient of restitution depended on the normal component of the relative colliding velocity of the particles vn [7], such that the collisions become elastic as vn approaches zero: 1 - (1 - e0 ) vn / g e = e0 44 for vn < g, for vn g, 3/4 e (vn ) = (3.29) We use this same model for coefficient of restitution, except in the study of shock formation, where we investigate the effects of velocity-independent restitution coefficient (e = e0 ) and compare to the results from this velocitydependent model (cf Chapter 4). The particles are constrained between the bottom plate which oscillates sinusoidally between height z = 0 and z = 2A, and a ceiling which is fixed at z = 200. These plates use the same collision rules and coefficient of restitution as that used between particles, with the plate treated as a plate of infinite mass. Periodic boundary conditions are used in the horizontal directions x and y. In each case in which MD results are compared to results from continuum simulations, the total mass of the layer matches with that used in continuum simulation. To be consistent with the continuum model, surface friction is neglected between particles as well as between the particles and the plate. The effect of surface friction will be discussed in Chapter 7. Previous simulations using this code but including surface friction have quantitatively reproduced standing wave patterns found experimentally in oscillated layers of lead spheres (see Fig. 3.5.) For more details on the algorithm used for this program, please see Ref. [6]. 45 Experiment Simulation Experiment Simulation (a) =3.00 f*=0.27 (e) =5.00 * f =0.44 (b) =3.00 f*=0.44 (f) =5.79 f*=0.47 (c) =4.00 f*=0.38 (g) =6.00 f*=0.84 (d) =4.00 f*=0.38 (h) =7.00 f*=0.75 Figure 3.5: Standing wave patterns from experiments using lead spheres and from molecular dynamics simulations for different nondimensional frequencies f and accelerational amplitude relative to gravity ; showing (a) square patterns (b) stripes, (c) and (d) alternating phases of hexagons, (e) a flat layer, (f) squares, (g) stripes, and (h) hexagons. Patterns (a)-(e) oscillate at f /2, while patterns at higher (f)-(h) oscillate at f /4. From Ref. [7]. 46 Chapter 4 Shocks in Shaken Granular Layers We use the continuum simulations to study the behavior of shock-like behavior in shaken granular layers, and compare these results to those found in MD simulations of the same system. When the accelerational amplitude of the plate is greater than the acceleration of gravity, the layer leaves the plate during each cycle and collides with the plate later in the cycle. These collisions produce shock-like waves which we investigate using a hydrodynamics approach. The majority of this study has appeared in ref. [10] The MD simulations in this study were performed by Sung Joon Moon. 4.1 Background on Shock Waves The speed of sound is the speed at which an infinitesimal pressure wave propagates through a material with respect to the local velocity of the material itself. The relative Mach number M a between a fluid and an obstacle is the ratio of the the relative speed of the obstacle (with respect to the local mean fluid velocity) to the local speed of sound in that fluid. When fluids encounter an obstacle with M a greater than unity at the point of collision, a compression wave front is formed near the object and steepens to form a shock. 47 In an ideal fluid with no viscosity or heat conduction, the wavefront steepens until the fields develop a mathematical discontinuity between the compressed region and the undisturbed region. In a viscous, conducting fluid, the fields vary continuously, but the front steepens until the shock thickness (the width of the region in which the fields change from the undisturbed values to the compressed values) is on the order of a mean free path. When a supersonic fluid impinges perpendicularly onto a flat, impenetrable plate, a shock front parallel to the plate forms and propagates normal to the plate. With respect to the plate, the normal component of the mean velocity of the fluid is zero at the plate. Far from the plate, however, the undisturbed fluid is still supersonic, flowing towards the plate. The Mach number decreases abruptly from greater than unity to near zero as the plate is approached and the density, temperature, and pressure also change rapidly in the region of the shock front. Two major regimes of the flow are divided by this relatively thin (on the order of a mean free path) shock region: the undisturbed region far from the plate, and the compressed region near the plate. This normal shock propagates back upstream in a direction opposite that of the upstream flow velocity. In a shaken layer, when the maximum acceleration of the plate is greater than g, the layer leaves the plate during each cycle. When the layer is off the plate, it is cooled by inelastic collisions between particles, while the particles are simultaneously accelerated towards the plate by gravity; this leads to a large mean velocity compared to the sound speed (See Sec. 4.5). We show in 48 Sec. 4.5 that when the layer contacts the plate later in the cycle, the Mach number is much greater than unity with respect to the plate, i.e., the flow is supersonic. Here we examine the normal shocks that are formed as the supersonic layer collides with the bottom plate during each oscillation cycle. 4.2 Properties of the Shaken Layer Shocks similar to those presented here were found for a wide range of and f . Here we report results obtained from a particular set of parameter values: = 3 and f = 0.095 g/, where is the diameter of a particle. This would correspond to a frequency of approximately 30 Hz for particles with a diameter of 0.1 mm. The particles have a coefficient of restitution e = 0.9, except in Sec. 4.6, where e is varied to investigate the effect of inelasticity on shock propagation. We examine shock formation and propagation in layers that are approximately 9 particles deep as poured except in Sec. 4.7, where layer depth is varied. Throughout this paper, the phrase "the layer" refers to the dense region in which the volume fraction of particles is greater than 4% of the random close-packed volume fraction, max = 0.65. The number of particles in the container per unit area of the bottom plate is 10/ 2 , except in Sec. 4.7, where the number of particles is varied to investigate shocks in layers of different depths. That is, the layer would be 10 deep if the particles were arranged in a simple cubic lattice. In actual packings seen experimentally, 10 particles/ 2 corresponds to an average depth 49 of approximately H = 9 as poured [7]. For consistency with previous investigations, we use H (the depth of the layer as poured) to express the depth of the layer throughout this paper. The values of and f in this study were examined previously in experiments and MD studies of spatial pattern formation in shallow granular layers [7]. We study simulations exhibiting pattern formation in Chapter 5, and find that shocks may coexist with patterns for these parameters. For the present study, however, we wish to isolate the effects of shock formation and propagation from the effects of standing wave patterns formation. Therefore, in the present simulations, pattern formation was intentionally suppressed by considering a container smaller than one full wavelength in either horizontal direction (with periodic boundary conditions at the side walls); the absence of patterns permits better visualization of the dependence of the fields on height z. 4.3 4.3.1 Simulation Parameters Molecular Dynamics Simulations The particles are constrained between a bottom plate which oscillates sinusoidally between height z = 0 and z = 2A, and an upper plate which is fixed at z = 200. In the two horizontal directions, the width is 20, and the boundary condition is periodic. Except in Sec. 4.7, 3937 particles were used; the number of particles was chosen to correspond to H = 9 and the total mass of the layer matches with that in the continuum simulations. In Sec. 4.7, 50 1969 particles were used to correspond to H = 4.5, and 5906 particles were used to correspond to H = 13.5. To compare directly to continuum simulations, we characterize the dissipation by a single parameter: the normal coefficient of restitution e = 0.9, except in Sec. 4.6, where we test the effects of the velocity-dependent restitution model described in Eq. (3.29). The container has the same e as the particles, and the mass of the container is assumed to be infinitely large compared to that of the granular layer. Fields are calculated by dividing the container height into small bins of size . The number density is the number of particles found in a bin divided by the volume of the bin. The number density n from MD simulations is averaged over the same phase angle of the oscillatory state for ten cycles of the plate. The velocity v of each particle is used to calculate the mean velocity in each bin, u = v , where the brackets represent the average over all the particles found in a bin at the same phase angle of the oscillatory state for ten cycles of the plate. Granular temperature is then defined as T = 1 3 |v - u|2 . Throughout this paper, we use units such that the particles have mass unity. 4.3.2 Continuum Simulations The granular fluid is contained between an impenetrable bottom plate which oscillates sinusoidally between height z = 0 and z = 2A, and an upper plate which oscillates with the bottom plate. We initially conducted the simulation in a cell with a ceiling that oscillates with the bottom plate at a 51 height 200 above the bottom plate. For H = 9, where most of these simulations were conducted, an average of less than one particle per oscillation reached a height of z 80 in the MD simulation. Changing the height of the cell from 200 to 80 in the continuum simulation resulted in a significant speedup in computation time, but did not change the shock profile. The fractional root mean square difference over one cycle between the shock location produced in the tall cell (200 high) and the shorter cell (80 high) was less than 1%. Therefore, for computational efficiency, the continuum simulations were conducted with a ceiling located at a height 80 above the bottom plate, oscillating with the bottom plate. In MD simulations of shallow layers, H 5, some particles reach heights of z > 80, so we used a cell height of 200 in continuum simulations for H 5 (see Sec. 4.7). We use periodic boundary conditions in the horizontal directions. Initial simulations were conducted with a container 20 wide in the horizontal directions. For reasons of computational efficiency, we reduced the horizontal width of the cell to 11, which produced no significant changes in the shock behavior. The width in the horizontal directions is 11 for continuum simulations throughout this chapter. 4.4 4.4.1 Shock Formation and Propagation Tracking the Layer, the Shock, and the Shock Width Both simulations start with a slightly randomized flat layer near the plate which is then oscillated (cf Chapter 3). A periodic state is reached after 52 a transient state of several cycles of the plate. Although most of the particles move together as a high-density layer, some particles are always found above and below the high-density region. These particles are found both in the MD simulation, where with sufficient time, a small number of particles are found even in very low density regions; and in the continuum simulation, which finds non-zero densities throughout the entire container. We define "the layer" as the region in which the volume fraction is greater than 4% of the random close-packed volume fraction max . Below this volume fraction, there are too few particles for accurate averaging in the MD simulation. With each cycle of the plate, a normal shock wave is formed and propagates through the layer. This shock wave separates a compressed region near the plate from the undisturbed region which is still falling towards the plate. Both the undisturbed region and the compressed region vary in space and time. In Sec. 4.4.2 we find that when the shock is in the layer, granular temperature increases with height in the compressed region, drops quickly in the shock, and then increases again in the undisturbed region (Fig. 4.1). The shock is identified in simulations as the region in which the temperature falls from its maximum at the top of the compressed region to the local minimum at the bottom of the undisturbed region. The shock width d is defined as the distance between the locations of this maximum and this minimum. For purposes of finding the propagation speed of the shock, the "shock location" is defined as the steepest region of temperature decrease, which is located near the center of the shock region (shaded regions in Fig. 4.1). 53 4.4.2 Shock Profile During Formation and Propagation During each oscillation, the shock forms and propagates at the same phase angle with respect to the plate's oscillation, so the dynamics of the system are fully characterized by the time interval between time f t = 0 and one cycle later, f t = 1. We now investigate the formation and propagation of the shock by examining the behaviors found in MD and continuum simulations for various times f t during the cycle (Fig. 4.1). 4.4.2.1 f t = 0: The layer falls towards the plate At f t = 0 the container is at its minimum height. The particles, thrown off the plate in the previous cycle, are now falling towards the plate. The temperature is nearly zero for most of the material, increasing with height only in the low density region above the layer (Fig. 4.1). The effects of the artificial dissipation which is used in the continuum simulation are most pronounced at this time in the cycle, while the layer is dilated and above the plate. The shock formed in the previous cycle has propagated out through the top of the layer and is not present in this picture. The difference between the MD and continuum simulation at large heights has little effect on the shock behavior at later times in the cycle. 4.4.2.2 f t = 0.15: Shock forms as the layer strikes plate The layer now begins to strike the rising plate. The layer compresses and particle velocities are randomized by collisions, increasing and T dras- 54 0 40 z/ 30 20 10 0 40 z/ 30 20 10 0 40 z/ 30 20 10 0 40 z/ 30 20 10 0 0 10 T/(g) 20 30 T/(g) ft = 0 plate T/(g) ft = 0.15 shock region plate T/(g) ft = 0.22 plate T/(g) plate ft = 0.40 0.2 0.4 0.6 Figure 4.1: Dimensionless temperature T /(g) and volume fraction as functions of dimensionless height z/ at four times f t in the cycle. For each time, a picture from MD simulation is shown in the left column, with particles colorcoded according to temperature: high T in red, low T in blue, and the bottom plate of the container shaded solid gray. The right column shows horizontally averaged (blue, labeled on the bottom graph) and T /(g) (red, labeled on the top graph) as functions of z/ (ordinate) for the same four times. The plate is shown as a horizontal black solid line, results from MD simulation are dots, and continuum results are solid lines. At times when a shock is present in the graph, the shock region (see Sec. 4.4.1) is shaded in gray. 55 tically near the plate. This results in the formation of a shock region where and T change rapidly from the compressed, high temperature region near the plate to the dilute, low temperature undisturbed region which is still falling towards the plate. At f t = 0.15 the shock is just beginning to form and is broad and not fully developed (see Sec. 4.5). During this time and throughout the formation and propagation of the shock, the shock profile from the continuum model shows remarkable agreement with the profile found in the MD simulation. 4.4.2.3 f t = 0.22: Shock propagates through the layer The layer continues to compress on the plate, resulting in the propagation of the shock up through the layer. Here the shock is fully developed and has propagated through much of the layer, but there is still a significant number of particles falling towards the plate. The shock has steepened since f t = 0.15, with and T changing from the upstream to the downstream values within a distance of about 2. As will be discussed in Sec. 4.5, the width of the shock continues to change throughout the cycle as the shock propagates through parts of the layer with different densities. A discontinuity in the derivatives of and T appears at the leading edge of the shock, showing the boundary of the undisturbed region. Collisions cause the layer to gradually cool behind the shock, creating a lower temperature near the plate. 56 4.4.2.4 f t = 0.40: The layer begins to leave the plate The shock has propagated through the bulk of the layer and now enters the very dilute region above the layer. In this region the temperature becomes uncertain in the MD simulation because there are too few particles to use in averaging. Simultaneously, the plate is approaching its maximum height, and the layer begins to leave the plate as the downward acceleration of the plate becomes larger than g. The layer continues to cool, setting the stage for the next oscillation. 4.4.3 Shock Dynamics The formation and propagation of the shock is shown in Fig. 4.2 for three values of the restitution coefficient e. The effect of varying inelasticity will be discussed in Sec. 4.6. In each case, when the bottom of the falling layer hits the plate, the center of mass of the layer is still falling towards the plate, and the layer is compressed. At this time, a shock is formed near the plate. The material below the shock is visibly compressed as compared to the region above the shock, as can be seen from the shading in Fig. 4.2. The shock propagates up through the center of mass and out into the low density region above the layer. Then the layer leaves the plate, dilates and cools during its free flight, and once again falls towards the plate. As the shock propagates through the layer, it goes through regions of differing and T , and the shock velocity shows some slight variation with height, as seen in the curvature of the shock location (Fig. 4.2). 57 60 z/ (a) e=0.99 shock location 50 40 center of mass 30 20 10 0 (b) e=0.90 50 40 center of mass 30 20 10 0 (c) e=0.85 50 40 center of mass 30 20 10 0 0.2 0 0.2 ft plate shock location plate shock location plate 0.6 0.5 0.4 0.3 0.2 0.1 0.6 0.5 0.4 0.3 0.2 0.1 0.6 0.5 0.4 0.3 0.2 0.1 0.4 0.6 Figure 4.2: Location of the shock (solid lines for continuum, triangles for MD) and the center of mass of the layer (dashed lines for continuum, circles for MD) as a function of time f t during one cycle of the plate (thick solid line) for particles with (a) e = 0.99, (b) e = 0.90, and (c) e = 0.85, starting from the same initial conditions at f t = -0.33. The plot is shaded according to the volume fraction from the continuum simulation, so that high volume fraction is dark and low volume fraction is light. The "top" and the "bottom" of the layer from MD (when the volume fraction drops to less than 4% of random close-packed) are shown as +'s and 's, respectively. 58 4.5 Mach Number and Shock Width A prerequisite for shock formation is that the local Mach number of the flow be greater than unity with respect to the object causing the disturbance. We calculate the speed of sound from MD and continuum simulations using a relation derived from the continuum equation of state Eq. (3.6) [101]: c= where = 1 + 2(1 + e)G(). 2 T 1 + + , 3 (4.1) As the layer falls towards the plate, it is accelerated by gravity and becomes supersonic in the reference frame of the plate. When the layer hits the plate, it is stopped by the impenetrable boundary, and the Mach number falls to nearly zero with respect to the plate (Fig. 4.4). The undisturbed region is still falling towards the plate at supersonic speeds, and the shock is formed at the boundary between the compressed (low M a) region and the undisturbed (high M a) region. The Mach number in the undisturbed region is largest near the shock where the layer has been falling the furthest under gravity, and gradually decreases with increasing height. The Mach number is directly calculated from the values of , u, and T, so the apparent increased discrepancy between the two simulations as compared to the discrepancy in previous figures is due to propagation of errors in these quantities in calculating M a. While the shock is in the layer, the Mach number always reaches its 59 40 z/ 30 10 Ma 20 10 plate 0 0 2 4 6 Ma, 10 8 10 Figure 4.3: Mach number (M a) (dashed line for continuum and triangles for MD) and rescaled volume fraction (10) (solid line for continuum and circles for MD) as functions of z/ when f t = 0.22. At this time the shock is fully developed and is propagating through the layer (see Fig. 4.1). 60 40 z/ 30 / 20 Ma 10 plate 0 0 2 4 6 Ma, / 8 10 Figure 4.4: Mach number (M a) (dashed line for continuum and triangles for MD) and estimated mean free path (/) (solid line for continuum and circles for MD) as functions of z/ when f t = 0.22. At this time the shock is fully developed and is propagating through the layer (see Fig. 4.1). 61 20 15 10 d/ Ma, d/ Ma 5 1 0.1 1 0.2 ft Figure 4.5: The maximum Mach number (M a) (dashed line for continuum and triangles for MD) and the ratio of shock width d to mean free path (solid line for continuum and circles for MD) as functions of time f t. The shock begins forming at f t = 0.12 and leaves the layer at f t = 0.35. maximum at the lower boundary of the undisturbed region. Thus the maximum value of M a at a given time represents the value at the leading edge of the shock. The flows in this system are hypersonic, with M a > 10 just above the plate immediately before the shock is produced (Fig. 4.5). As the shock moves through the layer, the maximum M a decreases since the fastest flow has been stopped by the shock, but the flow just above the shock remains supersonic until after the shock leaves the layer at time f t = 0.35, and M a always changes from supersonic to subsonic as the shock is crossed moving towards the plate. 0.3 62 Since the density of the layer varies in space and time, the mean free path of the particles in the layer changes as a function of height and time. The mean free path may be estimated as [49]: = /8 . n ( 2 ) (4.2) When the shock is fully developed, the mean free path is less than a single particle diameter ( < ) in the compressed region behind the shock (Fig. 4.5). The mean free path increases in the shock region, and continues to increase as height increases away from the shock. From f t = 0.12 to f t = 0.17 the shock is forming but is not fully developed. The dilute bottom of the layer hits the plate first at f t = 0.12, so that while the initial shock is forming, a higher density region is still falling toward the plate. As higher and higher density material hits the plate, the shock width d broadens (Fig. 4.5) as it forms from f t = 0.12 to f t = 0.17. After f t = 0.17, the shock steepens until the fully formed shock is about one mean free path in width, where it remains until the shock leaves the layer at f t = 0.35. After the shock leaves the layer, d remains approximately one mean free path, but the mean free path itself increases in the low density region above the layer, leading to a broadening of the shock front as the shock leaves the layer (compare f t = 0.40 to f t = 0.22 in Fig. 4.1). 63 4.6 Effect of Inelasticity Previous studies have investigated shocks in elastic and inelastic mate- rials in the asymptotic limit of infinite time [43] and in inelastic materials at finite times [61] for a piston moving into a uniform granular gas at a constant speed in the absence of gravity. We wish to investigate the effect of inelasticity on shock dynamics at finite times in the granular shaker system. However, the oscillatory state of the layer is not appropriate for comparing a gas of elastic hard spheres to one of inelastic hard spheres. An elastic gas, which has no mechanism for removing energy, will never reach an oscillatory state since the container injects energy into the layer at each cycle of the plate. In order to investigate the effect of inelasticity on shocks in granular materials, we slightly modify the system to look at the propagation of a single shock from identical initial conditions rather than looking at oscillatory behavior. We run simulations with different coefficients of restitution e, starting from identical initial conditions of a dilated layer falling towards the plate. These initial conditions are taken from the oscillatory state of the system for e = 0.9 using the same parameters as those used in previous sections. At time f t = -0.33, when the layer is near the peak of its flight above the plate, the coefficient of restitution of the particles is suddenly increased or decreased. Changing e significantly changes properties of the shock and the layer following the next collision with the plate (Fig. 4.2). Since the initial conditions are the same for each case and the layer is in free fall, very few differences result from the change in e until the layer collides with the plate. When the 64 layer hits the plate it is compacted and many collisions occur between particles. The lower the coefficient of restitution, the more energy the particles lose through these collisions, and the more compact the layer becomes. In the case of high e, the layer rebounds and dilates quickly after the collision, while in the case of lower e, the layer remains very compact and leaves the plate almost as a solid body. The center of mass of the high e layer flies very high before beginning to fall back towards the plate, while the lower e value leads to a more shallow flight, with the center of mass beginning to fall earlier in the cycle. The shock trajectory exhibits curvature similar to that of the center of mass, remaining mostly straight for high e, but bending downward significantly for low e. After the collision with the plate, we track the propagation of the shock through the layer with a range of new values of e. Starting from the same initial conditions, the shock starts at the same time and height regardless of e. However, the velocity of the shock depends on e, as Fig. 4.6 shows. Although the shock speed changes somewhat as it propagates through the layer (see Fig. 4.2) , the average shock speed vshock increases monotonically with e (see Fig. 4.6). The special case of elastic particles appears to match with the limit e 1; thus there is no qualitative difference between elastic and inelastic gases. Thus far in this chapter, both MD and continuum simulations were based on a collision model in which e was independent of impact velocity. As discussed in Chapter 3, a more realistic hard-sphere collision model could 65 150 continuum MD vshock / g 100 50 0 0.75 0.80 0.85 e 0.90 0.95 1 Figure 4.6: Average dimensionless shock speed vshock / g in the reference frame of the plate as a function of coefficient of restitution e. vshock is calculated as the average speed of the shock from when the shock is formed until it leaves the layer. 66 include a coefficient of restitution e = e (vn ) that depends on the normal component of the relative colliding velocity of the particles, vn , such that the collisions become elastic as vn approaches zero. To test the effect of such a velocity-dependent restituition coefficient, we ran MD simulations including e (vn ) of the same form (Eq. 3.29) as that used in studies of pattern formation in vertically oscillated layers [7], with e0 = 0.9. The velocity-dependence of the coefficient of restitution produces negligible changes in the shock profiles and trajectories, and less than a 1% change in the average velocity of the shock (Fig. 4.7). For the range of restitution values used in this paper (0.8 e 1.0), very few particles collide with low enough relative velocity for the velocity-dependence to significantly affect the outcome. However, for velocity independent e 0.8, the MD simulation was unsuccessful due to inelastic collapse, while velocity dependent e (vn ) (Eq. 3.29) may be used to avoid inelastic collapse with more inelastic particles. 4.7 Shocks in Layers of Different Depths Thus far, we have restricted our investigation of shocks in oscillated layers to layer depths of H = 9. We now investigate how changing H affects shock formation and propagation for particles with e = 0.9. For fluidized layers, increasing the layer depth increases the pressure compressing the layer on the plate during collision with the plate, resulting in tighter packing of the particles near the plate. For 4 H 15, increasing the layer depth increases the volume fraction of the layer near the plate, so 67 40 z/ 30 center of mass 20 10 shock 0 -0.2 0 0.2 ft plate 0.4 0.6 Figure 4.7: Location of the shock ( for constant e = 0.9, for velocitydependent e = e (vn )) and the center of mass of the layer ( for constant e = 0.9, + for velocity-dependent restitution e = e (vn )) as a function of time f t during one cycle of the plate (thick solid line) in the MD simulation. 68 that the trajectory of the center of mass of a layer of depth H = 13.5 is nearly identical to that for H = 9, even though the former has 1.5 times as many particles in the layer (Fig. 4.8(a)). Both trajectories are actually lower than the center of mass trajectory for H = 4.5. For H = 4.5, the layer is almost entirely gaseous, so much so that many particles reach heights greater than z = 80 in each oscillation of the plate in MD simulations with a ceiling at z = 200. Therefore, we use a cell of height 200 for continuum simulations of H = 4.5 (see Sec. 4.2). Changing layer depth also affects the formation and propagation of shocks in the layer. For more dilute layers, the shock forms earlier in the cycle, when the lower edge of the dilute layer collides with the plate. Thus, increasing layer depth causes the shock to form later in the cycle, at which point the shock then propagates faster through the higher densities found in deeper layers (Fig. 4.8(b)). For H = 13.5, the pressure on the layer during collision with the plate is so high that in some small regions near the plate, particles begin to form ordered lattices and the volume fraction increases beyond the random closepacked volume fraction. For deep layers H 15, a significant portion of the layer is compressed into a solid lattice. The equation of state Eq. (3.6) does not allow volume fractions > max , so continuum simulations cannot reproduce densities greater than random close-packed found in H H 15. For 4, the particles do not move as a coherent layer, but rather act as a gas filling the container for the entire cycle of the plate. 69 40 z/ 30 (a) 20 center of mass 10 plate 0 -0.2 (b) 0 0.2 ft 0.4 0.6 z/ 30 shock location 20 10 plate 0 0 0.2 0.4 ft Figure 4.8: (a) The height of the center of mass of the layer, and (b) the location of the shock as a function of time f t for various layer depths H. The trajectory of the center of mass (a) is shown for a complete cycle in the oscillatory state from the MD simulation with H = 4.5 ( ), H = 9 ( ), and H = 13.5 ( ). The center of mass location from continuum simulations agrees with MD to within a root mean square difference of 3% over one cycle in each of the three cases. The shock location (b) is shown for the fraction of the cycle in which the shock is in the layer, from the oscillatory state of H = 4.5 ( for MD, - - for continuum simulation), H = 9 ( for MD, -- for continuum), and H = 13.5 ( for MD, - - - for continuum ). 70 4.8 Discussion Shocks form in the vertically oscillated granular system with each col- lision of the layer with the plate. Properties of these shocks captured in both MD and continuum simulations show good agreement even though the system exhibits a strong time dependence and large spatial gradients. Results from the continuum simulation agree with the MD results even for ranges where the validity of the continuum equations is in question, including when e is significantly less than one (e = 0.8) and when the volume fraction approaches close-packed. These results demonstrate that continuum equations can describe even a complicated three-dimensional, time-dependent granular system. These simulations provide important information about how the properties of shocks in granular layers are affected by the properties of the particles and of the layer itself. Using a velocity-independent coefficient of restitution did not significantly change the behavior of the system when compared to results from the simple velocity-dependent model in Eq. 3.29. In addition, we found that increasing the coefficient of restitution increases the of speed shocks in the layer but causes no singular behavior in the shock velocity as e approaches unity. Finally, we found that deeper layers exhibit denser packing near the plate and higher shock speeds than shallow layers. Finally, understanding these shocks lays the basis for understanding other behavior in shaken layers. In shaken granular layers, the mechanism for fluidization of the granular media is through contact with the shaking plate. In 71 these small cells, the behavior of the layer is dominated by the formation and propagation of shocks. These shocks seem to play a vital role in the process by which oscillations of the plate provide the energy to drive granular flows. In the next chapter, we investigate larger shaken cells in which shock-like behavior coexists with standing wave patterns. 72 Chapter 5 Patterns in Shaken Granular Layers We use continuum simulations to study the behavior of standing wave patterns in shaken granular layers, and compare these patterns to those found in MD simulations of the same system. In both simulations, standing waves form stripe patterns above a critical acceleration of the cell. We compare the wavelengths found in these simulations to a dispersion relation fit to previous experiments and investigate the dynamics of these waves using continuum simulations. The MD simulations in this chapter and in Chapter 6 were conducted by Jennifer Kreft. 5.1 Patterns Near Onset in Shaken Granular Layers Above a critical acceleration of the plate, standing wave patterns form spontaneously in shaken layers. Just above onset, these patterns are subharmonic with respect to the plate, repeating every 2/f , where f is the frequency of oscillation of the plate [75]. Various subharmonic standing wave patterns, including stripe, square, and hexagonal patterns, have been found experimentally, depending on the nondimensional frequency f = f H/g and the nondimensional accelerational amplitude = A (2f )2 /g, where H is the 73 depth of the layer as poured, and g is the acceleration due to gravity [75]. Experiments and MD simulations indicate that inter-particle friction plays an important role in the standing wave patterns. MD simulations with friction between particles have quantitatively reproduced the stripe, square, and hexagonal subharmonic standing wave patterns seen experimentally for a wide range of parameters [7]. However, MD simulations using frictionless particles do not yield stable square or hexagonal patterns, but only yield stripe patterns, and exhibit the onset of patterns at lower than the critical value for frictional particles [80]. This result is consistent with experiments which show that adding graphite to reduce the friction between particles and between particles and the plate can prohibit the formation of square patterns with long range order [39]. Despite the significance of this system, hydrodynamic equations have not yet reproduced the standing wave patterns found in such layers. Here we investigate the onset of ordered standing wave patterns using continuum equations and compare to results from MD simulations. We consider layers of frictionless particles and investigate the onset of stripe patterns in these layers. However, to investigate other patterns such as squares or hexagons, simulations would need to include friction between particles. We will discuss the effects of friction more fully in Chapter 7. 74 5.2 Properties of the shaken layer Patterns are found in experiments, continuum, and MD simulations for a wide range of and f . Here we investigate patterns formed at = 2.2 for a range of frequencies 0.15 f 0.42. In Chapter 6 we vary to investigate the onset of patterns. In experiments without graphite and in MD simulations with frictional particles, = 2.2 is below the critical acceleration needed to produce patterns. However, reducing friction decreases the critical value of [80] and patterns form in simulations of frictionless particles for < 2.2 (cf Chapter 6). We simulate standing wave patterns in layers that are approximately 5.4 deep as poured. The number of particles in the container per unit area of the bottom plate is 6/ 2 . That is, the layer would be 6 deep if the particles were arranged in a simple cubic lattice. In actual packings seen experimentally, 6 particles/ 2 corresponds to an average depth of approximately H = 5.4 as poured [7]. We simulate particles with restitution coefficient e = 0.70 throughout this chapter. 5.3 5.3.1 Simulation Parameters Molecular Dynamics Simulations As in Chapter 4, the particles are constrained between a bottom plate which oscillates sinusoidally between height z = 0 and z = 2A, and an upper plate which is fixed at z = 200. We examine patterns in cells of varying horizontal lengths, so that Lx and Ly are varied throughout this chapter. In 75 each case, 6 Lx Ly / 2 particles are used, corresponding to an average depth of approximately H = 5.4 as poured, regardless of cell size. We use a velocity-dependent coefficient of restitution (Eq. 3.29), with e0 = 0.7. There is no friction between particles, and the top and bottom plates are modeled as frictionless plates of infinite mass. Boundary conditions in the horizontal directions are periodic to avoid sidewall effects. 5.3.2 Continuum Simulations As in the previous chapter, the granular fluid is contained between an impenetrable bottom plate which oscillates sinusoidally between height z = 0 and z = 2A, and an upper plate which oscillates with the bottom plate at a height 80 above the bottom plate. We vary the horizontal dimensions of the cell Lx and Ly , and use periodic boundary conditions in the horizontal directions. In each case, we compare continuum simulations to MD simulations in cells with the same horizontal dimensions, and the total mass of the layer is the same as that used in MD. Throughout this chapter, grid spacing in the horizontal direction x = y = 2 is used, except in sections Sec. 5.6.2 and Sec. 5.7, where we use x = y = for higher resolution pictures. Throughout, the vertical grid size is = .2. The dissipation is characterized by the coefficient of restitution e = 0.7. 76 5.3.3 Visualizing Patterns To visualize peaks and valleys formed by standing wave patterns, we calculate the height of the center of mass of the layer as a function of horizontal location in the cell at various times in the cycle zcm (x, y, t). At a given time t0 and horizontal location (x0 , y0 ), zcm (x0 , y0 , t0 ) is the center of mass of all particles whose horizontal coordinates lies within a bin of size xbin ybin centered at (x0 , y0 ). For continuum simulations, we use the simulation grid size to define the bins: xbin = x and ybin = y. For MD simulations, we use bins of size throughout this chapter to compare to the higher resolution pictures from continuum. Peaks in the pattern correspond to maxima of zcm , and valleys correspond to minima. 5.4 Standing Wave Stripe Patterns in Simulations of Granular Layers Continuum simulations produce standing wave patterns for = 2.2 and f = 0.4174. Alternating peaks and valleys form a stripe pattern which oscillates at f /2 with respect to the plate oscillation; a location in the cell which represents a peak during one cycle will become a valley the next cycle, and then return to a peak on the following cycle. Snapshots of an overhead view of the cell, showing zcm from continuum and MD simulations are shown in Fig. 5.1 for cell of size Lx = Ly = 126 in the horizontal directions. These snapshots are shown at a time in the cycle f t = 0.25, when the plate is at its equilibrium position, moving upward. At 77 this time in the cycle, the layer is off the plate, and the peaks and valleys in the height of the layer can be seen plainly. Both MD (Fig. 5.1(a)) and continuum (Fig. 5.1(b)) simulations produce stripe patterns starting from an artificial flat state (see Sec. 3.2.3). In each case, stripe patterns eventually align parallel to the cell orientation, although whether the stripes align parallel to the x- or y-direction depends on the initial conditions. In each case, a wavelength is selected which fits six times into the linear dimensions of the box, yielding a dominant wavelength of 21 4. Note that due to the finite size of the box, the wavelength selected must fit into the box an integer number of times in the orientation selected. For some cell sizes, stripes are found oriented diagonal to the box orientation depending on the wavelength of the pattern selected. The snapshot from MD simulations looks significantly less smooth to the eye than the continuum simulation, and the stripes appear less straight. These fluctuations may be understood as a result of the finite number of particles in the MD simulation, as opposed to the hydrodynamic equations which treat the material as a continuum. These fluctuations are investigated in more detail in Chapter 6. 5.5 Dispersion Relations in Continuum, MD, and Experiment Previous studies have shown that the wavelength of standing wave pat- terns in shaken granular layers depends on the frequency of the plate os78 126 a) y z cm 7 6 0 126 b) 4 5 3 0 0 x 126 Figure 5.1: An overhead view of a layer of grains, showing the center of mass height zcm as a function of horizontal position (x, y) in a cell with horizontal dimensions Lx Ly = 126 126, from (a) MD simulations, and (b) continuum simulations. Peaks of the layer corresponding to large center of mass height zcm are shown in white; valleys corresponding to low zcm are shown in black. 79 cillation [20, 74, 113]. For a range of layer depths and oscillation frequencies, experimental data for frictional particles near the onset of patterns were found to be fit by the function = 1.0 + 1.1f -1.32 0.03 , where = /H [113]. We investigate the frequency dependence of standing waves in continuum simulations and in MD simulations of frictionless particles. Dimensionless accelerational amplitude = 2.2 was held constant while dimensionless frequency f was varied. Simulations were conducted in a box of horizontal extent Lx = 168 and Ly = 10. This orientation causes stripe patterns to consistently form parallel to the y- axis. The dominant wavelength in each case was calculated from zcm (x, y, t) by finding the wavenumber kx in the x- direction which exhibited the maximum power during one cycle of the oscillatory state of the pattern. Due to the periodic boundary conditions and finite box size, wavelengths must fit in the box an integer number of times. This finite size effect of quantized wavelength yields inherent uncertainty in the wavelength that would be selected in an infinite box. Wavelengths found in continuum and MD simulations are compared to the dispersion relation fit to experimental data: = 1.0 + 1.1f -1.32 0.03 in Fig. 5.2. Investigation is limited to f > 0.15 by the box size, as only two wavelengths fit in the box in continuum simulations at this frequency. Neither simulation produced patterns for this box size for f 0.45. Both simulations agree quite well with the experimental fit throughout the range 0.15 ft 0.45. The dispersion relations indicate very good agreement between contin80 uum simulations and MD simulations. Comparison to the experimental fit shows that both MD and continuum simulations produce wavelengths consistent with experimental results for frictional particles. Finally, these data indicate that friction seems to be unimportant in wavelength selection through this parameter range. 5.6 Dynamics of Standing Wave Patterns We examine the dynamics of the oscillatory state of standing wave patterns. It is difficult to visualize the internal dynamics of 3D patterns in experiment, but MD simulations give us access to position and velocity data for the particles in the system. Density and velocity may be calculated from appropriate averaging over several cycles of the cell [6]. However, with stripe patterns present, many more cycles are needed to obtain sufficient statistics when compared with the shock study in Chapter 4, where there was no horizontal variation and we were able to average over all particles at a given height. Temperature is even more difficult as it is calculated from the distribution of velocities. A previous study was forced to assume no variation along stripes to calculate the temperature from MD simulations [6]. The continuum simulations, however, calculate flow fields such as density, velocity and temperature at all points in the flow directly as a part of the simulations. For a cell of size 126 126 and a cell of size 168 10 in the horizontal directions, the dominant wavelength was = 21 for f = 0.4174 (see Fig. 5.1 and Fig. 5.2). Here we examine the dynamics of waves in a cell 81 16 (a) 14 * 12 10 8 6 4 2 0 0 continuum only MD only continuum & MD experiment 0.1 0.2 0.3 0.4 f* 0.5 Figure 5.2: Dispersion relations for stripes which form perpendicular to the long dimension of cells with horizontal dimensions 168 10. Data for continuum simulations are shown as triangles, MD simulations as circles, and points where continuum and MD simulations yield the same wavelength are shown as squares. In both continuum and MD simulations, the dominant wavelength of the final oscillatory state fits very well to the dispersion relation found in experiments = 1.0 + 1.1f -1.32 0.03 (solid line) [113]. Error bars in both cases are calculated exclusively from discretization due to periodic boundary conditions in a finite size box. 82 2 2 = 42 42 in the horizontal directions, with periodic boundary conditions in the horizontal directions. This size was selected to minimize computational time by using a small box, while also minimizing the effect of boxsize by using an integer multiple of selected in larger cells. We use f = 0.4174 and = 2.2 for the remainder of this chapter. 5.6.1 Wave Patterns in MD Simulations The time-dependence of the standing wave patterns is evident in snapshots of the container at various points in the cycle, viewing the container from the side. In this MD simulation, the stripes orient parallel to the y-axis, with little variation in the y- direction. We view the cell from the side, with the x - z plane visible and particles with various values of y projected into the x - z plane (Fig. 5.3). Particles are shown as open circles, so that particles with different locations in the y-direction (in front or in back of each other in the figure) appear to overlap. We examine the wave patterns at various times in the cycle f t for two cycles of the plate (Fig. 5.3): 5.6.1.1 f t = 0: The layer collides with the plate At f t = 0 the container is at its minimum height. Since the driving frequency and amplitude are different from those studied in Chapter 4, the layer contacts the plate at a different time in the cycle. In this case, the lowest particles in the layer are just beginning to contact the plate at f t = 0. The layer is dilute and two peaks and two valleys are evident. As in Chapter 4, 83 20 y ft=0 ft=1 10 20 ft=0.25 ft=1.25 10 20 ft=0.5 ft=1.5 10 20 ft=0.75 ft=1.75 10 0 0 20 40 20 x 40 Figure 5.3: Views of MD simulations from the side of the cell, looking down the length of stripes oriented parallel to the y-axis at various times f t during two cycles of the plate. All particle locations within the range (0 y 10) are projected into the x-z plane, regardless of their y-coordinate, and particle positions are represented by open blue circles. Particles that appear to overlap are at different y-coordinates (in front of or behind each other). The plate height is represented as a thick solid red line. The full cell is 42 42 in horizontal extent, so that the entire span of the cell is shown in the x-direction. The cell extends to a height of z = 200, but the figures show only z 20 since no particles fly higher than z = 20 during these two cycles. 84 the "top" of the layer is not unambiguous, as evidenced by the individual particles which can be seen separate from the bulk of the layer, especially near the peaks. As in Chapter 4, the phrase "the layer" refers to the dense region in which the volume fraction of particles is greater than 4% of the random close-packed volume fraction, max = 0.65. Averaged over many cycles, the density approaches zero smoothly (cf Chapter 4) both above the layer (as in this picture) and below the layer when the layer is off the plate. 5.6.1.2 f t = 0.25: The layer is compressed on the plate At f t = 0.25, the plate is moving upwards and the layer is being compressed on the plate. The peaks and valleys can still be seen in this picture, although the amplitude has decreased due to compaction on the plate. The part of the layer near the plate is very compact, as evidenced by the large amount of overlap between particles in the projected image. The tip of the peaks are still falling toward the bulk of the layer at this time in the cycle. 5.6.1.3 f t = 0.5: The layer begins to leave the plate At f t = 0.5, the plate is at the peak of its flight and the layer just begins to leave the plate, although the gap between the lowest particles in the cell and the bottom plate is less than one particle diameter. Through its contact with the plate, the layer has become very flat and compact, and peaks and valleys are obscured. 85 5.6.1.4 f t = 0.75: The layer begins to leave the plate At f t = 0.75, the plate is falling away from the layer, leaving a gap between the layer and the plate. New peaks and valleys are forming, and the layer has begun to expand. The peaks and valleys have now reversed locations; a peak now appears where a valley was before the plate collision, and vice versa. 5.6.1.5 1.0 f t < 2.0: The cycle repeats At f t = 1.0, the plate has undergone one full oscillation, and has returned to its lowest point in the cycle. The features of this point in the cycle mimic those found at f t = 0, except that the peaks and valleys have reversed location. Aside from the reversing of peaks and valleys, this next cycle of the plate exhibits the same features as the previous cycle at all times during the cycle. At f t = 1.75, the peaks and valleys once more reverse, setting the stage for another cycle identical to the first, with the peaks and valleys returned to their original location. This demonstrates the subharmonic behavior of the patterns with respect to the plate. 5.6.2 Wave Patterns in Continuum Simulations The continuum simulations display a very similar time-dependence to that found in MD simulations. In this simulation, stripes have again formed parallel to the y-axis. In this case, we examine snapshots of the standing wave patterns at various times in the cycle from a side view of a slice through the 86 middle of the cell y = Ly /2, displaying variation parallel to the x - z plane. The dynamics of the standing waves are visualized in the volume fraction (Fig. 5.4), the horizontal component of velocity ux (Fig. 5.5), and the vertical component of velocity uz (Fig. 5.6). All velocities in these figures are shown in the "laboratory" reference frame in which the particle is moving. For these plots, a higher resolution grid x = y = was used in the horizontal direction for better visualization of the peaks and valleys. At all times in the cycle, the maximum magnitude of the horizontal component of velocity in the y- direction (uy ) is less than 10% of the maximum magnitude of the velocity in the x- direction, indicating little motion parallel to the stripes. We examine the wave patterns at various times in the cycle f t for two cycles of the plate (Fig. 5.4, Fig. 5.5, and Fig. 5.6): 5.6.2.1 f t = 0: The layer collides with the plate At f t = 0 the container is at its minimum height. As in MD simulations, the layer is just beginning to contact the plate at f t = 0. The layer is dilute and two clear peaks and two valleys are evident (Fig. 5.4). The bulk of the layer is falling towards the plate, although in a very thin layer near the plate, the material has been slowed to nearly zero velocity (Fig. 5.6). Near the peaks, the net velocity in the horizontal direction is towards the peaks (Fig. 5.5), reinforcing the strength of the peaks at this time. Compared to other times in the cycle, however, the horizontal velocity is relatively low in magnitude. Above the layer, the horizontal velocity switches direction; how- 87 20 y ft=0 ft=1 10 20 ft=0.25 ft=1.25 10 20 ft=0.5 ft=1.5 10 20 ft=0.75 ft=1.75 10 0 0 20 40 20 x 40 0.1 0.2 0.3 0.4 0.5 0.6 Figure 5.4: Volume fraction from continuum simulations at a slice y = 21 parallel to the x - z plane at various times f t during two cycles of the plate. Empty space ( = 0) is white, random close-packed volume fraction ( = 0.65) is black, and the color increases from white to black through shades of blue as volume fraction increases. The plate height is represented as a thick solid red line. The full cell is 42 42 in horizontal extent; the entire span of the cell is shown in the x-direction. The cell extends to a height of 80 above the lower plate, but the figures show only z 21, since the density is quite low above this height. No particles fly higher than z = 20 during this time in the MD simulations. 88 20 y ft=0 ft=1 10 20 ft=0.25 ft=1.25 10 20 ft=0.5 ft=1.5 10 20 ft=0.75 ft=1.75 10 0 0 20 40 20 x 40 -1.5 -1 -0.5 0 0.5 1 ux/ g 1.5 Figure 5.5: Dimensionless horizontal velocity in the x- direction ux / g from continuum simulations, through a slice y = 21 parallel to the x - z plane for various times f t during two cycles of the plate. Flow to the right ux > 0 is represented by red colors, flow to the left (ux < 0) is represented by blues, and the darkness of the color increases from white (ux = 0) to black as the magnitude of ux increases. The plate height is represented as a thick solid black line. The full cell is 42 42 in horizontal extent; the entire span of the cell is shown in the x-direction. The cell extends to a height of 80 above the lower plate, but the figures show only z 21. 89 20 y ft=0 ft=1 10 20 ft=0.25 10 ft=1.25 20 ft=0.5 10 ft=1.5 20 ft=0.75 10 ft=1.75 0 0 20 40 20 x 40 -5 0 uz/ g 5 Figure 5.6: Dimensionless Vertical velocity uz / g from continuum simulations, through a slice y = 21 parallel to the x - z plane for various times f t during two cycles of the plate. Upwards flow uz > 0 is represented by red colors, downwards flow (uz < 0) is represented by blues, and the darkness of the color increases from white (uz = 0) to black as the magnitude of uz increases. Velocity is shown in a stationary reference frame in which the plate is oscillating sinusoidally between z = 0 and z = 2A. The plate height is represented as a thick solid black line. The full cell is 42 42 in horizontal extent; the entire span of the cell is shown in the x-direction. The cell extends to a height of 80 above the lower plate, but the figures show only z 21. 90 ever, there are too few particles in this area to affect the behavior of the layer itself. Also in the low density region above the layer, some particles are still moving upwards even though the bulk of the layer is falling. 5.6.2.2 f t = 0.25: The layer collides with the plate At f t = 0.25, the plate is moving upwards and the layer is compacting on the plate (Fig. 5.4). The peaks and valleys can still be seen in this picture, although the amplitude has decreased due to compaction on the plate. The layer is very compact near the plate, nearly reaching random close-packed volume fraction (see colorbar). The bulk of the layer is moving upwards with the plate, although the very tip of the peaks have nearly zero velocity, and the low density region above the layer is still moving downwards. (Fig. 5.6). In the last quarter cycle, as the layer has compacted on the plate, the horizontal velocity has switched directions in the layer (Fig. 5.5). At this time, the net flow is away from the peaks and towards the valleys. This could be expected to flatten the layer further, as mass is carried by the flow out of the peaks and into the valleys. 5.6.2.3 f t = 0.5: The layer begins to leave the plate At f t = 0.5, the plate is at the peak of its flight. The layer is very flat and compact, although the volume fraction next to the plate is less than it was one quarter cycle previously (Fig. 5.4). The bulk of the layer is still moving upwards, although it is no longer doing so homogeneously (Fig. 5.6). 91 The maximum vertical velocity is found in the locations which represented valleys at f t = 0. Near the plate, the vertical velocity is still slightly positive, even though the plate is not moving, indicating that the layer is beginning to detach. Some of the low density region above the layer is still falling towards the plate. Finally, the horizontal velocities have grown in magnitude, and they are still moving away from the old peaks and towards the old valleys, even though the layer is mostly flat at this time (Fig. 5.4). This motion acts to produce new peaks in the location of the old valleys. 5.6.2.4 f t = 0.75: The layer is leaving the plate At f t = 0.75, the plate is falling away from the layer, leaving a gap between the layer and the plate (Fig. 5.4). New peaks and valleys are forming, and the layer has begun to expand. The peaks and valleys have now reversed locations that is, a peak now appears where a valley was before the plate collision, and vice versa. The layer is already beginning to fall, although the plate is also moving downwards at its fastest rate at this time in the cycle (Fig. 5.6). Some of the low density region above the layer is still moving upwards at this time. The magnitude of the horizontal velocity has decreased somewhat, but the flow is still towards the new peaks and away from the new valleys, acting to increase the amplitude of the patterns (Fig. 5.5). 92 5.6.2.5 1.0 f t < 2.0: The cycle repeats At f t = 1.0, the plate has undergone one full oscillation, and has returned to its lowest point in the cycle. The features of this point in the cycle mimic those found at f t = 0, except that the peaks and valleys have reversed location, and the direction of the horizontal velocities have reversed at each point in the cell. Aside from the reversing of peaks and valleys, this next cycle of the plate exibits the same features as the previous cycle at all times during the cycle. In addition, the "sloshing" motion evidenced in the horizontal velocities once again reverses, producing peaks and valleys that by f t = 1.75 are back in their original locations f t = 0. This repetitive sloshing back and forth of peaks and valleys demonstrates f /2 behavior of these patterns for = 2.2. 5.7 Shocks and Patterns The previous section demonstrated the existence of patterns in both MD and continuum simulations, as well as the similarity between the dynamics of these patterns. Continuum simulations also demonstrated a sloshing motion in the horizontal velocity which produced the f /2 behavior of the standing wave patterns. This sloshing motion is consistent with that observed in studies of the velocity fields from MD simulations [6]. In addition to these properties of the pattern, however, these simulations demonstrate a dramatic compression of the layer upon contact with the plate, coupled with a sharp boundary between upwards and downwards 93 flow near the boundary between compressed and dilute areas of the cell. This behavior is very similar to that found in shocks in a smaller cell (cf Chapter 4). The dynamics of the temperature field (Fig. 5.7) confirm the presence of shock behavior coexisting with the standing wave patterns. At f t = 0, as the layer is just contacting the plate (Fig. 5.4), a thin region of very high temperature is formed near the plate (Fig. 5.7). As in the shock investigation (Chapter 4), the shock moves very quickly through the layer, so that by f t = 0.25, the high temperature region is already located at the boundary between the compressed layer and the low density region above the layer. At later times in the cycle (f t = 0.5 and f t = 0.75), the high temperature region propagates through the low density region, losing its sharp shock characteristic as it does so (cf Chapter 4). This shock behavior repeats itself through the next cycle of the plate. Thus shock behavior coexists with pattern forming behavior in this system. 5.7.1 Shock Formation and Propagation To investigate the properties of these shocks in more detail, we focus on the time in the cycle in which the shocks form and propagate through the layer. We visualize volume fraction (Fig. 5.8), granular temperature T (Fig. 5.8), and vertical velocity uz (Fig. 5.9), as indicators of shocks in the layer. In addition, we visualize horizontal velocity ux (Fig. 5.10) as an indicator of the sloshing motion which produces standing wave patterns. Finally, we visualize pressure P (Fig. 5.9) and the derivative of P in the x- direction 94 20 y ft=0 10 ft=1 20 ft=0.25 10 ft=1.25 20 ft=0.5 10 ft=1.5 20 ft=0.75 10 ft=1.75 0 0 20 40 20 T/(g) x 40 0 0.5 1 1.5 Figure 5.7: Dimensionless granular temperature T / (g) from continuum simulations, through a slice y = 21 parallel to the x - z plane for various times f t during two cycles of the plate. Large T is white, T = 0 is black, and the color increases from black to white through shades of red as horizontal velocity increases. The plate height is represented as a thick solid green line. The full cell is 42 42 in horizontal extent; the entire span of the cell is shown in the x-direction. The cell extends to a height of 80 above the lower plate, but the figures show only z 21, since the density is quite low above this height. No particles fly higher than z = 20 during this time in the MD simulations. 95 x P = P x (Fig. 5.10) to investigate how the pressure field produced by the shock influences pattern formation. We use these fields to investigate the behavior of the layer during the time 0.875 f t 1.125 in which a single shock is formed and propagates through the layer: 5.7.1.1 f t = 0.875: The layer is off the plate At time f t = 0.875, the plate is near the lowest point of its cycle. At this time, the layer is off the plate and dilute , with a very low granular temperature (Fig. 5.8). The layer is falling towards the plate, although some of the low density region above the plate is still rising, and the pressure throughout the container is low (Fig. 5.9). The horizontal velocities are flowing towards peaks and away from valleys, and there is very little horizontal variation in pressure (Fig. 5.10). 5.7.1.2 f t = 1: The layer collides with the plate At time f t = 1, the plate is at the bottom of its trajectory and the bottom of the layer is colliding with the plate, increasing the density of the layer slightly near the plate. This produces a shock near the plate at the boundary between very hot region right next to the plate and the rest of the layer which is still cold (Fig. 5.8). The flow below the shock has been stopped by the plate and has nearly zero vertical velocity, while the rest of the layer is still falling rapidly towards the plate (Fig. 5.9). A region of high pressure is developing behind the shock, and a small horizontal pressure gradient is 96 developing behind the shock; the pressure gradient points towards the peaks and away from the valleys, in a direction which will oppose the horizontal velocity at this time (Fig. 5.10). The magnitude of the horizontal velocities has decreases slightly compared to f t = 0.975, but there is little change in ux . 5.7.1.3 f t = 1.125: The shock propagates through the layer At time f t = 1.125, the plate is moving upwards and the shock has propagated through most of the layer. Most of the layer is compacted on the plate, although the tips of the peaks are still dilute and above the shock, as can be seen from the volume fraction and temperature fields (Fig. 5.8). Nearly the entire layer is moving upwards with the plate, although the tips of the peaks are still falling towards the plate. The entire part of the layer behind the shock is at a high pressure; peaks and valleys have developed in the pressure profile, reflecting the peaks and valleys in the density (Fig. 5.9). This produces strong horizontal gradients in the pressure, again pointing away from the peaks and toward the valleys. This pressure gradient has begun to reverse the direction of ux in the regions behind the shock, although the horizontal velocities are very small at this time (Fig. 5.10). The horizontal velocity above the shock has not changed significantly; the shock front creates the pressure gradient which reverses the direction of the flow. 97 y 10 ft=0.875 ft=0.875 ft=1 10 ft=1 ft=1.125 10 ft=1.125 ft=1.25 10 ft=1.25 0 20 40 20 x 40 T/(g) 0.1 0.2 0.3 0.4 0.5 0.6 0 0.5 1 1.5 Figure 5.8: Volume fraction (left column) and granular temperature T (right column) from continuum simulations at a slice y = 21 parallel to the x - z plane at times 0.875 f t 1.25 when the shock is forming and propagating through the layer. The colorbar below each column represents the scale for all pictures in that column. The plate is a thick red line in the left column and a thick green line in the right column. 98 20 y ft=0.875 ft=0.875 10 20 ft=1 ft=1 10 20 ft=1.125 ft=1.125 10 20 ft=1.25 ft=1.25 10 0 0 20 40 20 x 40 arb. -5 0 uz/ g 5 0 2 4 6 8 P( units ) 10 Figure 5.9: Vertical velocity uz (left column) and granular pressure P (right column) from continuum simulations at a slice y = 21 parallel to the x - z plane at times 0.875 f t 1.25 when the shock is forming and propagating through the layer. The colorbar below each column represents the scale for all pictures in that column. The plate is a thick black line in the left column and a thick green line in the right column. 99 20 y ft=0.875 ft=0.875 10 20 ft=1 ft=1 10 20 ft=1.125 ft=1.125 10 20 ft=1.25 ft=1.25 10 0 0 20 40 ux/ g 1.5 20 x 40 P(arb. ) units x -1.5 -1 -0.5 0 0.5 1 -1.5 0 1.5 Figure 5.10: Horizontal velocity ux (left column) and the horizontal derivative of granular pressure x P = P (right column) from continuum simulations at x a slice y = 21 parallel to the x - z plane at times 0.875 f t 1.25 when the shock is forming and propagating through the layer. The colorbar below each column represents the scale for all pictures in that column. The plate is a thick black line in both columns. 100 5.7.1.4 f t = 1.25: The shock leaves the layer At time f t = 1.25, the shock has propagated through the layer and is now leaving the layer. The entire layer is compact on the plate, and the amplitude of the peaks and valleys are starting to decrease (Fig. 5.8). The layer itself has cooled rapidly due to inelastic collisions; a portion of the low density region above the layer is quite hot, and the shock front is still visible as the boundary between the hot region and the cold region above it. The entire layer is moving upwards with the plate, and the shock front acts as a boundary between low density regions still falling and low density regions now moving upwards (Fig. 5.9). The pressure is high in the layer, with a very nonuniform horizontal profile. The horizontal velocity has now changed so that there is a strong net flow pointing from high pressure to low pressure away from peaks and towards valleys (Fig. 5.10). This velocity profile sets up the sloshing motion which will flatten the layer and then produce peaks and valleys in their new locations (see Fig. 5.5). 5.8 Discussion Standing wave patterns form in oscillated granular systems above a critical value of . These patterns oscillate subharmonically with the driving frequency of the plate. Continuum simulations and MD simulations without friction both yield stripe patterns for a range of driving frequencies. Both simulations yield patterns with good agreement between the dominant wavelength of patterns for the same frequencies of oscillation. These wavelengths are con101 sistent with dispersion relations found experimentally for frictional particles. These simulations provide important information about the dynamics of standing wave patterns in shaken layers. The peaks and valleys are produced by a sloshing motion indicated by alternating directions of flow in the horizontal direction perpendicular to the stripes [6]. However, the energy needed to fluidize the system and drive these waves is input into the layer from a vertically oscillating plate with no horizontal component of oscillation. The transmission of energy into the system and the creation of horizontal sloshing motion is apparently linked to shock creation and propagation. Each collision with the plate produces a shock which creates a sharp boundary between regions which differ in density, temperature, pressure, and both vertical and horizontal velocities. The sharp pressure gradients produced by the travel of the shock through the layer drive horizontal flow from high pressure to low pressure regions. This behavior creates the sloshing motion which drives the f /2 patterns in the system. Thus standing wave patterns seem to be intrinsically linked to shocks propagation in shaken layers. In the next chapter, we investigate how increasing through a critical value affects the properties of pattern formation in the shaken layer. In addition, we discuss the effect of finite particle number on the onset of patterns in MD simulations. 102 Chapter 6 Onset of Standing Wave Patterns and the Role of Fluctuations We investigate the onset of patterns with increasing in continuum simulations, and explore the effect of changing on patterns formed above this onset. We compare these results to the effect of changing on onset and patterns in MD simulations. Finally, we describe the effect of fluctuations due to finite particle number in the MD simulation and use Swift-Hohenberg theory to investigate how these fluctuations produce differences between MD and continuum simulations. 6.1 Onset of Patterns in Shaken Granular Layers Experiments and simulations have shown that the dimensionless accel- erational amplitude plays an important role in pattern formation in shaken granular layers [7, 75, 76]. For < 1, the maximum acceleration of the plate is less than that of gravity, and the layer rides up and down on the plate as a static pile. As increases past unity, the layer is thrown off the plate at some point in the cycle, and collides with the plate later in the cycle. Initially, the layer remains relatively flat throughout its flight, but as is increased through 103 a critical point c , standing wave patterns form. We investigate the onset of patterns for increasing in MD and continuum simulations. We also investigate how changing affects layer properties both above and below the critical value of onset in these simulations. 6.2 Simulation Parameters We investigate patterns formed as is increased through the range 1.7 2.3. Frequency is held constant at f = 0.4174, except in Sec. 6.6, where frequency is changed to f = 0.25 to explore the relationship between noise and the wavelength of the pattern. As in Chapter 5, we simulate standing wave patterns in layers in which the number of particles in the container per unit area of the bottom plate is 6/ 2 , corresponding to layers that are approximately 5.4 deep as poured. 6.2.1 Molecular Dynamics Simulations As in Chapter 5, the particles are constrained between a bottom plate which oscillates sinusoidally between height z = 0 and z = 2A, and an upper plate which is fixed at z = 200. For f = 0.4174, the wavelength found in larger cells was = 21 (see Fig. 5.2). We use cells of size 42 42 in horizontal extent, such that two wavelengths fit in the box, except in Sec. 6.6, where we explore the relationship between noise and the wavelength of the pattern. The number of particles in the cell is 42,000 except in Sec. 6.6, where it varies with container size. 104 We use a velocity-dependent coefficient of restitution (Eq. 3.29), with e0 = 0.7. There is no friction between particles, and the top and bottom plates are modeled as frictionless plates of infinite mass. Boundary conditions in the horizontal directions are periodic to avoid sidewall effects. 6.2.2 Continuum Simulations As in the previous chapter, the granular fluid is contained between an impenetrable bottom plate which oscillates sinusoidally between height z = 0 and z = 2A, and an upper plate which oscillates with the bottom plate at a height 80 above the bottom plate. For direct comparison to MD simulations, we use cells of the same horizontal extent as those used in MD simulations, and the total mass of the layer is the same as well. Throughout this chapter, grid spacing in the horizontal direction is x = y = 2, and the vertical grid size is z = 0.2. The dissipation is characterized by the coefficient of restitution e = 0.7, and boundary conditions in the horizontal directions are periodic to avoid sidewall effects. 6.3 6.3.1 Onset of Patterns for Increasing Layers Above and Below Onset of Patterns Both continuum and MD simulations exhibit pattern formation above a critical acceleration of the plate (cf Chapter 5). Below the critical value of , standing wave patterns are not observed (Fig. 6.1). For = 2.2, both MD (Fig. 6.1(b)) and continuum (Fig. 6.1(d)) simulations show well defined peaks 105 and valleys which form stripe patterns with two wavelengths fitting in the box. The only difference between this system and that investigated in Sec. 5.4 is the horizontal size of the cell; these patterns look very similar to a section of the patterns formed in the larger cell (Fig. 5.1). Reducing the accelerational amplitude to = 1.9 while keeping all other parameters constant, yields no ordered waves in either MD (Fig. 6.1(a)) or continuum (Fig. 6.1(c)). Thus both continuum and MD simulations appear to have a critical value of somewhere in the range 1.9 c 2.2, such that no patterns are formed for < c , and patterns are formed for > c . 6.3.2 Growth and Ordering of Patterns With Increasing Despite the similarities, differences between MD and continuum simulations are observable. For = 1.9, the continuum simulation yields a very smooth, flat layer (Fig. 6.1(c)), while MD exhibits visible fluctuations (Fig. 6.1(a)). For = 2.2, the continuum simulations produce stripes (Fig. 6.1(d)) which are much smoother than those found in MD simulation (Fig. 6.1(b)). These results duplicate those found in larger cells (see Fig. 5.1). To explore the differences between the two simulations, we investigate the onset of patterns in more detail in continuum simulations and MD simulations separately. 106 42 zcm a) b) 7 6 5 4 0 3 42 c) d) 7 6 5 4 0 0 42 0 42 3 Figure 6.1: An overhead view of the layer of grains, showing the center of mass height zcm (x, y) of the layer as a function of location in the box, for (a) MD simulations with a plate acceleration with respect to gravity = 1.9, (b) MD simulations with = 2.2, (c) continuum simulations with = 1.9, and (d) continuum simulations with = 2.2. Peaks corresponding to large zcm are shown in white, while valleys corresponding to small zcm are shown in black. The grayscale for all four images use the same units as shown by the colorbars on the right of the figure. 107 6.3.3 Characterizing patterns To measure the size of patterns and fluctuations, we examine the nondi- mensional deviation of the height of the center of mass of the layer as a function of horizontal location in the cell from the center of mass height averaged over the cell at that phase in the cycle: (x, y, t) = (zcm (x, y, t) - zcm (x, y, t) ) /, where x and y are the horizontal coordinates, t is the time in the cycle, and the brackets represent an average over all horizontal locations in the cell at time t. Throughout this chapter, we characterize the patterns at a phase one quarter of the way through a sinusoidal oscillation cycle, such that the height of the plate zp = A, and the plate is moving upwards. Using this definition, 2 (t) represents the mean square deviation of the height of the layer from the mean height of the layer at that phase of the plate. 2 is large for layers with high amplitude patterns or fluctuations, and goes to zero as the layer becomes perfectly flat. In addition to 2 , we distinguish between ordered patterns (stripes) and disordered fluctuations by characterizing the long range order of the pattern. To characterize the long range order of the patterns, we first calcu~ late the power spectrum of the pattern: S (kx , ky , t) = (kx , ky , t) , where ~ (kx , ky , t) = Lx 0 Ly 0 2 (x, y, t) e-ikx x e-iky y dxdy. We then transform to polar 2 2 kx + ky , k = tan-1 (ky /kx ) to find S(kr , k , t) coordinates in k space: kr = in the range 0 k < . We integrate radially and time-average to find the angular orientation of the power spectrum: S(k ) = 108 K 0 S (kr , k ) kr dkr , where K = 2xbin . Lx We bin k into 21 bins between k = 0 and k = , and characterize the long range order of the patterns by the fraction of the total integrated power that lies in the bin with the maximum power: Pmax = 0 Smax , S (k ) dk (6.1) where Smax is the integrated power within an angle /21 of the maximum value of S (k ). Thus Pmax is the fraction of the total power that lies within approximately /21 of the angular location of the maximum power. For a perfectly disordered state with equal power in all directions, Pmax would approach 1 21 0.05, while Pmax = 1 for a state with all power in a single bin. Thus Pmax provides a measure of order when stripes form. 6.3.3.1 Onset of patterns in continuum simulations We investigate the onset of patterns in continuum simulations by comparing 2 of standing waves for different values of . Each simulation begins with a flat layer above the plate with small amplitude random fluctuations. The simulation is run until it reaches a periodic state, at which point 2 is calculated as an average over ten cycles of the same phase of the plate. For 1.95, the initial fluctuations decay rapidly until the layer is quite flat, as represented by negligible values of 2 (Fig. 6.2). As increases, there is a sudden onset to large amplitude waves, as seen by the sudden jump in 2 in Fig. 6.2. This onset occurs at the critical value c = 1.955 0.005. For < c , initial fluctuations decay until it the layer is very flat, while for all 109 layers above onset ( > c ) , these waves produce ordered patterns of stripes similar to those in Fig. 6.1(d). 6.3.3.2 Onset of patterns in molecular dynamics simulations We examine the onset of patterns in MD simulations using the same methods as for continuum equations. Fig. 6.2 shows the mean square height deviation 2 as a function of for MD simulations as well as for continuum simulations. For each value of , the simulation was run for 400 cycles of the plate until the layer has reached a periodic state, then and S (k ) are calculated from an average of the next 100 cycles. As in continuum simulations, 2 increases as is increased. Unlike the continuum results, 2 is non-negligible in MD simulations even for < 1.95. There is still a definite increase in the slope of the curve, but it is delayed until > 2.1, and there is no sharp transition from completely flat layer to standing wave patterns as there is in continuum simulations. Although 2 changes relatively smoothly, making it difficult to pinpoint an onset of patterns from the amplitude alone, there is a distinct onset of long range order in the system (Fig. 6.3). For low , the fluctuations are disordered (cf Fig. 6.1(a)), while for higher , standing wave stripe patterns are observed (cf Fig. 6.1(b)). A clear transition from disordered fluctuations to an ordered stripe pattern is demonstrated by the sharp increase in Pmax as crosses the critical value for long range order determined from Fig. 6.3(b) as LR = 2.15 0.01. c 110 0.8 <2> 0.6 MF=1.955 c continuum 0.4 MD 0.2 0 1.7 1.9 2.1 2.3 Figure 6.2: The mean square deviation 2 of the local center of mass height from the average center of mass height of the entire layer as a function of changing accelerational amplitude for MD (triangles) and continuum (circles) simulations. The vertical dotted line represents the onset of stripe patterns in the continuum simulations (c = 1.955). For this figure, 2 was calculated using bins of size xbin = ybin = 2 for both continuum and MD simulations. 111 1.7 (a) 0.8 1.9 2.1 2.3 <2> 0.6 MF=1.965 C 0.4 LR=2.15 C 0.2 0 0.3 (b) MF=1.965 C P max 0.2 0.1 LR=2.15 C 0 0.1 0 0.1 0.2 Figure 6.3: Comparison of MD simulations (triangles) to the Swift Hohenberg model (solid lines) for (a) 2 , and (b) Pmax as a function of control parameter (bottom axis) for SH, and (top axis) for MD. This plot is obtained by using the parameters for SH simulations which give the best fit between MD and SH simulations. This fit yields a noise strength F = (1.2 0.2) 10-2 and a delayed onset of long range order M F = 0.094 corresponding to LR = 2.15 0.01 in c c MD (the vertical dotted line in the figure). The global ordering jumps sharply at M F = 0.10 (corresponding to LR = 2.15 in MD), representing a transition c c to stripe patterns, while 2 increases smoothly through that transition. This fit predicts a mean field onset value of M F = 1.965 0.007, corresponding to c MF = 0 (the vertical dashed line in the figure.) c 112 6.4 Fluctuating Hydrodynamics The MD simulations differ from continuum simulations in a few signifi- cant ways. First, noisy fluctuations are visible below onset in MD simulations (Fig. 6.1(a)), while continuum simulations are flat to within less than a particle diameter below onset (Fig. 6.1(c)). In addition, even above onset of stripes, patterns appear noisy in MD (Fig. 6.1(b)), while continuum simulations yield much smoother stripes (Fig. 6.1(d)). These behaviors are reflected in a difference in 2 between continuum and MD simulations (Fig. 6.2), and a delayed onset of long range order in MD simulations when compared to continuum simulations (Fig. 6.3(b)), where the onset of long range order coincides with the onset of non-negligible 2 . The delayed onset of long range order, coupled with disordered coherent fluctuations below that onset, are features characteristic of noise-driven damped hydrodynamic modes close to a bifurcation. Near the onset of convection patterns in Rayleigh-B nard convection of fluids, fluctuations caused e by thermal noise create deviations from dynamics predicted from the NavierStokes equations without a noise source. These fluctuations are described by the addition of noise terms to the Navier-Stokes equations [66, 123]. Fluctuating hydrodynamics theory accurately describes the dynamics of fluids near the onset of convection [91, 108, 121]. Due to the relatively small number of particles per wavelength in granular systems, fluctuations caused by noise may be quite significant in the granular shaker. Since the hydrodynamic model used in the continuum sim113 ulations does not include a stochastic noise term characteristic of fluctuating hydrodynamics, the differences between the continuum and MD simulations may be consistent with the presence of strong noise in the MD simulations due to the small number of particles per wavelength. We test whether discrepencies between continuum and MD simulations are consistent with fluctuating hydrodynamics by comparing results from MD simulations to those from the Swift-Hohenberg model. 6.4.1 Swift-Hohenberg Simulations The Swift-Hohenberg (SH) model was developed to describe noisedriven phenomena near the onset of long range order in Rayleigh-B nard cone vection [108]. The SH model describes the time evolution of a scalar field (x, t): = t where - 1 + 2 2 - 3 + (x, t) , (6.2) is the bifurcation parameter, and is a stochastic noise term of strength F , such that (x, t) (x , t ) = 2F (x - x ) (t - t ) and (x, t) = 0. In the absence of stochastic noise (F = 0), called the mean field (MF) approximation, there is a sharp onset of stripe patterns with long range order at = MF c = 0 [102, 108]. For F = 0, the nonlinearity acting on the noise will LR c delay the onset of long range (LR) order to a new critical value: delay in onset is characterized by c > 0. The = LR c - MF . c In addition, the presence of noise creates disordered fluctuations below the onset of long range order ( < LR c ). 114 We numerically solve the SH equation Eq. 6.2, using the scheme described in [22], on a 42 42 grid with periodic boundary conditions, using integration timesteps of 0.5. The box is chosen to be size Lx = Ly = 4, so that two wavelengths of the resulting pattern fit in the box, to match MD and continuum simulations. The simulation was allowed to run for 8000 timesteps to reach a final pattern, then 2 and S (k ) were calculated from an average of the next 2000 timesteps. 6.4.2 Comparing Swift-Hohenberg and molecular dynamics simulations To find the strength of the noise and the mean field onset corresponding to our MD simulations, we fit the SH model to the data from MD simulations (Fig. 6.3) by varying three parameters corresponding to F , c , and an overall scale factor, as in [82]. Of the three parameters, only the noise strength F changes the overall shape of the curve. For a given F , the SH simulation is run for a range of -0.2 0.2, and SH and Pmax are calculated from the steady state solution for each value of . 2 Note that SH in SH simulations is found as a function of control parameter -0.2 SH 2 0.2, where in MD simulations, M D is found as a function of control parameter 1.7 2.3. To compare the onset of the SH model to the onset in MD simulations, we define MD = - M F /M F , C C SH where M F is the mean field onset of patterns, comparable to C = 0. How- 115 ever, we do not know a priori the value of M F . C The onset of long range order is used to establish a correspondence between and . For MD simulations, we measure the onset of long range order as the point of sharpest increase in Pmax (Fig. 6.3(b)). In SH simulations, c represents the onset of long range order. We match the single point of c steepest increase of Pmax between the two curves. The measured value SH then predicts the mean field onset M F corresponding to C Once the relationship between and MD in = 0. is determined, the overall scale 2 2 factor for a given F is found by a least squares fit between SH and M D for the range 1.7 LR (see Fig. 6.3(b)). This minimization procedure c gives the best possible fit for a given value of F . This entire procedure is repeated for varying F , minimizing the squared residual R2 = 2 2 ( M D - SH ) /N , where N is the number of bins. The 2 best fit yields an onset of long range order at c =0.94, corresponding to LR = 2.15. Figure 6.3(a) shows < 2 > as a function of for SH simulations, c and as a function of for MD simulations. The fit shows good agreement in 2 below = 0 (Fig. 6.3). Although the parameters are fit only in the range 1.7 LR , agreement is reasonable even for > LR . c c The three parameter fit not only allows for agreement in 2 , but also matches the measure of order Pmax in the SH model to that found in MD simulation (Fig. 6.3(b)). In both MD and SH simulation, below the critical value of long range order, the fluctuations are disordered, leading to a small 116 6 R2 5 4 3 2 1 x 10 -3 0 0 0.01 0.02 0.03 0.04 F 0.05 2 2 Figure 6.4: The squared residual R2 between M D and SH as a function of the noise strength F used in SH simulations. The best least squares fit is given by F = (1.2 0.2) 10-2 . 117 value in Pmax . Crossing the critical value, Pmax jumps up significantly, and the observed patterns are ordered stripes. Below the onset of stripes, when the fluctuations are constantly shifting and changing, there is significant uncertainty in finding the value of Pmax , as seen by the noisy curve on the plot. Above this onset, however, the standing waves produce stable stripes, and Pmax plateaus and remains quite constant, with good agreement between MD and SH simulations. Finally, the mean field onset M F = 1.965 0.007 predicted c by this fit agrees remarkably well with the critical value c = 1.955 0.005 found in our simulations of Navier-Stokes order continuum equations. 6.5 The Effect of Noise on Onset In the MD simulations shown thus far in this chapter, there are only 42, 000 particles in a system which is two wavelengths square in horizontal extent, and roughly one quarter of a wavelength deep. The hydrodynamic equations, however, treat the system as if it were a real continuum. In fluids, thermal noise due to finite particle number is accounted for by fluctuating hydrodynamic models with extra terms added to the Navier-Stokes equations [66, 123]. For this study, we used no stochastic forcing term in the continuum equations, so the only noise in the system is due to numerical methods, and should be negligible. Continuum simulations should be compared to results from MD in the absence of noise. The continuum simulation shows a sharp bifurcation from a flat layer to a layer with ordered stripes, as fits the prediction of the SH model with118 out noise. In addition, SH simulations with noise are able to fit both the dependence of 2 on control parameter , and the delay in long range order Pmax found in MD simulations . Finally, the value c = 1.955 0.005 found in continuum simulations agrees remarkably well with the mean field onset M F = 1.965 0.007 from MD simulations. c This agreement supports the hypothesis that the differences between the global behavior at onset for the continuum and MD simulations are due to the presence of noise in the MD simulation. The strength of the noise is such that it plays a strong role in the overall behavior of the system for a large range both above and below the mean field onset. 6.6 The Effect of Changing Frequency on Noise Strength Throughout this chapter up to this point, investigation was limited to a single frequency, f = 0.4174. We now examine the effect of fluctuations on layers shaken at a different frequency, f = 0.25. 6.6.1 Growth of Patterns and Fluctuations for Varying Frequency We investigate the effect of changing frequency on the onset of patterns in MD simulations. For cells of horizontal extent 168 10, layers shaken with a frequency f = 0.25 form peaks with a dominant wavelength = 42, which is twice the wavelength found for patterns investigated throughout this chapter at f = 0.4174 (see Fig. 5.2). We examine layers shaken at f = 0.25 in cells of size Lx = Ly = 2 = 119 84, while holding layer depth H = 5.4 and restitution coefficient e = 0.70 constant. We vary through the same range 1.7 2.3 investigated for f = 0.4174 in cells of size Lx = Ly = 2 = 42 earlier in this chapter, and compare the two frequencies. For frictional particles, square patterns are formed for f = 0.25, although in the absence of friction, peaks and valleys remain disordered, and no regular square lattice forms in experiments or MD simulations [39, 80]. We examine the onset of patterns in MD simulations using the same methods as in Sec. 6.3.3.2. Fig. 6.5 shows the growth of 2 normalized by the mean center of mass height of the layer squared 2 2 / zcm (zcm - zcm )2 / zcm 2 2 = for MD simulations with frequencies f = 0.25 and f = 0.4174. In each case, the simulation is run for 400 cycles of the plate at each value of until the layer has reached a periodic state, then is calculated from an average of the next 100 cycles. The lower frequency (f = 0.25) exhibits a much sharper jump in 2 than that seen at f = 0.4174. Below this onset, the curve is much flatter for f = 0.25, while at f = 0.4174, the curve increases much more gradually through onset. This sharper onset is consistent with lower noise strength for f = 0.25 than that found for f = 0.4174. In addition, the rapid growth of peaks and valleys occurs at a much smaller value of for f = 0.25, corresponding to an onset even below the mean field onset M F for the larger c cell. 120 0.1 2<2> 2 <zcm > f *=0.25 0.05 * f =0.4174 0 1.7 1.9 2.1 2.3 Figure 6.5: Comparison of growth of 2 normalized by the mean center of mass height of the layer 2 2 / zcm 2 = (zcm - zcm )2 / zcm 2 for MD simulations with f = 0.25 (squares) and f = 0.4174 (triangles). Both frequencies use a layer depth H = 5.4, with a cell of size 2 2 in the horizontal direction, where the dominant wavelength = 21 for f = 0.4174 and = 42 for f = 0.25. 121 6.6.2 Noise Strength Variation with Frequency We follow the same procedure as in Sec. 6.4.2 to fit the data from MD simulation to the Swift-Hohenberg model. Unlike the stripe patterns formed for f = 0.4174, only disordered peaks are found for f = 0.25. Thus the onset of long range order is ill defined in this case. However, this lower frequency exhibits a much sharper onset in the growth of 2 , which is used to find c . The best fit yields a noise term F = (4 3) 10-4 , and a mean field onset of M F = 1.85 0.01. c Note that the noise strength at f = 0.4174 is approximately 30 times larger than the noise strength at f = 0.25. This leads to qualitatively different behavior of 2 near onset, yielding a smoother curve for the higher frequency and a sharper onset for lower frequency. Finally, the onset is barely delayed for the lower frequency, with for f = 0.4174. How can this wide variation in F be explained? The noise strength found in Rayleigh-B nard (RB) systems provides an analogy for this problem. e The noise strength in Rayleigh-B nard systems FRB is given by: e FRB = kB T d ( /)2 2P r q0 = 0 0 Rc q0 T d d 2kB , 0 0 Rc (6.3) c = 0.01 for f = 0.25, as opposed to c = 0.10 where kB is the Boltzmann constant, q0 is the dominant nondimensional wavenumber at onset, P r = d is the Prandtl number, d = /cp is the thermal dif- fusivity, cp is the specific heat at constant pressure, is density, d is the depth 122 of the layer, and 0 , 0 , and Rc are dimensionless numbers characterizing the Rayleigh-B nard system and associated boundary conditions [48, 115]. e A functional form for the noise strength in shaken granular layers has not yet been found as it has for the Rayleigh-B nard system. Finding such a e formula could allow calculation of the noise strength for a variety of parameters for the granular shaker. However, in the granular shaker, , T , , and vary significantly both spatially and temporally throughout the cell during one oscillation of the cell. Thus it is non-trivial to find such a relation. For an estimate of F for the granular system, we naively apply the analogous noise strength from Rayleigh-B nard convection (henceforth referred to e as the "RB noise") to various points in the cycle. We use the functional form FRB = q0 C0 FRB , where C0 = 2kB 0 0 R c is a nondimensional constant (note that granular temperature is defined such that kB = 1 for granular systems and we use units such that the mass of a particle is unity), and FRB = nT cp d is a nondimensional term which contains the entire functional dependence on hydrodynamic fields n,u, and T . We take the depth d = 5.4 to be the layer depth at rest, and q0 represents the dimensionless wavenumber of patterns near onset. We calculate the specific heat cp using a relation derived from the continuum equation of state Eq. (3.6) [101]: cp = 2 3 + , 2 + d d (6.4) where = 1 + 2(1 + e)G(). We then calculate FRB at various times and locations in the cell, from n, T , u calculated from one cycle of the continuum 123 ft 0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 FRB (0.25, f t) 0.0005 0.0002 0.0002 0.0011 0.0089 0.0190 0.0177 0.0214 FRB (0.4174, f t) 0.0351 0.0019 0.0002 0.0007 0.0060 0.0202 0.0414 0.0622 (0.25)/ 42 42 42 42 42 42 42 42 (0.4174)/ 21 21 21 21 21 21 21 21 Table 6.1: Coefficient of RB noise strengths FRB at the densest region of the cell for f = 0.25 (FRB (0.25, f t)), and for f = 0.4174 (FRB (0.4174, f t) are tabulated for various times f t in the cycle (MD simulation). Also shown is the wavelength in units of particle diameter ((f )/) for each frequency. simulation at = 2.2. The value of FRB calculated at the most dense location in the cell at a time f t during an oscillation of frequency f is defined as FRB (f , f t). We examine the simulations in cells of horizontal extent 126 10 previously investigated in Sec. 5.5 to determine FRB (f , f t). After starting from an artificial flate state, the simulations are run for 30 cycles of the plate, then FRB (f , f t) is calculated for the next cycle of the plate. We display FRB at various times f t in the oscillation cycle for both f = 0.25 and f = 0.4174 (Table 6.1). The value of q0 (0.4174) FRB (0.4174, f t) ranges from 2 10-4 to 6 10-2 ; the noise found by fit to MD simulations for the same frequency was F = 1.2 10-2 . In addition, the ratio of RB noise strengths q0 (f2 )FRB (f2 , f t)/q0 (f1 )FRB (f1 , f t) is shown for f2 = 0.4174 and f1 = 0.25 124 ft 0.0000 0.125 0.25 0.375 0.5 0.625 0.75 0.875 q0 (f2 ) q0 (f1 ) 2 2 2 2 2 2 2 2 69.1015 9.1250 0.9085 0.6399 0.6715 1.0639 2.3329 2.9091 FRB (f2 ,f t) FRB (f1 ,f t) 138.2031 18.2499 1.8170 1.2799 1.3430 2.1277 4.6658 5.8181 q0 (f2 )FRB (f2 ,f t) )F q0 (f1 RB (f1 ,f t) Table 6.2: The ratio of the RB noise strength for f2 = 0.4174 to the RB noise strength for f1 = 0.25 at the densest region of the cell q0 (f2 )FRB (f2 , f t)/q0 (f1 )FRB (f1 , f t) is tabulated for various times f t in the cell. Shown separately are the ratio of nondimensional wavenumbers q0 (f2 )/q0 (f1 ), and of the RB dependence on n, u, and T : FRB (f2 , f t)/FRB (f1 , f t). in Table 6.2. The RB ratio varies from approximately 1.3 to 140, depending on the time in the cycle f t. The equivalent ratio of measured noise strengths obtained by fit of the SH equation to MD simulations is 1.2 10-2 /4 10-4 = 30. The contributions to this ratio are shown separately for the ratio of dimensionless wavenumbers and for the functional dependence on the flow fields n, T , , cp , and . The ratio of the wavenumbers above onset is always 2, since the wavelength selected for f = 0.25 is twice that selected for f = 0.4174. This term alone predicts a larger noise strength in the higher frequency, which corresponds to the observed noise strengths, although it does not account for the size of the ratio. The ratio FRB (f2 , f t)/FRB (f1 , f t) varies greatly through- out the cycle, from 0.6 to 69. This is due to the quickly changing flow fields 125 in the problem. This behavior is contrasted to Rayleigh-B nard systems, in e which these terms are usually relatively constant in both space and time. At all times in the cycle, combining these terms into the RayleighB nard formula for noise strength correctly predicts that the noise strength e should be larger for f = 0.4174 than for F = 0.25. However, the RB ratio of noise strengths varies over two orders of magnitude, with the largest ratio measured from the point in the cycle when the layer just begins to collide with the plate, and the lowest ratio measured from the time in the cycle when the layer is off the plate completely. The measured ratio of 30 falls within the range set by the RB noise term. It may be hypothesized that the large value is due to the fact that the velocity fields which create the peaks and valleys are created when the layer contacts the plate, when the largest ratio is found from the RB noise. However, we are not able to make any quantitative predictions about the noise strength. Much more data from MD simulations is needed in order to propose and test a theory for granular shakers of the type used in Rayleigh-B nard systems. e 6.7 Discussion For the parameters under consideration, both MD and continuum sim- ulations of granular materials form stripe patterns of comparable wavelength above a critical value > c , and display no stripes for < c . The onset of stripe patterns in MD simulations is considerably delayed when compared with continuum simulations, and the amplitude of fluctuations below onset is 126 significant in MD, and absent in continuum simulations. These differences can be accounted for by the presence of noise in the system, and Swift Hohenberg simulations allow us to predict the mean field onset from MD simulations, which agrees well with the actual onset in continuum simulations. We find the strength of the noise to be F = (1.2 0.2) 10-2 for stripes at f = 0.4174, and F = (4 3) 10-4 for disordered squares at f = 0.25. The noise strength in fluid systems is dependent on fields such as temperature, density, and viscosity, as well as the wavelength of the pattern formed. In shaker patterns, all of these fields tend to change when oscillation frequency is altered. Finding and testing a functional form of these fields or of the shaker parameters may be possible with extensive further investigation using MD simulations or experiments. As in ordinary fluids, such a functional form for granular systems may well depend on the wavelength of the pattern as well as other variables such as temperature, density, and viscosity. The value discovered experimentally for a slightly shallower granular shaker system in a previous study F = 3.5 10-3 [40] lies within the range of noise values demonstrated in this investigation (F = 4 10-4 to 1.2 10-2 . The smallest noise strength found for granular shakers is comparable to the largest noise strength found thus far in experiments in ordinary fluids, which is F = 7.1 10-4 , while typical values are closer to 10-9 [82,121]. In our case, the noise is strong enough to delay onset of long range patterns by 10% in MD simulation, and influences strongly the behavior of the system well below (more than 20% below) this onset. Noise should therefore play an important 127 role in the behavior of real systems involving granular materials. 128 Chapter 7 Future Work: Modifying the Continuum Model Continuum simulations of the Jenkins-Richman equations are able to capture important aspects of shocks (cf Chapter 4) and pattern formation (cf Chapter 5 and Chapter 6) in oscillated granular materials. However, this study has also shown some important ways in which this model falls short of capturing the entire range of behavior shown experimentally; these problems will need to be addressed in order to use the continuum simulations to directly model experimental systems. Including noisy fluctuations and friction between particles stand out as two important ways in which continuum simulations could be modified to more accurately reflect the behavior of grains in experiments. In this chapter, we briefly outline some possible methods which can be used to incorporate such corrections into continuum equations similar to the Jenkins-Richman equations. 7.1 Granular Fluctuating Hydrodynamics Experiments have shown [40] the existence of fluctuations below the onset of patterns, which may be accounted for by fluctuating hydrodynamic 129 theory. We have presented evidence in Chapter 6 that MD simulations of shaken grains yield results which are consistent with the presence of noise that is absent in continuum simulations, and that the noise strength in granular shakers is large compared to noise found in ordinary fluid convection. Thus while the continuum equations were able to reproduce wavelengths seen experimentally, the onset of patterns is strongly influenced by the presence of noise in MD simulations. This behavior cannot be reproduced in continuum equations without modification; finding a way to include fluctuating hydrodynamics in continuum simulations would be a step towards directly reproducing experimental results. 7.1.1 Further study of noise in the granular shaker For Rayleigh-B nard convection in ordinary fluids, the functional dee pendence of the thermal (kB T ) noise on hydrodynamic variables (Eq. 6.3) is known [48, 115], allowing noise to be included in model equations such as the SH equations with stochastic forcing of the appropriate strength. However, the equivalent dependence is not known for shaken granular layers. We were able to find the required noise strength for a single frequency and layer depth of shaken grains only by extensive MD simulations and comparison to the SH model equations (Chapter 6). To make matters worse, the noise strength appears to have a strong dependence on shaking parameters (noise strength changes by a factor of 30 when frequency changes by a factor of 2). Creating a granular model able to simulate wide parameter ranges would require better 130 understanding of the noise strength for more than just two sets of parameters. Here we propose a method for investigating the strength of noise for granular systems. One procedure for investigating the noise strength in shaken layers would be to empirically find the amount of noise present for various shaker parameters, and attempt to extract a functional dependence from the data. This could be accomplished with either MD simulations or experiments; MD simulations would perhaps be a more productive starting point since the particle properties can be more precisely controlled. We propose the following procedure (a study similar to this procedure is currently being conducted by Jennifer Kreft): 1. Conduct several MD simulations of layers shaken at different frequencies. 2. For each frequency, repeat the basic procedure in Chapter 6: measure pattern and fluctuation properties as is varied near onset, then fit SH simulations to the MD data. Extract the noise strength needed for each frequency. This procedure will yield an empirical curve of noise strength vs. frequency which would be a useful starting point for our investigations. However, the results from this survey would still be quite limited. For instance, if there is a dependence of noise strength on layer depth, this procedure would have to be repeated for different layer depths. This limitation 131 indicates the basic shortcoming of the procedure: the SH model equation would require the fitting procedure to be repeated to find the noise strength for each new parameter. It would be much more useful if we could find a functional dependence of noise strength on the shaker parameters, as has been done for Rayleigh-B nard convection, where the parameters (such as critical e Rayleigh number, temperature and density in the cell) are combined into a single dimensionless number (see Eq. 6.3) for the noise strength, proportional to kB T . We could attempt to extract this information from our basic procedure by the addition of further steps: 3. For each frequency, calculate the hydrodynamic variables n, u, T , and the wavelength of the pattern near onset. 4. Once data has been taken for a wide range of frequencies and layer depths, attempt to find a nondimensional combination of variables which can fit the noise strength for a variety of parameters (similar to the noise strength dependence in Eq. 6.3). If such a functional form could be found, the SH equation (or a similar continuum model equation) could be used to model bifurcations in shaken granular layers without first having to re-fit the noise strength to each new set of parameters. However, there is a fundamental difficulty with applying this method to the granular shaker system. Unlike Rayleigh-B nard systems, e 132 oscillated granular layers exhibit strong spatial and temporal variation of hydrodynamic variables, both above and below the onset of patterns (see the shocks in Chapter 4 for an example of strong variation even in the absence of patterns). Thus it may be difficult to find a single nondimensional number to characterize the noise strength. Where and at what time do you measure the variables? Does the noise depend on average over the cycle, or is there a time or place in the cycle which sets the noise strength for the rest of the cycle? These are questions which are difficult to answer in the granular shaker, where time and spatial dependence are inseparable and it is difficult to isolate the effects of one from the other. 7.1.2 Investigating noise in other systems In addition to applications in model equations such as the SH equation, fluctuations can be incorporated directly into the governing equations of fluid dynamics such as the Navier-Stokes equations [67, 68]. The average noise strengths and correlations may be found as functions of , , and T . I will not attempt to derive equivalent relations for granular systems at this time. However, if such relations could be found, fluctuations could be directly incorporated at the level of the JR equations. The new equations would need to be tested by comparison to MD simulations or experiment. Because of the difficulties involved in shaken granular layers (specifically the time and spatial dependence previously mentioned), other granular systems may be useful testbeds for investigating the effects of noise. Specifi- 133 cally, systems which limit or eliminate either temporal or spatial dependence (or both) would be good candidates for such an investigation. One such possible system would be a granular analogue to a RayleighB nard (RB) system, where a granular material is heated from below. Une fortunately, it is experimentally difficult to enforce a set temperature at a boundary for granular systems. It has been proposed that a layer vibrated with a period small compared to the mean collisions time between particles may be considered analogous to a thermal bath for granular materials [88]. MD simulations, however, have no such difficulty a temperature may be set by randomizing particle velocities upon collision with a boundary [107]. Convection rolls similar to those found in RB systems have been observed experimentally in experiments of high frequency vibrofluidized layers and in MD simulations heated from below [88, 107]. Molecular dynamics simulations of granular layers heated from below may be a good test bed for investigating noise in granular systems, as they are closely analogous to RB convection in ordinary fluids. They have the additional advantage of developing steady-state convection rolls which would eliminate temporal dependence, although spatial gradients will still exist. This system could be a good candidate for testing the effect of fluctuations at both the model equation level (since the appropriate noise for SH simulations is known for RB fluid systems), as well as for fluctuations at the level of the JR governing equations. Finally, both temporal and spatial dependence could be eliminated in 134 MD simulations of a homogeneously heated granular system in the absence of gravity. Such a system has proven useful in investigations of granular materials due to the fact that it is possible to create steady-state solutions at a uniform temperature [8, 120]. In a homogeneously heated simulation it would be possible to set density and temperature and measure noise strength dependence directly as the hydrodynamic variables are systematically varied. One drawback is the possible gap between the behavior in such an artificial system as compared to experimentally realizable systems. 7.2 Adding Friction to the Continuum Model Adding friction to continuum simulations distinguishes itself as a top priority for modeling realistic grains. In this work, our investigation of patterns was limited by the absence of friction in our continuum model. We were limited to studying stripe patterns, because we are unable to get square or hexagon patterns without the presence of friction. In addition, MD simulations show that frictionless particles display some important differences (e.g. onset of patterns at lower ) than particles which include friction [80]. Due to these complications, we have been forced to use MD simulations as an intermediary between continuum simulations and experiment; MD simulations with friction reproduce experimental results well, so we compare continuum results to frictionless MD results to evaluate the effectiveness of the continuum model. Accounting for the effects of friction could allow continuum simulations to be directly compared to experiments and may go a long way 135 towards creating more realistic continuum models. 7.2.1 Continuum Models with Friction The Walton collision model [119] used in our frictional MD simulations models collisions between particles using both their relative normal velocities (which are then modified by the normal coefficient of restitution), and the relative velocities of the particle surfaces at contact (which are then modified by frictional parameters). The normal velocities may be calculated by knowing only the translational velocities of the centers of mass of the particles. The surface velocities, however, include a contribution due to the tangential part of the translational velocities plus a contribution due to the rotational velocities (spin) of the particles. As previously discussed, the Walton model is sufficient to allow accurate modeling of experiments using MD simulations. However, a similar model has not yet been successfully tested for continuum equations. Researchers have attempted to account for friction in continuum theory in various ways. Below we outline the basic approaches which have been used, and discuss advantages and drawbacks of each, and how these issues inform the project of developing continuum models with friction. 7.2.1.1 Including Rotations in Continuum Model The Jenkins-Richman equations were derived from the BoltzmannEnskog equation by using a particle collision model in which the relative post- 136 collisional normal velocities of colliding particles are modified by the normal coefficient of restitution e. However, the tangential velocities are unmodified from the collision rules for elastic particles, and particle rotations are neglected. The starting point for most attempts to model friction in granular continuum theories is to include both particle rotations and modified tangential velocities in the collision model [42, 57, 72]. At the macroscopic level, including particle spin yields additional hydrodynamic variables: the average angular velocity due to particle rotations , and the rotational temperature Tr , accounting for the average energy due to particle rotations. In three dimensions, this approach yields four new partial differential equations (PDE's) based on balance laws for rotational temperature and for the three components of angular momentum. These equations are coupled to the equations for translational motion due to the possibility of energy and momentum being transferred between translational and rotational degrees of freedom through collisions. The downside of this model is that it increases the difficulty of the numerical problem immensely. Integrating Navier-Stokes type equations for compressible, viscous fluids is already a difficult numerical problem for threedimensional flows; the addition of four new coupled partial differential equations makes direct solution to these equations impractical with current computational capabilities. To this date, I know of no attempts to try to calculate the nine coupled PDE's required for this model of granular materials including friction. 137 7.2.1.2 Additional approximations Due to the difficulty involved, nearly all of the previous attempts to add frictional dependence to continuum equations have made additional approximations to simplify the equations. One approach to simplify these equations is to assume that the mean particle spin is equal to half of the vorticity of the mean velocity [42,57] that is, that there is no local collisional supply of angular momentum. Another approach is to consider specific cases in which the ratio between rotational and translational energy is constant [57, 72]. These approximations attempt to reduce the full equations, once derived, back into systems in which the translational motion may be solved uncoupled from rotations. Using these approximations, Jenkins and Zhang have reduced the equations to a model in which the only effect including rotations is a modification to the normal coefficient of restitution. This approach yields an "effective coefficient of restitution" which includes a dependence on the coefficient of friction [57]. This model has the advantage of simplicity it effectively tries to account for rotations based only on viewing friction as an additional mechanism for energy loss. However, MD simulations of patterns have shown that simply changing the coefficient of restitution cannot reproduce stable square and hexagon patterns [80]. Jenkins and Richman tried a different approach which modeled disks in two dimensions [55] by defining a "generalized velocity" which includes particle spin as a component of velocity, so that only one additional equation 138 is needed to account for particle rotation. However, this approach yielded two distinct theories: one for nearly frictionless particles, and a separate one for very frictional particles which nearly reverse rotations during collision. In a realistic system, either behavior may be possible depending on relative colliding velocities of the particles. Thus the incorporation of friction into continuum equations presents difficulties models have been derived which account for friction at the particle level but which are very difficult to solve; meanwhile, assumptions which attempt to reduce the difficulty of the equations may or may not be relevant to physical flows of grains. 7.2.1.3 A possible method for modeling friction These attempts highlight an important question: "Are there approximations which may be used to simplify frictional continuum equations, while maintaining the important physics of granular materials?" One possibility is to include the effects of friction resulting from relative translational motion between colliding particles, but to neglect the effects of particle spin. The JR equations neglect friction entirely, only including the coefficient of restitution, but no frictional parameters. The various models discussed in this section try to incorporate friction by including the rotation of particles in the particle collision. However, MD simulations of pattern formation in oscillated granular layers show that the amount of energy in particle rotations is on the average two orders of magnitude smaller than the amount of translational 139 energy [78,80]. Thus the effect of rotations may be negligible on average when compared to the effect of translational motion. On the other hand, the relative tangential translational motion (which contributes to the frictional terms in the Walton model) and the relative normal translational motion (where the coefficient of restitution is used) may both be calculated from the relative translational motion between two colliding particles. A collision model could incorporate friction between particles by only including translational components when calculating the relative surface velocities of colliding particle, but still neglecting the rotational motion of the particles when calculating the relative surface velocities. In this case, continuum equations which account for energy loss due to friction as the particles slide past each other may be able to reproduce experiments without requiring the extra equations necessary to model rotations [78]. 7.2.2 Proposed procedure for incorporating friction A possible procedure for testing the above hypothesis is as follows: 1. Simulate oscillated granular layers with MD simulations including surface friction between particles, but neglecting particle rotations entirely (that is, calculate frictional effects only based on the relative tangential translational motion of the two colliding particles). Test whether these simulations can accurately reproduce square and hexagon patterns seen in experiment (the normal coefficient of restitution and the coefficients 140 of friction are already treated as fit parameters for comparison to experiment these may need to be re-fit to account for additional dissipation). 2. If step one is successful in reproducing experimental results, derive new continuum equations from the Enskog-Boltzmann equation, including the effect of friction in calculating tangential velocity changes during collisions, but neglecting particle rotations. This should result in five equations similar to the JR equations, but with additional terms reflecting the effect of friction on the tangential velocities. 3. Numerically solve the new continuum equations (or add the new terms to the existing computer program) to simulate pattern formation in vertically oscillated granular materials. Test whether the new equations are able to capture stable square and hexagonal patterns. This procedure is not guaranteed to work. However, if it is possible to account for the dominant effects of friction without including particle rotations, the new continuum equations could be much more powerful tools for modeling real experiments, while remaining not significantly more difficult than the current JR equations. 141 Chapter 8 Conclusions and Future Work Continuum simulations of granular systems were conducted to investigate various phenomena in vertically oscillated granular materials. These simulations were compared with hard-sphere molecular dynamics simulations to shed light on the potential of such continuum models to accurately reproduce phenomena in granular materials. These simulations indicate great potential for such models; they were able to accurately reproduce many aspects of flows seen in MD simulations, and were able to reproduce experimentally observed pattern wavelengths, despite the lack of an accurate model of friction between particles. The existence of a tested, accurate theory of granular hydrodynamics would open up great new opportunities for investigation. Such a theory could be used to develop powerful new techniques for analyzing and modeling granular flows. In addition, it could provide a new testing ground for a study of hydrodynamics in novel situations or in situations which are experimentally difficult for elastic fluids, such as investigations of high Mach number flows, inelastic and frictional effects, and systems far from equilibrium. This research represents a first attempt to test the validity of hydrodynamics-type 142 equations in time-dependent, three-dimensional inhomogeneous rapid flows of granular materials. The results of these simulations highlight the potential for such models, while also illustrating shortcomings of the current model and suggesting important questions and future research. 8.1 Discussion of major results These findings represent progress in ascertaining the applicability of continuum theory to granular materials. In particular they show that important aspects of rapid granular flows may be captured by a hydrodynamics approach, even in inhomogeneous, time-dependent flows in three-dimensions. Quantitative comparison to experiment is possible for some aspects of these systems (e.g. pattern wavelengths) despite known shortcomings (no friction or stochastic forcing) of our particular model. Our simulations of the JenkinsRichman equations have demonstrated the ability to reproduce important aspects of experimental investigations of vertically oscillated granular materials, including: The Jenkins-Richman simulations qualitatively reproduce the behavior of shaken granular layers as a function of time. The layers leave the plate during an oscillation cycle and collide with the plate later in the cycle, resulting in compression near the plate. This basic dynamic which is seen experimentally in layers oscillated at > 1 is captured in continuum simulations. 143 The continuum simulations produce standing wave patterns above a critical onset of the plate acceleration. These waves oscillate subharmonically with respect to the plate oscillations, reproducing common experimental findings near the onset of patterns. The patterns which develop in continuum simulations exhibit wavelengths which are consistent with a dispersion relation previously measured in experiment. This is true despite the lack of friction in the continuum model. Continuum simulations also quantitatively reproduced results from frictionless MD simulations in cases where direct comparison to experiment was not carried out. In these cases, MD simulations act as a kind of bridge between frictionless continuum simulations and experiments (in which particles have frictional properties). Frictional MD simulations reproduce experimental results well for oscillating granular layers, so when continuum simulations can reproduce results from frictionless MD simulations, it indicates that addition of friction to continuum simulations may allow comparisons directly to experiment. Comparisons to MD simulations include: Continuum simulations and MD simulations produce shock profiles which compare well in density, temperature, and velocity fields throughout the cell. Shocks in continuum simulations quantitatively reproduce time dependence found in MD simulations; the center of mass of the layer and shocks 144 found in the layer follow trajectories which reproduce those found in MD simulations. Continuum simulations reproduce stripe patterns for the same frequencies as MD simulations. Their wavelengths compare well (as well as being consistent with experiment, as noted above). The presence of noise in MD simulations means that the onset of patterns and the amplitude of fluctuations below onset differs between continuum and MD simulations. However, the onset of patterns in continuum simulations are consistent with the mean field onset of patterns for the noise strengths found in MD simulations. A hydrodynamics approach also has promise for extending our understanding of granular materials. Throughout this paper, continuum simulations were used to investigate the physics of granular materials, and in addition, continuum variables n, u, and T were used to analyze shocks, even in results from MD simulations. This approach yielded important insight into the physics of oscillating granular layers: Shocks form in vertically oscillated granular layers with each collision of the layer with the plate. Properties of the particles affect the shock dynamics: increasing the coefficient of restitution increases the velocity of the shock. Similarities between shocks in granular layers and in ordinary gases are highlighted by the fact that the limit of elastic particles (e = 1) 145 is not singular. Finally, deeper layers exhibit denser packing near the plate and higher shock speeds than shallow layers. Continuum and MD simulations exhibit standing wave patterns above a critical acceleration of the plate. These patterns demonstrate the same wavelength as those found in experiments with frictional particles, suggesting that friction does not play significant a role in wavelength selection. The presence of shocks in the system creates pressure gradients which drive the sloshing motion of patterns in the system. Differences between the onset of patterns in continuum and MD simulations are consistent with the presence of noisy fluctuations in MD simulations which are absent in hydrodynamic simulations. We use the Swift-Hohenberg model equation to find the strength of the fluctuations to be orders of magnitude larger than those often found in convection in ordinary fluids. Finally, despite these successes, this study has also highlighted some shortcomings of the Jenkins-Richman model: The absence of friction between particles limits the investigation of patterns in the continuum simulations. It may be assumed that friction plays an important role in other aspects of rapid granular flows as well. The lack of a granular fluctuating hydrodynamic model means that noisy fluctuations are not accounted for in present continuum theory. This 146 limits the ability of granular hydrodynamics to model certain aspects of granular patterns near onset. These two shortcomings immediately suggest directions for future work incorporating modifications to the continuum equations designed to address these problems. I have outlined a proposal for such investigations in Chapter 7. The potential shown by these equations also suggests that they could be very useful in understanding other aspects of granular flows, including: Vibrofluidized granular materials could be studied in more detail using continuum equations. Investigations of the stability of the initial shock front in vibrated granular layers could shed light on the mechanism for initial wavelength selection in the standing wave patterns. With appropriate friction models, other patterns such as squares, hexagons, and f /4 patterns could be investigated. Granular gases at high plate accelerations exhibit steady-state regions which would be interesting to explore using continuum theory. Continuum equations could be used to investigate other systems, including convection in layers vibrated at high frequency (the Rayleigh-B nard e analogue), sheared flows (such as granular Couette flows), or horizontally shaken layers. Some progress has already been made in the Center for Nonlinear Dynamics using hydrodynamic approaches to investigate shocks formed when grains fall onto a solid wedge, and to investigate 147 wake flows created by a cylinder dragged through a fluidized granular layer [46, 93]. Appropriate simplifications to the full continuum equations could yield powerful tools for understanding granular systems. For instance, granular Rayleigh-B nard convection could be studied using linear stability e analysis derived from equations similar to the Jenkins-Richman equation. The Swift-Hohenberg equation has already demonstrated the potential for model equations to illuminate granular flows. Amplitude equations or other models could be derived specifically for granular materials from knowledge of the full governing equations. There are many possible extensions of continuum theory to granular materials. Many granular systems, such as deep vibrated layers or granular Couette flow, exhibit transitions between solid-like and fluid-like behavior. Continuum equations could be used to investigate transition between these granular phases [70]. Two-fluid models are already used to model grains interacting with an interstitial fluid [34]. And continuum equations could be extended to polydisperse media as well as to collections of identical particles. Some of these projects are ongoing, some of them are mere possibilities to be counted among many more possible projects which will bear fruit in the future. There is much we still do not know about the behaviors of granular materials in general, and specifically about granular hydrodynamics. However, 148 one thing is clear there is much exciting work yet to be done and much we can learn from this dynamic field. 149 Bibliography [1] G. Ahmadi and M. Shahinpoor. A kinetic model for rapid flows of granular materials. Int. J. Non-linear Mechanics, 19(2):177 186, 1983. [2] M.P. Allen and D. J. Tildesley. Computer Simulation of Liquids. Oxford University Press, Oxford, 1987. [3] Keiko M. Aoki and Tetsuo Akiyama. Simulation studies of pressure and density wave propagations in vertically vibrated beds of granules. Phys. Rev. E, 52:3288, 1995. [4] Robert P. Behringer, Eric Cl ment, Junfei Geng, Dan Howell, Ljubinko e Kondic, Guy Metcalfe, Corey O'Hern, Guillaume Reydellet, Sarath Tennakoon, Loic Vanel, and Christian Veje. Science in the sandbox: Fluctuations, friction, and instabilities. In D. Reguera, L. L. Bonilla, and J.M. Rub editors, Coherent Structures in Complex Systems, page 351. i, Springer-Verlag, 2001. [5] E. Ben-Naim, S. Y. Chen, G. D. Doolen, and S. Redner. dynamics of inelastic gases. Phys. Rev. Lett., 83:4069, 1999. [6] C. Bizon. Simulations of Wave patterns in Oscillated Granular Media. Ph.D. thesis, University of Texas at Austin, 1998. Shocklike 150 [7] C. Bizon, M. D. Shattuck, J. B. Swift, W. D. McCormick, and H. L. Swinney. Patterns in 3D vertically oscillated granular layers: Simulation and experiment. Phys. Rev. Lett., 80:57, 1998. [8] C. Bizon, M.D. Shattuck, J.B. Swift, and H.L. Swinney. Transport coefficients for granular media from molecular dynamics simulations. Phys. Rev. E, 60:4340, 1999. [9] L. Bocquet, D. Schalk W. Losert, T.C. Lubensky, and J.P. Gollub. Granular shear flow dynamics and forces:experiment and continuum theory. Phys. Rev. E, 65:011307, 2001. [10] J. Bougie, S.J. Moon, J.B. Swift, and H.L. Swinney. Shocks in vertically oscillated granular layers. Phys. Rev. E, 66:051301, 2002. [11] J. J. Brey, M. J. Ruiz-Montery, and F. Moreno. Hydrodynamics of an open vibrated granular system. Phys. Rev. E, 63:061305, 2001. [12] J. Javier Brey, James W. Dufty, and Andr s Santos. Dissipative dye namics for hard spheres. J. Stat. Phys., 87(5/6):1051, 1997. [13] J.J. Brey, M.J. Ruiz-Montero, and F. Moreno. Boundary conditions and normal state for a vibrated granular fluid. Phys. Rev. E, 62:5339, 2000. [14] R.L. Brown. The fundamental principles of segregation. J. Inst. Fuel, 13:15, 1939. 151 [15] N. Burtally, P.J. King, and M.R. Swift. Spontaneous air-driven separation in vertically vibrated fine granular mixtures. Science, 295:1877, 2002. [16] Ludwig Boltzmann (Translated by S.G. Brush). Lectures on Gas Theory. University of California Presss, Berkely, 1964 (originally published in German in 1896 & 1898). [17] C. S. Campbell. Rapid granular flows. Annu. Rev. Fluid Mech., 22:57, 1990. [18] E. Cerda, F. Melo, and S. Rica. Model for subharmonic waves in granular materials. Phys. Rev. Lett., 79(23):4570 4573, 1998. [19] Giovanni Cicotti, Daan Frenkel, and Ian R. McDonald. Simulation of Liquids and Solids: Molecular Dynamics and Monte Carlo Methods in Statistical Mechanics. North-Holland Physics Publishing, The Netherlands, 1987. [20] Eric Cl ment, Loic Vanel, Jean Rajchenbach, and Jaques Duran. Pate tern formation in a vibrated granular layer. 1996. [21] R. Courant and K. O. Friedrichs. Supersonic flow and shock waves. Phys. Rev. E, 53:2972, Interscience Publishers, Inc., New York, 1948. [22] M.C. Cross, D. Meiron, and Y. Tu. investigation. Chaos, 4:607, 1994. 152 Chaotic domains: a numerical [23] I.G. Currie. Fundamental mechanics of fluids. McGraw-Hill, Inc., New York, 1974. [24] S. Douady, S. Fauve, and O. Thual. Oscillatory phase modulation of parametrically forced surface waves. 1989. [25] James W. Dufty. Kinetic theory and hydrodynamics for rapid granular flow a perspective. In Thomas Halsey and Anita Mehta, editors, Europhys. Lett., 10(4):309 315, Challenges in Granular Physics. World Scientific (also available at condmat/0108444), 2002. [26] J. Duran. Sands, Powders, and Grains: An Introduction to the Physics of Granular Materials. Springer-Verlag, Berlin, 2000. [27] Jens Eggers and Hermann Riecke. Continuum description of vibrated sand. Phys. Rev. E, 59(4):4476, 1999. [28] M. M. El-Kaissy and G. M. Homsy. Instability waves and the origin of bubble in fluidized beds, part 1: Experiments. Int. J. Multuphase Flow, 2:379, 1976. [29] P. Evesque and J. Rajchenbach. 307:223, 1988. [30] P. Evesque and J. Rajchenbach. Instability in a sand heap. Phys. Rev. Lett., 62:44 46, 1989. C. R. Acad. Sci. Ser. II (Paris), 153 [31] E. Falcon, R. Wunenburger, P Evesque, S. Fauve, C. Chabot, Y. Garrabos, and D. Beysens. Cluster formation in a granular medium fluidized by vibrations in low gravity. Phys. Rev. Lett., 83:440, 1999. [32] M. Faraday. On a peculiar class of acoustical figures; and on certain forms assumed by groups of particles upon vibrating elastic surfaces. Phil. Trans. R. Soc. London, 52:299, 1831. [33] S. Fauve, S. Douady, and C. Laroche. Collective behaviours of granular masses under vertical vibrations. J. Phys. Colloque (Paris), 50 C3:187 191, 1989. [34] I.K. Gamwo, Yee Soong, and R.W. Lyczkowski. Numerical simulation and experimental validation of solids flows in a bubbling fluidized bed. Powder Technol. (Switzerland), 103:117, 1999. [35] D. Geldart. Types of gas fluidization. Powder Technology, 9:285, 1973. [36] I. Goldhirsch. Scales and kinetics of granular flows. Chaos, 9:659, 1999. [37] I. Goldhirsch. Rapid granular flows. Annu. Rev. Fluid. Mech., 35:267 293, 2003. [38] I. Goldhirsch and G. Zanetti. Clustering instability in granular gases. Phys. Rev. Lett., 70:1619, 1993. [39] Daniel I. Goldman, M. D. Shattuck, Sung Joon Moon, J. B. Swift, and Harry L. Swinney. Lattice dynamics and melting of a nonequilibrium pattern. Phys. Rev. Lett., 90(104302), 2003. 154 [40] Daniel I. Goldman, J. B. Swift, and H.L. Swinney. Noise, coherent fluctuations, and the onset of order in an oscillated granular fluid. Phys. Rev. Lett., 92, 2004. [41] A. Goldshtein, A. Alexeev, and M. Shapiro. Shock waves in granular gases. In T. Poschel and N. Brilliantov, editors, Granular Gas Dynamics, page 188. Springer-Verlag, 2003. [42] A. Goldshtein and M. Shapiro. Mechanics of collisional motion of granular materials. Part 1. General hydrodynamic equations. J. Fluid Mech., 282:75, 1995. [43] A. Goldshtein, M. Shapiro, and C. Gutfinger. Mechanics of collisional motion of granular materials. Part 4. Expansion wave. J. Fluid Mech., 316:29, 1996. [44] A. Goldshtein, M. Shapiro, L. Moldavsky, and M. Fichman. Mechanics of collisional motion of granular materials. Part 2. Wave propagation through vibrofluidized granular layers. J. Fluid Mech., 287:349, 1995. [45] P. K. Haff. Grain flow as a fluid-mechanical phenomenon. J. Fluid Mech., 134:401, 1983. [46] Patrick Heil, E.C. Rericha, and Daniel I. Goldmanand Harry L. Swinney. Mach cone in a shallow granular fluid. Submitted to Phys. Rev. E, 2004. 155 [47] H. J. Hermann and Stefan Luding. Review article: Modeling granular media with the computer. Continuum Mechanics and Thermodynamics, 10:189, 1998. [48] P.C. Hohenberg and J.B. Swift. Effects of additive noise at the onset of Rayleigh-B nard convection. Phys. Rev. A, 46:4773, 1992. e [49] K. Huang. Statistical Mechanics, 2nd Edition. John Wiley and Sons, New York, 1987. [50] H. M. Jaeger, S. R. Nagel, and R.P. Behringer. Granular solids, liquids, and gases. Rev. Mod. Phys., 68:1259, 1996. [51] H. M. Jaeger, S. R. Nagel, and R.P. Behringer. The physics of granular materials. Physics Today, 49:32, 1996. [52] H.M. Jaeger and S.E. Nagel. Physics of the granular state. Science, 255:1523, 1992. [53] J. T. Jenkins. Boundary conditions for rapid granular flow: Flat, frictional walls. Transactions of the ASME, 59:120, 1992. [54] J. T. Jenkins. Kinetic theory for nearly elastic spheres. In H.J. Herrmann, J.-P. Hovi, and S. Luding, editors, Physics of Dry Granular Media: proceedings of the NATO Advanced Study Institute, page 353, 1998. 156 [55] J. T. Jenkins and M. W. Richman. Kinetic theory for plane flows of a dense gas of identical, rough, inelastic, circular disks. Physics of Fluids, 28:3485, 1985. [56] J. T. Jenkins and S. B. Savage. A theory for the rapid flow of identical, smooth, nearly elastic particles. J. Fluid Mech., 130:187, 1983. [57] J. T. Jenkins and C. Zhang. Kinetic theory for identical, frictional, nearly elastic spheres. Phys. of Fluids, 14:1228, 2002. [58] James T. Jenkins and Michel Y. Louge. On the flux of fluctuational energy in a collisional grain flow at a flat, frictional wall. Phys. Fluids, 9(10):2835, 1997. [59] J.T. Jenkins and M.W. Richman. Grad's 13-moment system for a dense gas of inelastic particles. Arch. Rat. Mech. Anal., 87:355, 1985. [60] L. P. Kadanoff. Built upon sand: theoretical ideas inspired by granular flows. Rev. Mod. Phys., 71:435, 1999. [61] V. Kamenetsky, A. Goldshtein, M. Shapiro, and D. Degani. Evolution of a shock wave in a granular gas. Phys. of Fluids, 12:3036, 2000. [62] J. B. Knight, E. E. Ehrichs, V. Y. Kuperman, J. K. Flint, H. M. Jaeger, and S. R. Nagel. An experimental study of granular convection. Phys. Rev. E, 54:5726, 1996. [63] W. Kroll. Chemie Ingenieur Technik, 27:33, 1955. 157 [64] A. Kudrolli and J. P. Gollub. Patterns and spatiotemporal chaos in parametrically forced surface waves: a systematic survey at large aspect ratio. Physica D, 97, 1996. [65] H. Lamb. Hydrodynamics. Dover Publications, NY, 1945. [66] L.D. Landau and E.M. Lifshitz. Ltd., Oxford, 1959. [67] L.D. Landau and E.M. Lifshitz. Fluid mechanics, 2nd edition. Pergamon Books Ltd., Oxford, 1987. [68] L.D. Landau and E.M. Lifshitz. Statistical Physics, 3rd edition. Butterworth-Heinemann, 1990. [69] S. Luding. Granular materials under vibration: simulations of rotating spheres. Phys. Rev. E, 52:4442, 1995. [70] S. Luding. Global equation of state of two-dimensional hard sphere Fluid mechanics. Pergamon Books systems. Phys. Rev. E, 63:042201, 2001. [71] C. K. K. Lun. A simple kinetic theory for granular flow of rough inelastic spherical particles. J. Appl. Mech., 54:47, 1987. [72] C. K. K. Lun. Kinetic theory for granular flow of dense, slightly inelastic, slightly rough spheres. J. Fluid Mech., 233:539, 1991. 158 [73] C. K. K. Lun, S. B. Savage, D. J. Jeffrey, and N. Chepurniy. Kinetic theories for granular flow: inelastic particles in Couette flow and slightly inelastic particles in a general flowfield. J. Fluid Mech., 140:233, 1984. [74] F. Melo and S. Douady. From solitary waves to static patterns via spatiotemporal intermittancy. Phys. Rev. Lett., 71:3283, 1993. [75] F. Melo, P. Umbanhowar, and H. L. Swinney. Transition to parametric wave patterns in a vertically oscillated granular layer. Phys. Rev. Lett., 72:172, 1994. [76] F. Melo, P. Umbanhowar, and H. L. Swinney. Hexagons, kinks, and disorder in oscillated granular layers. Phys. Rev. Lett., 75:3838, 1995. [77] Thomas H. Metcalf, James B. Knight, and Heinrich M. Jaeger. Standing wave patterns in shallow beds of vibrated granular material. Physica A, 236:202 210, 1997. [78] S. J. Moon. Simulations of Granular Materials: Kinetics and Hydrodynamic Phenomena. Ph.D. thesis, University of Texas at Austin, 2003. [79] S. J. Moon, M. D. Shattuck, C. Bizon, Daniel I. Goldman, J. B. Swift, and H. L. Swinney. Phase bubbles and spatiotemporal chaos in granular patterns. Phys. Rev. E, 65:011301, 2002. [80] Sung Joon Moon, J. B. Swift, and H. L. Swinney. Role of friction in pattern formation in oscillated granular layers. Phys. Rev. E, 69(031301):1, 2004. 159 [81] Nicol s Mujica and Fransico Melo. Solid-liquid transition and hydroa dynamic surface waves in vibrated granular layers. Phys. Rev. Lett., 80:5121, 1998. [82] J. Oh and G. Ahlers. Thermal-noise effect on the transition to RayleighB nard convection. Phys. Rev. Lett., 91:094501, 2003. e [83] J.S. Olafsen and Jeffery S. Urbach. Velocity distributions and density fluctuations in a granular gas. Phys. Rev. E, 60(3):4268, 1999. [84] H. K. Pak, E. Van Doorn, and R. P. Behringer. Effects of gases on granular materials under vertical vibration. Phys. Rev. Lett., 74:4643 4646, 1995. [85] T. Poschel and N. Brilliantov (Editors). Springer-Verlag, Berlin, 2003. [86] A. V. Potapov and C. S. Campbell. Propagation of elastic waves in Granular Gas Dynamics. deep vertically shaken particle beds. Phys. Rev. Lett., 77:4760, 1996. [87] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling. Numerical Recipes in C. Cambridge University Press, Cambridge, 1988. [88] R. Ram irez, D. Risso, and P. Cordero. Thermal convection in fluidized granular systems. Phys. Rev. E, 62:2521, 2000. [89] R. Ram and R. Soto. Temperature inversion in granular fluids under irez gravity. Physica A, 322:73, 2003. 160 [90] D. C. Rapaport. The Art of Molecular Dynamics Simulation. Cambridge University Press, Cambridge, 1997. [91] I. Rehberg, S. Rasenat, M. de la Torre Ju rez, W. Sch pf, F. H rner, a o o and H. R. Brand. Thermally induced hydrodynamic fluctuations below the onset of electroconvection. Phys. Rev. Lett., 67:596, 1991. [92] E. C. Rericha. Shocks in Rapid Granular Flows. Ph.D. thesis, University of Texas at Austin, 2004. [93] E. C. Rericha, C. Bizon, M. D. Shattuck, and H. L. Swinney. Shocks in supersonic sand. Phys. Rev. Lett., 88:014302, 2002. [94] Gerald H. Ristow, Gunther Stra burger, and Ingo Rehberg. Phase diagram and scaling of granular materials under horizontal vibrations. Phys. Rev. Lett., 79(5):833, 1997. [95] S. H rl ck and P. Dimon. Statistics of shock waves in a two-dimensional u granular flow. Phys. Rev. E, 60:671, 1999. [96] S. H rl ck and P. Dimon. Grain dynamics in a two-dimensional granular u flow. Phys. Rev. E, 63:031301, 2001. [97] P. J. Roache. Computational Fluid Dynamics, pp. 232-237. Hermosa Publishers, Albuquerque, 1976. [98] Daniel H. Rothman. Oscillons, spiral waves, and stripes in a model of vibrated sand. Phys. Rev. E, 57(2):R1239 42, 1998. 161 [99] F. Rouyer and N. Menon. Velocity fluctuations in a homogeneous 2d granular gas in steady state. Phys. Rev. Lett, 85(17):3676, 2000. [100] Hidetsugu Sakaguchi and Helmut R. Brand. Standing wave localized squares in pattern-forming nonequilibrium systems. 7(10):1325 30, 1997. [101] S. B. Savage. Streaming motions in a bed of vibrationaly fluidized dry granular material. J. Fluid Mech., 194:457, 1988. [102] M.A. Scherer and G. Ahlers. Deviations from linear theory for fluctuations below the supercritical primary bifurcation to electroconvection. Phys. Rev. Lett., 85:3754, 2000. [103] N. Sela and I. Goldhirsch. Hydrodynamic equations for rapid flows of smooth inelastic spheres - to Burnett order. 1998. [104] Troy Shinbrot. Competition between randomizing impacts and inelastic collisions in granular pattern formation. Nature, 389(6651):574 6, 1997. [105] L.E. Silbert, G.S. Grest, and S.J. Plimption. Boundary effects and selforganuzation in dense granular flows. Phys. Fluids, 14(8):2637 2646, 2002. [106] R. Soto and M.M. Mansour. Granular systems on a vibrating wall: the hydrodynamic boundary condition. Physica A, 327:88 93, 2003. J. Fluid Mech., 361:41, J. de Phys. II, 162 [107] P. Sunthar and V. Kumaran. Characterization of the stationary states of a dilute vibrofluidized granular bed. Phys. Rev. E, 64:041303, 2001. [108] J.B. Swift and P.C. Hohenberg. Hydrodynamic fluctuations at the convective instability. Phys. Rev. A, 15:319, 1977. [109] Harry L. Swinney and E.C. Rericha. Pattern formation and shock waves in granular gases. In F. Mallamace and E. Stanley, editors, The Physics of Complex Systems (New Advances and Perspectives)- The International School of Physics Enrico Fermi. Italian Physical Society, 2004. [110] M-L. Tan and I. Goldhirsch. Rapid granular flows as mesoscopic systems. Phys. Rev. Lett., 81(14):3022, 1998. [111] Sarath G. K. Tennakoon and R. P. Behringer. Vertical and horizontal vibration of granular materials: Coulomb friction and a novel switching state. Phys. Rev. Lett., 81(4):794 797, 1998. [112] Lev Tsimring and Igor Aranson. Localized and cellular patterns in a vibrated granular layer. Phys. Rev. Lett., 79(2):213 216, 1997. [113] P. Umbanhowar and H. L. Swinney. Wavelength scaling and square/stripe and grain mobility transitions in vertically oscillated granular layers. Physica A, 288:344, 2000. 163 [114] Paul Umbanhowar, Francisco Melo, and Harry L. Swinney. Localized excitations in a vertically vibrated granular layer. Nature, 382:793 796, 1996. [115] H. van Beijeren and E.G.D. Cohen. The effects of thermal noise in a Rayleigh-B nard cell near its first convective instability. J. Stat. Phys., e 53:77, 1988. [116] J. van Zon, J. Kreft, D. Goldman, D. Miracle, J.B. Swift, and H. Swinney. Crucial role of sidewalls in velocity distributions in quasi-2d granular gases. Submitted to Phys. Rev. E, 2004. [117] C. T. Veje, Daniel W. Howell, and R. P. Behringer. Kinematics of a twodimensional granular Couette experiment at the transition to shearing. Phys. Rev. E, 59(1):739, 1999. [118] Shankar C. Venkataramani and Edward Ott. Spatio-temporal bifurcation phenomena with temporal period doubling: patterns in vibrated sand. Phys. Rev. Lett., 80:3498, 1998. [119] O.R. Walton. Numerical simulation of inelastic, frictional particleIn M.C. Roco, editor, Particulate Two-Phase particle interactions. Flow, page 884. Butterworth-Heinnemann, 1993. [120] D.R.M. Williams and F. C. MacKintosh. Driven granular media in one dimension: correlations and equation of state. Phys. Rev. E, 54:R9, 1996. 164 [121] M. Wu, G. Ahlers, and D. S. Cannell. Thermally induced fluctuations below the onset of Rayleigh-B nard convection. e 75:1743, 1995. [122] X. Yan, Q. Shi, M. Hou, K. Lu, and C.K. Chan. Effects of air on the segregation of particles in a shaken granular bed. 91:014302, 2003. [123] V.M. Zaitsev and M.I. Shilomis. Soviet Physics JETP, 32:866, 1971. [124] Ya. B. Zel'dovich and Yu. P. Raizer. Physics of shock waves and hightemperature hydrodynamic phenomena. 1966. [125] T. Zhou and L. P. Kadanoff. Inelastic collapse of three particles. Phys. Rev. E, 54(1):623 628, 1996. Academic Press, New York, Phys. Rev. Lett., Phys. Rev. Lett., 165 Vita Jonathan Lee Bougie was born to Michael and Bonnie Bougie in Taranto, Italy on November 16, 1974. When he was two years old, his family moved to Vadnais Heights, MN where he spent the rest of his childhood. Jonathan enrolled in Carleton College in Northfield, MN in 1993, and left four years later with a B.A. in physics. He moved to Austin Texas in 1997 to begin his graduate studies at the University of Texas at Austin. He married Julie Fiorelli on May 29, 2004. Permanent address: 214 A W. 51st St. Austin, Texas 78751 A This dissertation was typeset with L TEX by the author. A LT EX is a document preparation system developed by Leslie Lamport as a special version of Donald Knuth's TEX Program. 166

Find millions of documents here - Study Guides, Homework Solutions, Papers, Exam Answer Keys and more. Course Hero has millions of course related materials that will enable you to learn better, faster and get an A in all your courses.
Below is a small sample set of documents:

manimalajc042.pdf
Path: Texas >> MANIMALAJC >> 042 Fall, 2009

doppmanngw026.pdf
Path: Texas >> DOPPMANNGW >> 026 Fall, 2009
guajardoma026.pdf
Path: Texas >> GUAJARDOMA >> 026 Fall, 2009
martinezvm029.pdf
Path: Texas >> MARTINEZVM >> 029 Fall, 2009
brownsonab029.pdf
Path: Texas >> BROWNSONAB >> 029 Fall, 2009
lambertg36961.pdf
Path: Texas >> LAMBERTG >> 36961 Fall, 2009
aljuaiedma042.pdf
Path: Texas >> ALJUAIEDMA >> 042 Fall, 2009
decastropj029.pdf
Path: Texas >> DECASTROPJ >> 029 Fall, 2009
cathrodl77285.pdf
Path: Texas >> CATHRODL >> 77285 Fall, 2009
mcglohenmk042.pdf
Path: Texas >> MCGLOHENMK >> 042 Fall, 2009
lansdellcp029.pdf
Path: Texas >> LANSDELLCP >> 029 Fall, 2009
stuberja80926.pdf
Path: Texas >> STUBERJA >> 80926 Fall, 2009
canterar35023.pdf
Path: Texas >> CANTERAR >> 35023 Fall, 2009
chatellemb042.pdf
Path: Texas >> CHATELLEMB >> 042 Fall, 2009
shackmanlc042.pdf
Path: Texas >> SHACKMANLC >> 042 Fall, 2009
complexity.txt
Path: CSU San Bernardino >> CS >> 330 Fall, 2009
okazakit51686.pdf
Path: Texas >> OKAZAKIT >> 51686 Fall, 2009
bamfordw82161.pdf
Path: Texas >> BAMFORDW >> 82161 Fall, 2009
russellr74662.pdf
Path: Texas >> RUSSELLR >> 74662 Fall, 2009
lab9.pdf
Path: CSU San Bernardino >> CS >> 201 Fall, 2009
mukadama15106.pdf
Path: Texas >> MUKADAMA >> 15106 Fall, 2009
kellerkm71167.pdf
Path: Texas >> KELLERKM >> 71167 Fall, 2009
oxfordwt32223.pdf
Path: Texas >> OXFORDWT >> 32223 Fall, 2009
bennettl81291.pdf
Path: Texas >> BENNETTL >> 81291 Fall, 2009
engelas504835.pdf
Path: Texas >> ENGELAS >> 504835 Fall, 2009
curranma71134.pdf
Path: Texas >> CURRANMA >> 71134 Fall, 2009
stanleyk74304.pdf
Path: Texas >> STANLEYK >> 74304 Fall, 2009
protsenkode026.pdf
Path: Texas >> PROTSENKOD >> 026 Fall, 2009
Chapter07.outline.pdf
Path: Concordia NE >> PHYS >> 110 Fall, 2009
rutherfordg022.pdf
Path: Texas >> RUTHERFORD >> 022 Fall, 2009
auerbachs13838.pdf
Path: Texas >> AUERBACHS >> 13838 Fall, 2009
dechapanyaw029.pdf
Path: Texas >> DECHAPANYA >> 029 Fall, 2009
shoemakerdb042.pdf
Path: Texas >> SHOEMAKERD >> 042 Fall, 2009
johnsonam71217.pdf
Path: Texas >> JOHNSONAM >> 71217 Fall, 2009
sampselld77810.pdf
Path: Texas >> SAMPSELLD >> 77810 Fall, 2009
complex.txt
Path: CSU San Bernardino >> CS >> 330 Fall, 2009
cadenheadjk046.pdf
Path: Texas >> CADENHEADJ >> 046 Fall, 2009
benjaminsmr042.pdf
Path: Texas >> BENJAMINSM >> 042 Fall, 2009
simpsonal13317.pdf
Path: Texas >> SIMPSONAL >> 13317 Fall, 2009
hamiltont84490.pdf
Path: Texas >> HAMILTONT >> 84490 Fall, 2009
kotrlaka518287.pdf
Path: Texas >> KOTRLAKA >> 518287 Fall, 2009
harrisont86130.pdf
Path: Texas >> HARRISONT >> 86130 Fall, 2009
brandonjc99738.pdf
Path: Texas >> BRANDONJC >> 99738 Fall, 2009
MATH107A46024536.doc
Path: MD University College >> ASIA >> 2092 Fall, 2009
crawforda65881.pdf
Path: Texas >> CRAWFORDA >> 65881 Fall, 2009
achacosom07761.pdf
Path: Texas >> ACHACOSOM >> 07761 Fall, 2009
jarroldwl86380.pdf
Path: Texas >> JARROLDWL >> 86380 Fall, 2009
sharyginany026.pdf
Path: Texas >> SHARYGINAN >> 026 Fall, 2009
goncalvesac026.pdf
Path: Texas >> GONCALVESA >> 026 Fall, 2009
zieglerkj47418.pdf
Path: Texas >> ZIEGLERKJ >> 47418 Fall, 2009
burtnerjc90760.pdf
Path: Texas >> BURTNERJC >> 90760 Fall, 2009
alvarezla07232.pdf
Path: Texas >> ALVAREZLA >> 07232 Fall, 2009
MATH012A46124534.doc
Path: MD University College >> ASIA >> 2092 Fall, 2009
bonningew86532.pdf
Path: Texas >> BONNINGEW >> 86532 Fall, 2009
CMIS141AA44024445.doc
Path: MD University College >> ASIA >> 2092 Fall, 2009
CMIS102AA42086692.doc
Path: MD University College >> ASIA >> 2088 Fall, 2009
STAT200A42186896.doc
Path: MD University College >> ASIA >> 2088 Fall, 2009
kulkarnis86095.pdf
Path: Texas >> KULKARNIS >> 86095 Fall, 2009
chapmanbg60287.pdf
Path: Texas >> CHAPMANBG >> 60287 Fall, 2009
slattonkc78713.pdf
Path: Texas >> SLATTONKC >> 78713 Fall, 2009
michalskylo026.pdf
Path: Texas >> MICHALSKYL >> 026 Fall, 2009
batemanmt33508.pdf
Path: Texas >> BATEMANMT >> 33508 Fall, 2009
lodowskid97061.pdf
Path: Texas >> LODOWSKID >> 97061 Fall, 2009
raichlend29983.pdf
Path: Texas >> RAICHLEND >> 29983 Fall, 2009
perkinsjd44616.pdf
Path: Texas >> PERKINSJD >> 44616 Fall, 2009
mehdiabadinj026.pdf
Path: Texas >> MEHDIABADI >> 026 Fall, 2009
borisovasa86653.pdf
Path: Texas >> BORISOVASA >> 86653 Fall, 2009
Abuhakema504399.pdf
Path: Texas >> ABUHAKEMA >> 504399 Fall, 2009
hw03_solution.doc
Path: Penn State >> ME >> 581 Fall, 2009
oestreichj19588.pdf
Path: Texas >> OESTREICHJ >> 19588 Fall, 2009