Transient granular shock waves and upstream motion on a staircase

A granular cluster, placed on a staircase setup, is brought into motion by vertical shaking. Molecular dynamics simulations show that the system goes through three phases. After a rapid initial breakdown of the cluster, the particle stream organizes itself in the form of a shock wave moving down the steps of the staircase. As this wave becomes diluted, it transforms into a more symmetric flow, in which the particles move not only downwards but also toward the top of the staircase. This series of events is accurately reproduced by a dynamical model in which the particle flow from step to step is modeled by a flux function. To explain the observed scaling behavior during the three stages, we study the continuum version of this model (a nonlinear partial differential equation) in three successive limiting cases. (i) The first limit gives the correct $t^{-1/3}$ decay law during the rapid initial phase, (ii) the second limit reveals that the transient shock wave is of the Burgers type, with the density of the wave front decreasing as $t^{-1/2}$, and (iii) the third limit shows that the eventual symmetric flow is a slow diffusive process for which the density falls off as $t^{-1/3}$ again. For any finite number of compartments, the system finally reaches an equilibrium distribution with a bias toward the lower compartments. For an unbounded staircase, however, the $t^{-1/3}$ decay goes on forever and the distribution becomes increasingly more symmetric as the dilution progresses.

Figure 1

(Color online) Four snapshots of a partial view of the staircase system, showing compartment 60–66 of the 1000 connected compartments used in the MD simulations [22]. The setup is vibrated in the vertical direction at frequency $f=100\text{Hz}$ and amplitude $a=0.2\text{mm}$, causing the particles to form a granular gas. At $t=0\text{s}$ (not shown), the experiment starts with 1000 particles distributed equally over compartments 50 and 51. (a) At $t=5.00\text{s}$, the downward-moving particle front has reached the first of the displayed compartments. (b) At $t=7.50\text{s}$, it is traveling further to the right. (c) At $t=10.00\text{s}$, the front has just passed through compartment 66. (d) Afterwards, as exemplified by the snapshot at $t=15.00\text{s}$, it leaves a trail of particles which becomes more dilute as time progresses.

Figure 2

(Color online) (a) The particle distribution $n_k\left(t\right)$ obtained from MD simulations with 1000 particles at 20 successive times starting from a cluster in compartments 50 and 51 [$n_{50}\left(0\right)=n_{51}\left(0\right)=0.50$, while the remainder of the $K=1000$ compartments is initially empty], at constant intervals of 2.0 s. The shaking parameters are the same as in Fig. 1. (b) The particle distribution $n_k\left(t\right)$ on a much longer time scale at intervals of 400 s. To reduce noise, a running average over 21 compartments has been taken (except for the initial distribution at $t=0\text{s}$). Note that a substantial—and increasing—amount of granular matter flows upwards, toward the left of compartment 50. (c) Doubly logarithmic plot of the maximum of the particle distribution $\left[n_{\text{max}}\left(t\right)\right]$ vs time. The three successive regimes (initial breakdown, shock wave, and symmetric diffusion) are indicated by dashed straight lines with slopes $-0.33$, $-0.63$, and $-0.31$, respectively.

Figure 3

(Color online) The flux functions toward the right [$F_R\left(n_k\right)$, Eq. (4)] and the left [$F_L\left(n_k\right)$, Eq. (5)] as function of $n_k$, i.e., the fraction of particles in the $k\text{th}$ compartment. In the dilute limit $n_k\to 0$, the fluxes start out equal, but for larger $n_k$, the difference between them rapidly diverges. In this plot, we have taken $A=1{\text{s}}^{-1}$ and $B_L=1000$.

Figure 4

(Color online) (a) The particle distribution $n_k\left(t\right)$ calculated with the flux model Eq. (6) at successive times $t=0,20,40,\dots ,400\text{s}$ starting from a cluster in compartments 50 and 51 $\left[n_{50}\left(0\right)=n_{51}\left(0\right)=0.50\right]$, while all the other of the $K=1000$ compartments are initially empty. The shaking parameters are $B_L=1000$ and $A=1{\text{s}}^{-1}$. (b) Doubly logarithmic plot of the maximum of the particle distribution $\left[n_{\text{max}}\left(t\right)\right]$ vs time. Indicated are the initial breakdown stage (characterized by a slope $-1/3$), the Burgers regime (with slope $-1/2$), and the diffusive regime (slope $-1/3$). The oscillatory structures visible for small $t$ reflect the filling and subsequent emptying of successive compartments that are, each in its turn, the best filled one; at the cusps, the maximum jumps from one compartment to the next.

