Soft-sphere simulations of a planar shock interaction with a granular bed

Here we consider the problem of shock propagation through a layer of spherical particles. A point particle force model is used to capture the shock-induced aerodynamic force acting upon the particles. The discrete element method (DEM) code liggghts is used to implement the shock-induced force as well as to capture the collisional forces within the system. A volume-fraction-dependent drag correction is applied using Voronoi tessellation to calculate the volume of fluid around each individual particle. A statistically stationary frame is chosen so that spatial and temporal averaging can be performed to calculate ensemble-averaged macroscopic quantities, such as the granular temperature. A parametric study is carried out by varying the coefficient of restitution for three sets of multiphase shock conditions. A self-similar profile is obtained for the granular temperature that is dependent on the coefficient of restitution. A traveling wave structure is observed in the particle concentration downstream of the shock and this instability arises from the volume-fraction-dependent drag force. The intensity of the traveling wave increases significantly as inelastic collisions are introduced. Downstream of the shock, the variance in Voronoi volume fraction is shown to have a strong dependence upon the coefficient of restitution, indicating clustering of particles induced by collisional dissipation. Statistics of the Voronoi volume are computed upstream and downstream of the shock and compared to theoretical results for randomly distributed hard spheres.

Figure 1

Representative simulation setup in the shock-attached frame of reference: Randomly packed particles enter on the left-hand side at the laboratory shock speed. At $x=0,$ they encounter a planar shock and experience a shock modulated force, causing collisions.

Figure 2

(Left) A sample random distribution of particles. (Right) The particle distribution encapsulated by a Voronoi tessellation. The edges of the Voronoi cells are shown in red.

Figure 3

The nondimensional drag force for a shock passing over an isolated particle. The point particle force model discussed presently is compared to experimental and direct numerical simulation data for a frozen particle. The point particle model captures the net momentum impulse well.

Figure 4

A pseudocolor plot of the Voronoi volume fraction for case 1. A green line at $x=0$ shows the location of the planar shock. Three coefficients of restitution are plotted: $e=1$ (top), $e=0.7$ (middle), and $e=0.4$ (bottom).

Figure 5

A pseudocolor plot of the Voronoi volume fraction for case 2. A green line at $x=0$ shows the location of the planar shock. Three coefficients of restitution are plotted: $e=1$ (top), $e=0.7$ (middle), and $e=0.4$ (bottom).

Figure 6

A pseudocolor plot of the Voronoi volume fraction for case 3. A green line at $x=0$ shows the location of the planar shock. Three coefficients of restitution are plotted: $e=1$ (top), $e=0.7$ (middle), and $e=0.4$ (bottom).

Figure 7

Ensemble-averaged $\varphi$ (a) and $x$ velocity (b) compared with the fixed-time, planar averages for case 1. The variation of $\overline{\varphi }$ is small compared to the variation seen in cases with lower ${\varphi }_e$.

Figure 8

Ensemble-averaged $\varphi$ (a) and $x$ velocity (b) compared with the fixed-time, planar averages for case 2. Variations in $\overline{\varphi }$ of almost $10\%$ are seen both above and below the ensemble-averaged values.

Figure 9

Ensemble-averaged $\varphi$ (a) and $x$ velocity (b) compared with the fixed-time, planar averages for case 3. Similar to case 2, variations in $\overline{\varphi }$ of almost $10\%$ are seen both above and below the ensemble-averaged values.

Figure 10

Histogram Voronoi volume data (green) resulting from case 1 is plotted against the $3\mathrm{\Gamma }$ distribution (blue). Values are shown upstream of the shock at the initial volume fraction ${\varphi }_0=0.1$ and downstream of the shock for various coefficients of restitution at the equilibrium volume fraction ${\varphi }_e=0.554$. Upstream of the shock, the histogram-obtained DEM data and the $3\mathrm{\Gamma }$ distribution are in excellent agreement. However, downstream of the shock, there is a departure from the $3\mathrm{\Gamma }$ distribution for all coefficients of restitution including the elastic limit $\left(e=1\right)$. As inelasticity is increased, the DEM data shift toward smaller Voronoi volumes.

Figure 11

