{'fraction_non_alphanumeric': 0.031085262434104973, 'fraction_numerical': 0.012548705019482008, 'mean_word_length': 4.505520504731861, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Continuum simulations of shocks and patterns in vertically oscillated granular layers

We study interactions between shocks and standing-wave patterns in vertically oscillated layers of granular media using three-dimensional, time-dependent numerical solutions of continuum equations to Navier-Stokes order. We simulate a layer of grains atop a plate that oscillates sinusoidally in the direction of gravity. Standing waves form stripe patterns when the accelerational amplitude of the plate's oscillation exceeds a critical value. Shocks also form with each collision between the layer and the plate; we show that pressure gradients formed by these shocks cause the flow to reverse direction within the layer. This reversal leads to an oscillatory state of the pattern that is subharmonic with respect to the plate's oscillation. Finally, we study the relationship between shocks and patterns in layers oscillated at various frequencies and show that the pattern wavelength increases monotonically as the shock strength increases.

In this paper, we model granular media by numerically solving a set of time-dependent continuum equations for the rapid flow of dissipative particles in three-dimensions. We use these simulations to investigate the relationship between shocks and standing-wave patterns in vertically oscillated granular layers.

B. Granular hydrodynamics
A successful theory of granular hydrodynamics would allow scientists and engineers to apply the powerful methods of fluid dynamics to granular flow. Despite experimental [7][8][9] and computational [8,[10][11][12][13] evidence demonstrating the potential utility of hydrodynamics models for grains, a general set of hydrodynamic governing equations is not yet recognized for granular media [14][15][16][17].
Several proposed rapid granular flow models use equations of motion for continuum fields: number density n, velocity u, and granular temperature T ( 3 2 T is the average kinetic energy due to random particle motion) [18][19][20]. In one approach, particle interactions are modeled with binary, inelastic hard-sphere collision operators in kinetic theory to derive continuum equations to Euler [21], Navier-Stokes [22], and Burnett [23] order. We use three-dimensional (3D) simulations of continuum equations to Navier-Stokes order to investigate shocks and standing-wave patterns in shaken granular layers.

C. Standing wave patterns in oscillated granular layers
A granular layer of depth H atop a plate that is oscillated sinusoidally in the direction of gravity with frequency f and amplitude A will leave the plate at some time during the oscillation cycle if the maximum acceleration of the plate a max = A (2πf ) 2 is greater than the acceleration of gravity g. In other words, the layer leaves the plate if the nondimensional accelerational amplitude Γ = a max /g exceeds unity.
After leaving the plate, the layer dilates above the plate, and then compresses when it collides with the plate later in the cycle. When the dimensionless accelerational amplitude Γ is larger than a critical value Γ C , standing-wave patterns spontaneously form in the layer. These patterns are subharmonic with respect to the plate, repeating every 2/f [6]. Depending on the nondimensional accelerational amplitude Γ and dimensionless frequency f * = f H/g, various subharmonic standing waves have been found, including stripe, square, and hexagonal standing-wave patterns [6].

D. Shocks in oscillated granular layers
If the Mach number M a (the ratio of the local mean fluid speed to the local speed of sound) is greater than unity at the point where a fluid encounters an obstacle, a compression wave front is formed near the object and steepens to form a shock. A distinguishing feature of granular materials is that granular flows reach supersonic speeds under common laboratory conditions [3,9,18]. Therefore, although shocks are formed in ordinary gases only under extreme conditions, shock formation is commonplace in granular media.
Experiments [3] and simulations [11,[24][25][26] demonstrate that shock waves form in shaken granular layers as the layer contacts the plate. Previous investigations have generally either focused on shock propagation [3,24,25] or pattern formation [6,12,13] or considered them as separate phenomena coexisting in shaken layers [11]. In this paper, we use continuum simulations to investigate the interactions between shocks and the standing-wave patterns formed in this system.

E. Model system
We simulate a layer of grains on an impenetrable plate which oscillates sinusoidally in the direction of gravity. The layer depth at rest is approximately H = 5.4σ, where the grains are modeled as identical, frictionless spheres with diameter σ, mass M , and coefficient of restitution e = 0.7. In this paper, we study patterns and shocks as a function of nondimensional frequency f * , while the dimensionless accelerational amplitude Γ = 2.20 is held constant. Previous simulations [12,13] have shown that this value exceeds the critical accelerational amplitude Γ C indicating that we should expect to see standing-wave patterns for a variety of frequencies.
Previous experiments [27] and molecular dynamics (MD) simulations [28] have shown that friction between grains plays a role in these patterns. Experimentally, adding graphite to reduce friction decreased Γ C and prevented the formation of stable square or hexagonal patterns found for certain ranges of f * and Γ in experiments without graphite [27]. Similarly, MD simulations with friction between particles have quantitatively reproduced stripe, square, and hexagonal subharmonic standing waves seen experimentally [29], but MD simulations without friction yield only stable stripe patterns and display a lower Γ C [28]. In this study, we investigate stripe patterns in continuum simulations of frictionless particles. To investigate other patterns such as squares or hexagons, simulations would have to include friction between particles.
Experiments [30] and simulations [12,13,31] indicate that fluctuations due to individual grain movement play a larger role in granular media than do thermal fluctuations in ordinary fluids. Fluctuating hydrodynamics (FHD) theory models these fluctuations by adding noise terms to the Navier-Stokes equations [32][33][34]. In previous simulations with FHD noise terms, the critical value of pattern onset Γ C was consistent with molecular dynamics simulations, while continuum simulations without these fluctuations exhibited pattern onset at Γ C approximately 10% lower than that found in MD simulations [12,13]. Above Γ C , however, simulations both with and without these fluctuations exhibited standing-wave patterns with wavelengths consistent with a dispersion relation [35] found experimentally for a range of shaking frequencies. In this paper, we do not include FHD noise terms; we investigate patterns for accelerational amplitude greater than Γ C determined from simulations with and without FHD.
We use continuum simulations to investigate the dynamics of this system including pattern formation and shock propagation. Section II describes the methods we use to simulate and analyze oscillated layers. Sec. III examines the dynamics of shocks and standing-wave patterns formed in this system at a fixed dimensionless oscillation frequency f * = 0.25, and Sec. IV examines how these patterns and shocks change when this frequency is varied. We present our conclusions in Sec. V.

A. Continuum equations
We use a continuum simulation previously used to model shock waves in a granular shaker [26]. Our simulation numerically integrates continuum equations of Navier-Stokes order proposed by Jenkins and Richman [22] for a dense gas composed of frictionless (smooth), inelastic hard spheres.
This model yields hydrodynamic equations for number density n (or equivalently, volume fraction ν = π 6 nσ 3 ), velocity u, and granular temperature T : where the components of the symmetrized velocity gradient tensor E are given by: . The components of the stress tensor P are given by the constitutive relation and the heat flux is calculated from Fourier's law: To calculate the pressure, we use the equation of state and radial distribution function at contact proposed by Goldshtein et al. [36] to include both dense gas and inelastic effects: where and the radial distribution function at contact, g 0 , is where ν max = 0.65 is the 3D random close-packed volume fraction. These equations differ from those for a compressible, dense gas of elastic particles by the energy loss term which arises from the inelasticity of collisions between particles. The bulk viscosity is given by the shear viscosity by and the thermal conductivity by (12) Other hydrodynamic models include modifications such as more accurate expressions for kinetic coefficients which incorporate high density corrections [37,38] and equations of state which allow for volume fractions greater than the random-close-packed limit used here [39]. However, the equations shown here have previously been used to separately examine shocks [26] and patterns [12,13] in shaken layers and have demonstrated quantitative agreement with experiments and molecular dynamics simulations. Therefore, we use these equations in this study to investigate the relationship between shocks and pattern formation; in principle, the equations could be modified to implement other forms of the constitutive relations.

B. Simulation method
We integrate these hydrodynamic equations to find number density, momentum, and granular temperature, using a second order finite difference scheme on a uniform grid in 3D with first order adaptive time stepping [26]. In these simulations, the granular fluid is contained between two impenetrable horizontal plates at the top and bottom of the container, where the lower plate oscillates sinusoidally between height z = 0 and z = 2A and the ceiling is located at a height L z above the lower plate. Periodic boundary conditions are used in the horizontal directions x and y to eliminate sidewall effects. Simulations were conducted in a box of size L x = 168σ, L y = 10σ, and L z = 160σ. This orientation causes stripes to form parallel to the y-axis. The numerical methods, boundary conditions at the top and bottom plate, and grid spacing are the same as used in previous studies of shocks [26] and patterns [12,13].

A. Standing wave patterns
Experimental investigations of shaken granular layers have shown that above a critical acceleration of the plate Γ C , standing wave patterns form spontaneously. These patterns oscillate subharmonically, repeating every 2/f , so that the location of a peak of the pattern becomes a valley after one cycle of the plate, and vice versa [6].
Continuum simulations produce standing wave patterns for Γ = 2.2 and f * = 0.25 (Fig. 1). Beginning with a flat layer above the plate with small amplitude random fluctuations, the simulation ran for 250 cycles of the plate until the layer reached a periodic state. Snapshots from various times during the next two cycles of the plate are shown in Fig. 1. 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. We examine the wave patterns at various times in the cycle f t for two cycles of the plate. At f t = 0 the container is at its minimum height. The bulk of the layer is dilute and four clear peaks and four valleys are evident in Fig. 1 (note that one of the peaks wraps around the edges of the cell since the boundary conditions are periodic). The material in each peak is falling towards the plate, although the bottom of the layer has begun to contact the plate and is beginning to form a pile on the plate. The thin layer colliding with the plate is characterized by a higher volume fraction ν than the material in the peaks. At f t = 0.25, the plate has begun move upward as the layer falls to meet it. At this time, most of the layer is piling on the plate, forming a high volume fraction region near the plate. Remnants of the peaks and valleys from the previous cycle are still visible, although the layer is flattening on the plate, with the height difference between peaks and valleys much smaller than at f t = 0 ( Fig. 1).
3. f t = 0.50: The layer begins to leave the plate.
As the container reaches its maximum downward acceleration at the top of its oscillation, the layer begins to leave the plate. Although the bottom of the layer is still in contact with the plate at this time, the layer Volume Fraction is expanding, as evidenced by the fact that it is visibly more dilute in this picture than in the previous snapshot ( Fig. 1). As the layer leaves the plate, new peaks and valleys develop, with peaks forming in the valleys from the previous cycle, and vice versa. 4. f t = 0.75: The layer is off the plate.
At f t = 0.75, the plate is moving downward leaving a gap between the layer and the plate. The mass of the layer has almost entirely left the plate, and the layer has expanded. The new peaks and valleys have become quite distinct, with a large portion of the layer in the peaks and very little material in the valleys connecting them (Fig. 1). 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 are similar to those at f t = 0, except the peaks and valleys have reversed location. Aside from the peaks and valleys reversing, the next cycle exhibits the same features as the previous cycle at all times during the cycle. By f t = 1.75 the peaks and valleys are back to their original location, demonstrating the subharmonic nature of these patterns.

B. Shocks
Previous experiments [3] and simulations [11,[24][25][26] indicate that shocks are formed in vertically oscillated granular layers during each collision of the layer with the plate. 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. Therefore, we calculate the speed of sound using a relation derived from the equation of state Eq. (6) [40]: where χ = 1 + 2(1 + e)G(ν). The Mach number M a is calculated as the local flow speed divided by the local speed of sound M a = |u| /c. From Fig. 1, we see that at f t = 1.00, a portion of the layer has begun to compress on the plate, while some of the material is still falling towards the plate. This is therefore a time at which we would expect to see shock formation. To examine the shock in detail, we plot onedimensional profiles of volume fraction ν and Mach number M a (here ν is scaled by a factor of 10 to fit on the same scale as M a) as functions of height z from the plate in Fig. 2.
Horizontal Location x = 32σ At f t = 1.00, the first peak of the pattern occurs at the horizontal location 16σ x 20σ (cf Fig. 1). We choose two horizontal locations to examine; the location y = 5σ, and x = 16σ is shown in the left column of Fig. 2, corresponding to the first peak of the pattern; the location y = 5σ, and x = 32σ is shown in the right column, corresponding to a location away from the peak.
At the horizontal location y = 5σ, x = 16σ corresponding to a peak of the wave pattern, the volume fraction approaches the close-packed value ν ≈ 0.6 near the plate at f t = 1.00 (Fig. 2a). As height increases, ν smoothly decreases for 0 ≤ z 7σ. The Mach number with respect to the plate M a is low throughout this region, as this part of the layer is compressed on the plate and thus matches the velocity of the plate at height z = 0. However, at z ≈ 7σ, there is a discontinuity in the derivative of ν, and a sharp increase in Mach number such that M a ≈ 3 for z 8σ. This sharp increase in Mach number from subsonic to supersonic, corresponding to a discontinuity in the derivative of volume fraction at the same height, indicates that this is the location of a strong shock front.
While the above description applies to x = 16σ (the horizontal location of the first peak), a shock also forms in the valleys. This can be seen in Fig. 2b, which shows 10ν and M a at the same time, but at horizontal location x = 32σ. Here we again see a high density, subsonic region near the plate, with a shock separating this region from a lower density, supersonic undisturbed region that is still falling towards the plate. However, the volume fraction near the plate is smaller than it was in Fig. 2a, and the shock profile is quite different. Also, the shock front at x = 32σ is closer to the plate than it is at x = 16σ.
By f t = 1.10, the shock has moved away from the plate and it continues to develop as it moves through different parts of the layer (Fig. 2c). The shock front in the valleys is also developing and moving away from the plate (see Fig. 2d), although the location and profile of the shock is different at these different horizontal locations.
Thus, shocks form with each collision of the layer with the plate, and the shock profile changes and develops as the shock moves through the layer. Shocks form at each horizontal location in the cell; however the shock front is non-uniform horizontally due to the differences in layer depth between the peaks and the valleys of the pattern.

C. Interaction between shocks and patterns
There is therefore a relationship between the shocks and the standing-wave patterns formed in the layer. If the layer were perfectly flat, the shock front would be uniform horizontally. However, as there is horizontal variation in layer depth, a non-uniform shock front forms.
While the system is driven by vertical motion of the plate, the patterns are characterized by horizontal variation between peaks and valleys. Periodic horizontal motion is required to produce the subharmonic oscillation. This horizontal flow can be thought of as a "sloshing" motion, in which material from the layer flows out of some regions and into others, creating peaks and valleys. With the next collision of the layer with the plate, the flow direction reverses, creating peaks where there were valleys, and vice versa. In this section, we show that strong pressure gradients created by the shock drive the sloshing motion of the layer.
To study this interaction, we investigate the layer properties near the time of shock formation. According to Fig. 1, most of the layer is off the plate at f t = 0.75, while the layer is compressed on the plate and the peaks and valleys have flattened significantly by f t = 1.25. In Fig. 3 we show snapshots of the the volume fraction ν (left column) and the dimensionless pressure p σ 2

Mg
(right column) at five different times 0.8 ≤ f t ≤ 1.2. We display side views of the layer at y = 5σ in (Fig 3), showing only the horizontal range 0 ≤ x ≤ 84σ and the vertical range 0 ≤ z ≤ 25σ to better view the interaction between the layer and the plate. This range shows two wavelengths of the standing wave pattern and corresponds to the bottom left quadrant of the range plotted in Fig. 1.
In order to examine the sloshing motion of the layer, we also show the dimensionless x-component of the pressure gradient ∂p ∂x σ 3 Mg (left column) and the dimensionless horizontal velocity u x √ gσ (right column) at the these same times and locations in Fig. 4. We now proceed to examine the horizontal variation in flow at these five times as the layer collides with the plate. At f t = 0.8, the layer is mostly off the plate, and peaks and valleys are clearly visible in the left column of Fig. 3. In the right column of Fig. 3, we show the dimensionless pressure p σ 2 Mg . At this time, the layer is more dilute than it will be at later times in the cycle, and the pressure is uniformly low throughout the layer. There are no strong pressure gradients visible, as can be seen in the left column of Fig. 4, where we show the dimensionless x-component of the pressure gradient ∂p ∂x σ 3 Mg . In the right column of Fig. 4, we can see that u x > 0 (the flow is moving towards the right) for 0 x 16σ; there is a small range of nearly zero horizontal flow rate near 16σ x 20σ; then u x > 0 (flow is moving towards the left) for 16σ x 34σ. Since the peak of the pattern is located at 16σ x 20σ at this time (cf Fig. 3), grains are flowing towards the peak from both the left and the right. Therefore, the flow is still removing material from the valleys and moving it towards the peak; the difference between peaks and valleys is growing at this time.
2. f t ≈ 0.9 : The layer begins to contact the plate.
At f t = 0.9, the bottom of the layer has just begun to contact the plate at the horizontal location of the two peaks, but the material between the peaks has not yet contacted the plate (Fig. 3). The layer looks similar to its appearance at f t = 0.8, except that in a very small region near the plate, the pressure has begun to slightly increase where the layer is starting to contact the plate.
This pressure increase at the two points of contact yields a horizontal pressure gradient with maximum pressure near the peak location. For instance, the reddish color found near the plate for 12σ x 18σ in the left column of Fig. 4 at this time indicates that the pressure is increasing from left to right; the bluish color for 18σ x 26σ in the same picture indicates that the pressure decreases from left to right.
The picture in the right column for this time shows that the flow is still moving towards the peak through the bulk of the layer at this time. However, the pressure gradient near the plate tends to drive flow in the opposite direction: from high pressure (near the peak) into low pressure (near the valley). Therefore, the horizontal flow velocity near the plate has begun to slow (leaving a white space near the plate where the pressure gradient is strongest) and is even beginning to reverse direction (note the light blue near the plate at the location of the second peak).
By f t = 1.0, more of the material is colliding with the plate, and the pressure has noticeably increased at the two main contact points between the layer and the plate (Fig. 3). This increase in layer density and pressure near the plate marks the formation of a shock (cf Fig. 2).
As discussed in Sec. III B, these shocks are not uniform horizontally; larger layer depth near the peaks leads to larger pressure at these locations. Thus, shock formation leads to pressure maxima near the plate at the horizontal location of the pattern peaks with strong pressure gradients as seen in the left column of Fig. 4.
These pressure gradients drive the flow in the direction opposite the pressure gradient; thus the direction of flow near the plate at this time is clearly reversing (see the right column of Fig. 4). For much of the layer, the flow is still away from the valleys and towards the peaks. Near the plate in the region of strong pressure gradients, In the left column, 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 gray as volume fraction increases. The plate is represented as a thick, horizontal black line. In the right column, low dimensionless pressure is white, high dimensionless pressure is black, and pressure increases from white to black through shades of gray.
however, the flow has reversed itself and is now flowing away from the peaks and towards the valleys. This separation of the flow into two distinct domains represents the formation of a shock near the plate. This shock then travels upward through the layer, away from the plate. By f t = 1.1, this shock is well-developed and has propagated through most of the layer, separating a high-density and pressure region near the plate from a low-density and pressure region above the shock front (Fig. 3).
The pressure gradient is quite strong near the plate ( Fig. 4) and most of the layer is behind the shock. The flow behind the shock has reversed direction, while the material ahead of the shock has not yet changed direction. There is still some material in the peaks that is still falling towards the plate ahead of the shock (cf Fig. 3), but its horizontal velocity is small, and the pressure is relatively low in this region (see Fig. 4). By f t = 1.2, nearly the entire layer is behind the shock and the layer has noticeably flattened (Fig. 3), although peaks and valleys are still visible. As the shock propagates into the low-density region above the layer, the dimensionless horizontal velocity ux √ gσ (right column) at a slice y = 5σ parallel to the x-z plane at various times f t.
Although the full cell has a length 168σ in the x-direction and a height 160σ in the z-direction, this figure only displays the horizontal range 0 ≤ x ≤ 84σ and the vertical range 0 ≤ z ≤ 25σ to show a closer view of the collision of the shock with the plate. In the left column, pressure decreasing left to right ( ∂p ∂x < 0) is shown in shades of blue, ∂p ∂x = 0 is white, and pressure increasing left to right ( ∂p ∂x > 0) is shown in shades of red. In the right column, flow towards the left (ux < 0) is shown in shades of blue, zero horizontal velocity (ν = 0) is white, and flow towards the right (ux > 0) is shown in shades of red. In both columns, the plate is shown as a thick, horizontal black line.
layer cools behind the shock due to inelastic collisions. This causes the pressure to decrease slightly at f t = 1.2 as compared to f t = 1.1.
The pressure gradient is still driving the flow away from the peaks, although this gradient is not as strong as in at f t = 1.1 (Fig. 4). At this time, nearly the entire layer has reversed direction and is now flowing away from the peaks and towards the valleys.
As this flow continues, the layer becomes flat, and then develops peaks where there were previously valleys, and vice versa. Thus, the horizontal pressure gradients created by the nonuniform shock front drive the flow to reverse itself. It is this reversal that creates the sloshing motion visible in the subharmonic oscillation.

IV. FREQUENCY DEPENDENCE OF SHOCKS AND STANDING WAVES
In Sec. III, we studied the interaction between shocks and patterns in layers oscillated with a particular nondimensional frequency f * = 0.25. In this section, we investigate frequency dependence by holding dimensionless accelerational amplitude Γ = 2.2 constant, while varying dimensionless frequency f * . Experiments have shown that wavelength λ depends on the frequency of oscillation [35,41]. For a range of layer depths and oscillation frequencies, experimental data for frictional particles near pattern onset were fit by the function λ * = 1.0 + 1.1f * −1.32±0.03 , where λ * = λ/H [35]. We examine the correlation between changes in shock properties with changes in pattern wavelength throughout this frequency range. For each frequency, we start from an artificial flat layer, then simulate 250 cycles of the plate to allow the layer to reach an oscillatory state. Data is taken from the next six cycles of the plate.

A. Maximum Mach number varies inversely with driving frequency
Shock properties depend on the Mach number of the layer with respect to the plate during collision. In Fig. 5, we plot the maximum Mach number M a max of the layer with respect to the plate as a function of dimensionless driving frequency f * . To calculate M a max , we find the Mach number M a (x, t) at each location x in the cell during an oscillation cycle. To ensure that we are looking at the Mach number of material in the layer, rather than in the low-density region above the layer, at each time t, we find the highest Mach number corresponding to at least 1% of the total mass of the layer. We then define M a max as the maximum Mach number found at any time during a cycle. Note that M a max monotonically decreases as f * increases (Fig. 5). For fixed Γ = 2.20, lower oscillation frequency corresponds to a higher maximum plate velocity and higher oscillation amplitude. Layers with lower f * will generally reach a higher maximum height when they leave the plate, and will have a larger speed relative to the plate when they collide later in the cycle. Although the relative Mach number between the layer and the plate depends on the speed of sound in the layer as well as the time during the cycle at which the collision takes place, higher frequencies correspond to lower M a at collision throughout this range.

B. Higher Mach number at collision produces stronger shocks
Pressure decreases rapidly when moving from the region behind a shock, across the shock front, to the undisturbed region ahead of the shock. We examine this pressure change by calculating the magnitude of the zcomponent of the pressure gradient ∂p ∂z (x, t) at each time t and location in the cell x. We then find ∂p ∂z max (t) as the highest value corresponding to at least 1% of the total mass of the layer. Finally, we calculate the time average of these layers ∂p ∂z max , nondimensionalize by a factor of σ 3 Mg , and plot this as a function of M a max in Fig. 6a. As the maximum Mach number M a max increases, the time-averaged maximum vertical pressure gradient increases monotonically. In other words, the higher the relative Mach number of the layer with respect to the plate, the sharper the pressure drop across the shock. Note that hydrostatic pressure would cause pressure to vary with depth even in static layers, consistent with the non-zero intercept of Fig. 6a.
As discussed in Sec. III, although the shock is produced by collision with a vertically oscillating plate, horizontal pressure variation develops as a result of the peaks and valleys in the layer. Shocks with stronger vertical pressure gradients correspond to stronger horizontal pressure gradients as well (Fig. 6b). A flat layer with no patterns would be expected to have pressure variation with depth but not with horizontal location. This is consistent with the apparent non-zero intercept in Fig. 6b.

C. Horizontal pressure gradients drive horizontal velocity
As strong pressure gradients develop in the layer, there will be a tendency for the material to flow from high pressure to low pressure; i.e. the direction of flow velocity will be opposite the pressure gradient (cf Fig. 4). We calculate the average flow speed in the x-direction at all times and all locations in the cycle |u x | . Figure 7 shows that the nondimensionalized average flow speed in the x−direction |u x | / √ gσ increases monotonically as the maximum dimensionless x-component of the pressure gradient found anywhere in the cell averaged over all times in the cycle ∂p ∂x max σ 3 Mg increases. Thus, the stronger the maximum horizontal pressure gradients produced across the shock are, the faster the average horizontal flow speed becomes. Faster flow in the x-direction corresponds to more horizontal motion of the particles throughout an oscillation cycle. Therefore, the wavelength (and therefore the distance between peaks) increases monotonically as the average horizontal flow speed increases. (see Fig. 8).

E. Dispersion relation
Summarizing the results from this section thus far, we find that higher frequencies of oscillation yield smaller Mach numbers between the layer and the plate during collision (Fig. 5). Larger values of the Mach number produce stronger shocks, as reflected in larger vertical and horizontal changes in pressure across the shock front (Fig. 6). Larger pressure gradients produce faster horizontal flow (Fig. 7) which in turn corresponds to larger horizontal distances between peaks of the pattern, and thus a larger wavelength (Fig. 8).
Given this chain of reasoning, it follows that as the Mach number of the layer with respect to the plate increases, wavelength of the pattern increases monotonically throughout the range. As Mach number increases with decreasing shaking frequency, wavelength should decrease monotonically as the oscillation frequency increases. Figure 9 shows that this is indeed the case throughout the range of frequencies we have simulated. Additionally, the wavelengths found in our simulation throughout this range are consistent with the dispersion relation previously found by fit to experimental data [35]. : Dimensionless wavelength λ * as a function of dimensionless frequency f * . Data from our simulations are shown as points, with error bars calculated exclusively from discretization due to periodic boundary conditions in a finite-size box. The dominant wavelength found in our simulations fits quite well the dispersion relation λ * = 1.0 + 1.1f * −1.32±0.03 (solid line) found as a fit to previous experimental data [35].

V. CONCLUSIONS
We have simulated vertically oscillated layers of granular media by numerically solving a proposed set of granular hydrodynamic equations to Navier-Stokes order. We have shown that these continuum simulations can describe important aspects of shocks and patterns in granular materials. In our simulations, layers spontaneously form subharmonic standing-wave patterns in which the wavelength λ depends on the frequency f of the oscillation of the plate. The wavelengths of the patterns formed in our simulation agree with the dispersion relations found in previous experiments.
As these standing waves oscillate subharmonically, shocks are formed with each collision of the layer with the plate. We have analyzed these shocks and found that variation in the layer depth produces a shock front that is not uniform horizontally. This horizontal variation in the shock front produces horizontal components of flow velocity, causing grains to move from high pressure to low pressure. Therefore, despite the fact that the plate oscillates vertically, the shocks produced during collision between the layer and the plate drive the horizontal sloshing motion that characterizes the standing-wave patterns.
We have simulated shocks in layers oscillated at various frequencies, and have shown that the strength of the shock varies depending on the oscillation frequency. We have also established that there is a relationship between pattern wavelength and the Mach number of the layer with respect to the plate at the time of collision. Specifically, the horizontal and vertical components of the pressure gradient both increase monotonically as the Mach number increases. The horizontal flow speed increases monotonically with the strength of these pressure gradients, leading to longer wavelengths.
Therefore, throughout the range we studied, the wavelength of the pattern increases monotonically with the maximum M a of the layer with respect to the plate. This analysis indicates that shocks play a significant role in the dynamics of the standing-wave patterns formed in these oscillating layers.