Figure 5

(Color online) Logarithm of the estimated ratio $R$ between the first- (shock wave) and second-order (diffusive) terms in Eq. (10), as a function of the profile height $n_{\text{max}}$ for ${\tilde{B}}_L=1000{\text{m}}^2$ and ${\gamma }^{\prime }=20$ (solid black line, corresponding to the MD simulations of Fig. 2 and the flux model results of Fig. 4). For comparison, the gray curve represents the case for ${\tilde{B}}_L=10000{\text{m}}^2$ and ${\gamma }^{\prime }=20$. At this higher ${\tilde{B}}_L$ value, the first-order term has gained in importance. The range of $n_{\text{max}}$ for which ${\text{log}}_{10}R\left(n_{\text{max}}\right)>1$ increases, as well as the magnitude of $R$ itself, i.e., the shock wave lives longer and is also more pronounced.

Figure 6

(Color online) [(a) and (b)] The distributions at $t=0.1$, 0.2, 0.3, and 0.4 s in the numerical simulations, with $B_L=1000$ and $A=1{\text{m}}^{-1}$ just as in Fig. 4, both (a) unscaled (where we have also included the initial condition at $t=0\text{s}$) and (b) rescaled as $0.5n_k\left(t\right)t^{1/3}$ vs $\left(k-k_0\right)t^{-1/3}$ [see Eq. (24)]. Here, $k_0=50$ indicates the position of the initial cluster. The cluster rapidly breaks down toward the right with the same exponent 1/3 as in the symmetric diffusion stage. (c) Solution $H\left(\xi \right)$ of Eq. (23) for $\Gamma \approx 0.5$, representing the distribution for small times $t$.

Figure 7

(Color online) (a) The rescaled distributions during the stage when the profile has the triangular form of a Burgers shock wave in the numerical simulations, $t=40–220\text{s}$, at time intervals of 20 s (i.e., at much later times than in Fig. 6). A running average over three compartments has been taken to smoothen the profiles to some extent. The density is rescaled as $n_k\left(t\right)t^{1/2}$ and the position along the staircase as $\left(k-k_1\right)t^{-1/2}$ [Eq. (30) with $A=1{\text{s}}^{-1}$], where $k_1=50$ indicates the compartment in which the tail of the triangular profile is positioned, which happens to coincide with the compartment in which we placed the initial cluster. (b) Idem calculated with the flux model Eq. (6), under the same conditions as in Fig. 4, for $t=100–250\text{s}$ at time intervals of 25 s.

Figure 8

(Color online) (a) The distributions at $t=420$, 910, 1410, 1910, and 2410 s [i.e., much later than in Fig. 7a] in the MD simulations, rescaled as $n_k\left(t\right)t^{1/3}$ vs $\left(k-k_2\right)t^{-1/3}$. To smoothen the profiles, a running average over 21 compartments has been taken. The fact that the best fit is obtained for $k_2=45$, i.e., five compartments to the left of the original cluster, reflects the transition from shock wave to diffusive regime, in which particles also move upwards. The upward twist at the left of the profiles shows that material is heaping up in the leftmost compartments. (b) Idem calculated with the flux model Eq. (6), under the same conditions as in Fig. 4, at times $t=3000–23000\text{s}$ at time intervals of 2500 s, and rescaled as $n_k\left(t\right)t^{1/3}$ vs $\left(k-k_2\right)t^{-1/3}$ with $k_2=43$. Although the diffusion is governed by the same Eq. (23) as the initial breakdown, the profiles differ considerably since the two regimes start out from very different situations (cf. Fig. 6).

Figure 9