Histogram Voronoi volume data (green) resulting from case 2 is plotted against the $3\mathrm{\Gamma }$ distribution (blue). Values are shown upstream of the shock at the initial volume fraction ${\varphi }_0=0.05$ and downstream of the shock for various coefficients of restitution at the equilibrium volume fraction ${\varphi }_e=0.25$. As inelasticity is increased, there is a much more pronounced shift than was observed in Fig. 10.

Figure 12

Histogram Voronoi volume data (green) resulting from case 3 is plotted against the $3\mathrm{\Gamma }$ distribution (blue). Values are shown upstream of the shock at the initial volume fraction ${\varphi }_0=0.1$ and downstream of the shock for various coefficients of restitution at the equilibrium volume fraction ${\varphi }_e=0.295$.

Figure 13

Sinusoidally varying plane average of the $3\mathrm{\Gamma }$ distribution for three volume fractions corresponding to the equilibrium volume fractions ${\varphi }_e$ downstream of the shock for cases 1(a), 2(b), and 3(c). A number of values for the amplitude are shown. As the amplitude of the wave increases the distribution shifts toward lower volumes similar to the downstream particles observed in Figs. 10, 11, and 12.

Figure 14

Case 1: Variation in granular temperature scaled by the peak granular temperature as well as the $x$ location of the peak. The elastic limit $\left(e=1\right)$ is shown to break the self-similar collapse.

Figure 15

Case 1: The peak nondimensional granular granular temperature (a) is shown to increase while the $x$ location of the peak (b) moves toward the shock as the coefficient of restitution approaches 1.

Figure 16

Case 1: The Voronoi volume fraction variance is plotted for three values of $e$, as well as the theoretical $3\mathrm{\Gamma }$ distribution. The variance is scaled by the initial upstream Voronoi variance present in the DEM simulations. For all observed restitution coefficients, the variance first increases for intermediate volume fractions during the compaction of the particles and settles to a lower value in the equilibrium regime. Decreasing the coefficient of restitution is shown to raise the final resting variance.

Figure 17

Case 1: The ensemble-averaged scaled velocity (a) and the ensemble-averaged volume fraction (b). The coefficient of restitution is shown to have little to no impact upon the mean momentum transfer in the streamwise direction.

Figure 18

Case 1: Individual components of the granular temperature are plotted showing that the temperature in the $x$ direction (a) carries between $56\%$ to $82\%$ of the velocity fluctuations, while the $y$ and $z$ temperatures [(b) and (c)] carry $9\%$ to $22\%$ of the velocity fluctuations.

Figure 19

Case 2: The mean velocity relaxation for ${\varphi }_0=0.05$ is shown to be the same for all coefficients of restitution.

Figure 20

Case 2: The scaled granular temperature variation is shown for ${\varphi }_0=0.05$. Unlike in case 1, all coefficients of restitutions collapse into a self-similar curve.

Figure 21

Case 2: Peak nondimensional granular granular temperature (a) and the nondimensional $x$ location of the peak (b).

Figure 22

The scaled Voronoi volume fraction variance for ${\varphi }_0=0.05$ for cases 2 (a) and 3 (b) plotted for three values of $e$ as well as the theoretical $3\mathrm{\Gamma }$ distribution. Unlike in case 1, the volume fraction variance is monotonically increasing for all values of $e$ as well as the $3\mathrm{\Gamma }$ distribution.

Figure 23

Case 3: The mean velocity relaxation for $\delta =0$ is shown to be the same for all coefficients of restitution.

Figure 24

Case 3: (a) The scaled granular temperature variation for $\delta =0$ shows a self-similar collapse under the prescribed scaling. The components of the granular temperature are scaled by the peak total values and plotted for the case where $\delta =0$. The trends in the $x$ temperature (b) and the $y$ temperature (c) are shown to follow that of case 1 where the relative contributions in the $x$ direction increase as the coefficient of restitution is lowered while the $y$ direction contributions increase as the decrease as the coefficient of restitution is lowered.

Figure 25

Case 3: For the $\delta =0$ case, the peak granular temperature (a) has an almost parabolic dependence upon $e$ first decreasing and then increasing as $e$ is increased. The location of the peak value (b) follows the trend of case 1, where it approaches the shock location as the coefficient of restitution is increased.