(Color online) (a) Apparent but untrue self-similarity. Solution $H\left(\xi \right)$ of Eq. (23) for $\Gamma =-0.5$. Although the shape bears a marked resemblance to that of the MD results of Fig. 8a (gray shaded curves), the depicted coincidence has only been obtained via a nonadmissible translation of the MD profiles to the left. (b) Nonconspicuous but true self-similarity. Rescaled flux model profiles obtained in a much larger system ($10000$ compartments with the cluster initially in compartment 4000 and its neighbor) and on a much longer time scale than in Fig. 4, namely, at $t={10}^3$, ${10}^5$, ${10}^7$, and ${10}^9\text{s}$. The values of $B_L=1000$ and $A=1{\text{s}}^{-1}$ are unaltered. These profiles reveal a slow convergence toward the symmetric solution of Eq. (23), i.e., the one with $\Gamma =0$ (see also the Appendix).

Figure 10

Approach toward the final state on a staircase of finite length. Starting out from a dilute homogeneous distribution over a staircase with $K=1000$ steps (and shaking parameters $B_L=1000$, $A=1{\text{s}}^{-1}$), we witness how the system slowly evolves toward a biased distribution in which the fluxes between all neighboring compartments are in equilibrium. In the last snapshot $\left(t={10}^9\text{s}\right)$, the convergence is complete, i.e., the density profile will not visibly change anymore. This profile is accurately described by Eq. (36) with $C=1.56{10}^6$.

Figure 11

(Color online) Decay of $n_{\text{max}}\left(t\right)$ calculated from the flux model Eq. (6) for four different values of $B_L$, each one starting from the same initial condition $n_{50}\left(0\right)=n_{51}\left(0\right)=0.50$ in a system of $K=1000$ compartments. (a) For $B_L=1$ (green curve), the breakdown of the cluster directly goes over into symmetric diffusion, without intervention of any shock wave, hence the slope is permanently close to $-1/3$. The deviation beyond ${\text{log}}_{10}t=4.3$ is caused by the left boundary, which is reached relatively soon at this strong shaking. [(b) and (c)] For $B_L=300$ and ${10}^4$ (blue and red curves) all three regimes are distinguishable [(i) initial breakdown with slope $-1/3$; (ii) shock wave with slope $-1/2$; (iii) diffusion with slope $-1/3$], with the length of the shock wave regime growing for increasing $B_L$. (d) For $B_L={10}^6$ (black curve), the shock wave regime extends all the way to the end of the depicted time interval, as the crossover takes place at a value $n_{\text{ii}\to \text{iii}}<{10}^{-3}$, in agreement with Eq. (38). Dashed lines with slopes $-1/3$, $-1/2$, and (again) $-1/3$ have been added for reference.

Figure 12

Starting out with a large heap on the steps $k=12,\dots ,17$ and a small one at $k=34,\dots ,39$ (see the snapshot at $t=0\text{s}$), one witnesses first a rapid breakdown and formation of two Burgers-like shock waves $\left(t=80\text{s}\right)$. The higher wave travels faster and overtakes the smaller one $\left(t=170\text{s}\right)$, leading to a single density profile that diffuses out into both directions $\left(t=1300\text{s}\right)$. Eventually, the upstream flow toward the left becomes almost equal to the flow toward the right. At $t=5000\text{s}$, the compartments 1–11 have already received a substantial amount of material. The shaking parameters are $B_L=2250$ [41] and $A=1{\text{s}}^{-1}$, which makes the dynamics of the large cluster directly comparable to Figs. 3, 4, 5, 6a, 6b, 7b, 8b, 9b. The velocity of the large shock wave at $t=80\text{s}$ is 0.100 compartments/s, against 0.077 compartments/s for the smaller one (indicated by the arrows).

Figure 13

(Color online) Solutions of Eq. (A1), or Eq. (23) in the main text, normalized on the interval $-10\le \xi \le 10$. Dashed black curve is the symmetric case for $\Gamma =0$, which can be solved analytically. Solid blue, red, and black curves that deviate increasingly from the symmetric solution correspond to $\Gamma \approx 0.01$, $\Gamma \approx 0.11$, and $\Gamma \approx 0.36$, respectively. These were determined numerically. Negative values of $\Gamma$ yield solutions that are the mirror images (with respect to $\xi =0$) of those for positive $\Gamma$.