Explosive dispersal of granular media

Abstract

Explosive dispersal of granular media widely occurs in nature across various length scales, also enabling engineering applications ranging from commercial or military explosive systems to the loss prevention industry. However, the complex particle–flow coupling makes the explosive dispersal behaviour of particles difficult to control or even characterize. Here, we study the central explosion-driven dispersal of dense particle layers using the coarse-grained computational fluid dynamics–discrete element method and present a comprehensive investigation of both macroscale dispersal behaviours and particle-scale pattern formation. Employing three independent dimensionless parameters that characterize the efficiency, homogeneity and completeness of explosive dispersal, we categorize the dispersal behaviours into ideal, partial, retarded and failed modes, and propose the corresponding thresholds. As the mass ratio of granular materials to central pressurized gases (M/C) spans four orders of magnitude, the dispersal mode transitions from ideal to partial, then to retarded and finally to failed mode. The transitions of dispersal modes correspond to the particle–flow coupling regime crossovers, which change from decoupling to weak, medium and finally to strong coupling as the dispersal mode undergoes corresponding transitions. We proceed to develop continuum models accounting for the shock compaction and the ensuing pulsation of the particle ring that are capable of identifying the ideal dispersal mode from various dispersal systems. We also provide insights into the origins of diverse particle-scale patterns that are strongly correlated with macroscale dispersal modes and critical for the accurate prediction of dispersal modes.

Explosive dispersal of granular media, whereby the detonation of a central explosive or sudden release of pressurized gases disperses densely packed particles to form dilute particle clouds, occurs in a wide range of natural phenomena and engineering processes (Formenti, Druitt & Kelfoun 2003; Eckhoff 2009; Aglitskiy et al. 2010; Kuranz et al. 2018; Frost 2018). In volcanic eruptions, the explosive expansion of pressurized gases through an initially concentrated dispersion of particles expels mixtures of pressurized gases and fragments of magma within volcanic conduits (Marjanovic et al. 2018). For various commercial and military explosive systems, such as fire extinguishing devices using explosive dispersed dry powders (Klemens, Gieras & Kaluzny 2007), thermobaric and fuel-air bombs (Frost, Goroshin & Zhang 2010; Frost et al. 2012), the explosive dispersal of granular media and the ensuing mixing of particulate matter with air are of particular importance to their reliable applications (Zhang et al. 2015; Bai et al. 2018; Posey et al. 2021). Another prominent application involves blast mitigation by surrounding the explosive with a layer of particles (Milne et al. 2014; Pontalier et al. 2018). The transfer of heat and energy from explosives to the particles during shock compaction and the ensuing gas–particle interaction significantly reduce the blast overpressure and impulse.

Considerable effort has been devoted to understanding the explosive dispersal process of granular media, which involves complex multiphase coupling and spans a range of scales (Frost et al. 2010; Rodriguez et al. 2013, 2014; Xue, Sun & Bai 2016; Osnes, Vartdal & Reif 2017; Frost 2018; Kun et al. 2018; Xue et al. 2018; Carmouze et al. 2019; Chiapolino & Saurel 2019; Fernandez-Godino et al. 2019; Mo et al. 2019; Koneru et al. 2020). In particular, most studies focus on the hierarchical jet-like structures observed in cylindrical and spherical geometries (Frost et al. 2010; Rodriguez et al. 2013, 2014; Xue et al. 2016; Osnes et al. 2017; Frost 2018; Kun et al. 2018; Xue et al. 2018; Carmouze et al. 2019; Chiapolino & Saurel 2019; Fernandez-Godino et al. 2019; Mo et al. 2019; Koneru et al. 2020). The origin and evolution of particle jetting have been extensively investigated and found to be closely associated with the spatial distribution of disseminated particles and the terminal velocities (Kandan et al. 2017; Loiseau et al. 2018; Pontalier et al. 2018; Li et al. 2022), two prominent attributes relevant to pertinent applications. However, some fundamental questions that are of much interest to engineering applications remain left unaddressed. Specifically, how fast is dispersal completed? Are the dispersed particles homogeneously distributed in space? Is any nontrivial proportion of particles not expelled successfully? These three questions are critical to assessing three fundamental attributes of explosive dispersal, namely, efficiency, homogeneity and completeness. Since the explosive dispersal processes for systems with vastly varying scales and structural parameters may occur on markedly different time and length scales, adequately quantifying the key attributes is quite challenging, let alone properly characterizing distinctly different explosive dispersal processes. Here, by identifying and understanding the most significant events contributing to each respective attribute, we construct three dimensionless parameters to categorize various dispersal processes into distinct dispersal modes.

The complex particle-flow interaction plays an important role in the explosive dispersal of particles, especially the formation of particle jetting (Xu et al. 2013; Osnes et al. 2017; Carmouze et al. 2019; Chiapolino & Saurel 2019; Fernandez-Godino et al. 2019; Mo et al. 2019; Koneru et al. 2020; Li et al. 2022). Based on particle resolved numerical investigations, Xu et al. (2013) and Mo et al. (2019) observed that the microgas jets that formed at the internal surface of the particle ring penetrated the particle layer through gaps created by the inelastic collision between particles, generating high-speed particle agglomerates. Recently, shock-driven multiphase instability (SDMI), a variant of the classic Richtmyer–Meshkov (RM) instability, was reported to be responsible for particle jetting in dilute systems (Osnes et al. 2017; Chiapolino & Saurel 2019; Koneru et al. 2020). In contrast, another interfacial instability, shock-driven granular instability (SDGI), is proposed to account for the jetting of shock-loaded dense granular interfaces (Li et al. 2022). Although these particle-flow interaction mechanisms provide important insights into the physics underlying particle dispersal, they fail to incorporate the coupling between the expanding ring as an entirety and the evolving central pressure field, which governs the dynamics of the particle ring. One striking phenomenon resulting from macroscale particle-flow coupling is the pulsation of the particle ring analogous to the pulsation of the gas bubble generated in the underwater explosion (Wang et al. 2021). The pulsation of the particle ring sustained by the central explosion, albeit rarely reported, is essential to understand the overall dispersal behaviour of the particle ring. The emergence, prevalence and waning of the particle ring pulsation must be understood based on the varying particle-flow coupling regimes, which can be quantified by calculating the characteristic time scale ratios of signature events.

As highlighted in many investigations involving shock/blast-driven particle-laden compressible flows, understanding or even predicting the behaviour of these flows requires knowledge of particle-scale phenomena, such as particle‒particle collisions, particularly in dense particle flows, shock-particle interactions and the coupling between particle and flow jetting (Osnes et al. 2017; Sundaresan, Ozel & Kolehmainen 2018; Carmouze et al. 2019; Chiapolino & Saurel 2019; Mo et al. 2019; Koneru et al. 2020; Tian et al. 2020). Thus, we conducted four-way coupled two-dimensional (2-D) Euler‒Lagrange simulations of the explosive dispersal of the particle ring wherein the denotation of the cylindrical burster was simulated by the sudden release of pressurized gases in the central gas pocket. A numerical protocol was established that allowed us to convert the dispersal process driven by the central detonation to the equivalent process driven by the pressurized gases and vice versa. Using this protocol, we performed a range of simulations for dispersal systems with varying structures, including the energy of pressurized gases, the mass of the particle ring, and the packing fraction of particles. Regardless, systems with the same mass ratio of particles to central gases, M/C = m_ring/m_gas, display similar dispersal behaviours, highlighting the role of M/C as one primary governing parameter. We developed theoretical models based on the continuum assumption to account for the shock compaction phase and the ensuing pulsation of the particle ring and to predict the dispersal mode for a particular dispersal system. The combined models are capable of identifying the ideal dispersal mode while falling short of recognizing the other three modes, which require knowledge of grain-scale particle loosening mechanisms.

The paper is organized as described below. In § 1, the numerical method is presented, followed by the description of the numerical setup in § 2. A variety of macroscale behaviours of dispersed particle rings are elaborated in § 3, which are characterized using three prominent attributes and categorized into different modes. The macroscale particle-flow coupling and its correlation with the dispersal mode are explored in § 4. The dominant mesoscale dispersal structures are discussed in § 5. Finally, the results are summarized in § 6.

1. Numerical method

Numerical simulations were performed based on coarse-grained compressible computational fluid dynamics–discrete parcel method (CCFD-DPM), a coarse-grained Euler–Lagrange numerical approach suitable for gas–particle flows in laboratory-scale systems (Sundaresan et al. 2018; Tian et al. 2020). The CCFD-DPM approach tracks and accounts for parcel–parcel contact interactions. Each parcel consists of multiple physical particles with the same physical and kinetic properties. The number of real particles that represent a computational parcel is quantified using a scaling factor called the super particle loading, χ, whose value is set based on the volume/mass fraction of the particles and computational memory available. For particle-gas systems, the value of χ reported in the previous literature ranges from O(10¹) to O(10³) (Osnes et al. 2017; Koneru et al. 2020). In the present study, χ is of the order of O(10¹).

For the gas phase, the volume-averaged governing equations ((1.1)–(1.3)) constructed in the Eulerian frame are based on a five-equation transport model, i.e. a simplified form of the Baer–Nunziato (B-N) model, which has been modified to account for compressible multiphase flows ranging from dilute to dense gas–particle flows (Baer & Nunziato 1986):

(1.1)\begin{gather}\frac{{\partial ({\phi _f}\langle {\rho _f}\rangle )}}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }({\phi _f}\langle {\rho _f}\rangle {\boldsymbol{\tilde{u}}_f}) = 0,\end{gather}

(1.2)\begin{gather}\begin{array}{ccccc} \dfrac{{\partial ({\phi _f}{{\boldsymbol{\tilde{u}}}_f}\langle {\rho _f}\rangle )}}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }({\phi _f}\langle {\rho _f}\rangle {{\boldsymbol{\tilde{u}}}_f}{{\boldsymbol{\tilde{u}}}_f} + {\phi _f}\langle {P_f}\rangle )\\ = \langle {P_f}\rangle \boldsymbol{\nabla }{\phi _f} + \sum\limits_i {\{ {\phi _{p,i}}{\rho _{p,i}}{D_{p,i}}({{\boldsymbol{\tilde{u}}}_{p,i}} - {{({{\boldsymbol{\tilde{u}}}_f})}_i})\} } , \end{array}\end{gather}

(1.3)\begin{gather}\begin{array}{ccccc} \dfrac{{\partial ({\phi _f}\langle {\rho _f}\rangle {{\tilde{E}}_f})}}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }({\phi _f}\langle {P_f}\rangle {{\tilde{E}}_f}{{\boldsymbol{\tilde{u}}}_f} + {\phi _f}\langle {\rho _f}\rangle {{\boldsymbol{\tilde{u}}}_f})\\ = \langle {P_f}\rangle \boldsymbol{\nabla }{\phi _f} \boldsymbol{\cdot} {{\boldsymbol{\tilde{u}}}_p} + \sum\limits_i {\{ {\phi _{p,i}}{\rho _{p,i}}{D_{p,i}}({\boldsymbol{u}_{p,i}} - {{({{\boldsymbol{\tilde{u}}}_f})}_i}) \boldsymbol{\cdot} {{\boldsymbol{\tilde{u}}}_{p,i}}\} } . \end{array}\end{gather}

The gas volume fraction, velocity, density, pressure and total energy of the gas are represented by ϕ_f, u_f, ρ_f, P_f and E_f, respectively, where E_f = ρ_fe_f + 0.5 u_f u_f. In (1.1)–(1.3), $\left\langle {} \right\rangle $ and $\widetilde {}$ denote phase-averaged and mass-averaged variables, respectively, ρ_p.i and u_p,i represent the density and velocity of parcel i, respectively, D_p,i is the drag force coefficient of parcel i and ϕ_p.i = w_i,fV_p,i/V_f is the contribution of parcel i to the weighted particle volume fraction (w_i,f is the distributed weight that parcel i contributes to the particle volume fraction in fluid cell f, and V_p,i and V_f are the volumes of parcel i and the fluid cell, respectively).

We employ the Di Felice model combined with Ergun's equation ((1.4)–(1.7)) to calculate D_p (Felice 1994), which considers the effects of both the particle Reynold number, Re_p, and the particle phase volume fraction, ϕ_p, and has been widely used in particle-laden multiphase flows:

(1.4)\begin{gather}{D_{p,i}} = \frac{3}{{8{\rho _p}/{\rho _f}}}{C_d}\frac{{|{\boldsymbol{u}_f} - {\boldsymbol{u}_{p,i}}|}}{{{r_p}}},\end{gather}

(1.5)\begin{gather}{C_d} = \frac{{24}}{{R{e_p}}}\left\{ {\begin{array}{@{}ll} {8.33\dfrac{{{\phi_p}}}{{{\phi_f}}} + 0.0972R{e_p}}&{\textrm{if}\;{\phi_f} < 0.8}\\ {{f_{base}} \phi_f^{ - \zeta }}&{\textrm{if}\;{\phi_f} \ge 0.8} \end{array}} \right.,\end{gather}

(1.6)\begin{gather}{f_{base}} = \left\{ {\begin{array}{@{}ll} {1 + 0.167Re_p^{0.687}}&{\textrm{if}\;R{e_p} < 1000}\\ {0.0183R{e_p}}&{\textrm{if}\;R{e_p} \ge 1000} \end{array}} \right.,\end{gather}

(1.7)\begin{gather}\zeta = 3.7 - 0.65\exp [ - {\textstyle{1 \over 2}}{(1.5 - {\log _{10}}R{e_p})^2}].\end{gather}

In (1.4) and (1.5), C_d is the dimensionless coefficient of the drag force and r_p is the particle radius. For dense particle flows (ϕ_f < 0.8), (1.4) reduces to the original Ergun equation. Otherwise, C_d takes the form of Stokes's law multiplied by a factor of f_base, which varies with Re_p, as indicated by (1.6) and (1.7).

The particle phase is represented by discrete parcels whose motion is governed by Newton's second law ((1.8) and (1.9)):

(1.8)\begin{gather}\frac{{\textrm{d}{\boldsymbol{u}_{p,i}}}}{{\textrm{d}t}} = {D_{p,i}}({\boldsymbol{u}_f} - {\boldsymbol{u}_{p,i}}) - \frac{1}{{{\rho _p}}}\boldsymbol{\nabla }\langle {P_f}\rangle + \frac{1}{{{m_p}}}\sum\limits_j {{\boldsymbol{F}_{C,ij}}} ,\end{gather}

(1.9)\begin{gather}\frac{{\textrm{d}\kern0.7pt{\boldsymbol{x}_{p,i}}}}{{\textrm{d}t}} = {\boldsymbol{u}_{p,i}},\end{gather}

where u_p,i and x_p,i denote the velocity and displacement of parcel i, respectively, and F_C,ij represents the collision force between parcels i and j.

A four-way coupling strategy was adopted to account for the momentum and energy transfer between gases and particles. Specifically, the particle drag force and the associated work were incorporated into the momentum and energy equations of the gas phase as the source terms. The particle parcels are driven by the pressure gradient force, drag force and collision force between parcels (1.8). A soft sphere model, which is represented by a linear-spring dashpot model, was employed to model the collision force between parcels. Hence, F_C,ij consists of a repulsive force and a damping force:

(1.10)\begin{equation}{\boldsymbol{F}_{C,ij}} = {k_{n,p}}{\boldsymbol{\delta}_n} - {\gamma _{n,p}}{\boldsymbol{u}_{n,ij}},\end{equation}

where k_n,p and γ_n,p are the stiffness and damping coefficients of parcels, respectively, and δ_n and u_n,ij are the overlap and difference in the normal velocity between parcels in contact, respectively. Here, γ_n,p is a function of the parcel restitution coefficient ε_p:

(1.11)\begin{equation}{\gamma _{n,p}} ={-} \frac{{2\ln {\varepsilon _p}}}{{\sqrt {{{\rm \pi}^2} + \ln {\varepsilon _p}} }}\sqrt {{m_p}{k_{n,p}}} .\end{equation}

The weighted essentially non-oscillatory (WENO) scheme was used to reconstruct the primary flow variables and solve the equations governing the gases. A Riemann solver proposed by Harten, Lax and van Leer was used to obtain the intercell fluxes. The third-order Runge–Kutta method was applied for time integration. The equations describing the parcel velocity and position were discretized using the velocity-Verlet algorithm. Bilinear/trilinear interpolation functions were adopted to calculate the particle volume fraction and source terms on the Eulerian grids, as well as the fluid variables on Lagrangian parcels (Liu, Osher & Chan 1994). Numerical details of CCFD-DPM are available in our previous publication (Tian et al. 2020). The CCFD-DPM framework introduced the variable A and has been validated in several benchmark experiments involving shock-driven particle-laden flows (Tian et al. 2020), such as experiments performed by Rogue et al. (1998), which investigate shock dissipation through particle curtains; the experiments conducted by Britan & Ben-Dor (2006), which assess both the gaseous and solid pressures inside particle columns impinged head on by shocks; and the experiments assessing shock-induced interfacial instability of granular media (Kun et al. 2018).

Notably, for dispersal systems where a central burster is enclosed by ductile or brittle particles, the deformation and/or fracture of particles and the sintering between particles are inevitable for particles adjacent to the burster. These phenomena may well affect the ensuing dispersal, particularly for systems with thin particle shells. Since the current CCFD-DPM framework does not account for these particle-scale phenomena, we must be aware of the scenario in which the simulations and the related discussions are valid, namely, either the explosive sources are so mild that the particle-scale phenomena are negligible, or particles are composed of materials with a combination of high compressive strength, high toughness and moderate/high hardness. Otherwise, the particle-scale phenomenon may cause the actual dispersal to deviate from that predicted by the framework in which these particle-scale phenomena are excluded.

2. Numerical set-up

The two-dimensional (2-D) configuration shown in figure 1(a) serves as an archetypical system to investigate the explosive dispersal behaviour of particle rings. Instead of the burster used in the explosive dispersal experiment, the numerical simulations employ a high-pressure gas pocket with a radius of R_{gas ,0} filled with air to disperse the enclosing ring consisting of 2-D disks with varying packing fractions, ϕ ₀. For the investigated dispersal systems described below, the initial radius and temperature of the gas pocket remain constant, R_{gas ,0} = 20 mm and T ₀ = T_amb = 298 K, respectively, while the initial pressure P ₀ varies. The radius expansion algorithm is used to generate the particle packing with a given ϕ ₀ in an annular domain with the inner and outer radii of R_{in ,0} and R_{out ,0}, respectively. A population of disk-like parcels with artificially small radii ensure no parcel overlap is randomly created within the specified domain. Then, all parcels are expanded until the specified parcel size distribution and desired packing fraction are both satisfied. The real 2-D particle has a diameter of 100 μm, while the diameter of the parcel uniformly ranges from 360 to 780 μm to avoid potential crystallization during shock compaction. Due to the use of computational parcels, the collision distance must be modified in the same manner as the parcel diameter. A random but homogenous arrangement of parcels is achieved, as shown in the zoomed-in inset of figure 1(a), where the parcels are coloured according to the local Voronoi packing fraction, ϕ_p,voro. The parcel has a density of 2500 kg m⁻³, which is the same as the real particles. The restitution coefficient, ε_p, is set to 0.6, which accounts for the energy dissipation inside the parcel; the normal stiffness of contacts between parcels is set to 2 × 10⁷ N m⁻¹.

Figure 1. (a) A quarter of the two-dimensional numerical configuration wherein the central gas pocket is enclosed by a particle ring consisting of polydispersed computational parcels. Inset shows the zoomed-in arrangement of particle packing. (b) Distribution of the structural parameters, P ₀, h, ϕ ₀ as well as the resulting M/C for dispersal systems which are to be numerically investigated. Here, h is the ring thickness, h = R_{out ,0} – R_{in ,0}.

Evidently, a host of structural variables characterize the dispersal system and play a role in the resulting dispersal behaviour with varying levels of relative importance. Instead of performing a conventional parametric study, we focus on their combined effects manifested by one single dimensionless parameter, the mass ratio, M/C, which defines the energy imparted to the confining materials per unit mass. Previous studies have accentuated the fundamental role played by the M/C in explosive dispersal by correlating the initial expansion velocity of disseminated particles with the M/C (Loiseau et al. 2018; Pontalier et al. 2018). The structural variables denoted in figure 1(a) all contribute to the mass ratio M/C, as expressed in (2.1):

(2.1)\begin{equation}M/C = \frac{{{\rm \pi}(R_{out,0}^2 - R_{in,0}^2){\phi _0}{\rho _p}}}{{{\rm \pi} R_{gas}^2{\rho _{gas}}}} = \frac{{(R_{out,0}^2 - R_{in,0}^2){\phi _0}{\rho _p}R{T_0}}}{{R_{gas}^2{P_0}}},\end{equation}

where R is the specific gas constant of air, $R=288.7\ {\rm J}/{\rm kg}/{\rm K}$. The mass ratio clearly does not reflect the nature of various explosive sources with vastly different specific energies. Hence, two dispersal systems with the same M/C but different explosive sources, for instance, the burster and the pressurized gas pocket, will undoubtedly exhibit distinct dispersal behaviours. In Appendix A, we propose a mass ratio conversion between systems with different explosive sources based on the energy equivalent principle whereby the dispersal behaviours exhibited by the systems with explosive sources other than pressurized gases are incorporated into the dispersal mode classification framework proposed in § 3.3.

As illustrated in figure 1(b), by varying the values of P ₀, R_{out ,0} and ϕ ₀, the mass ratio M/C for the simulated dispersal systems ranges from 9.7 to 6800, as marked in figure 1(b), spanning four orders of magnitude. Three different values of ϕ ₀, namely, 0.5, 0.6 and 0.65, were used to represent loose, dense and densest random packings, respectively. Notably, the packing fraction values presented throughout this work are those of the equivalent three-dimensional (3-D) packings derived from the conversion between the packing fractions in the 2-D and 3-D packings proposed by Borchardt-Ott (2012):

(2.2)\begin{equation}{\phi _{3\text{-D}}} = 0.2595 + \frac{{{\phi _{2\text{-D}}} - 0.0931}}{{0.2146 - 0.0931}}(0.4760 - 0.2595).\end{equation}

Hence, the 3-D particle packing fractions ϕ _3D = 0.5, 0.6 and 0.65 correspond to the 2-D particle packing fractions ϕ _2D = 0.77, 0.83 and 0.86, respectively. For each M/C, at least three systems with different combinations of P ₀, R_{out ,0} and ϕ ₀ were constructed to assess the variability of the results. In total, we performed explosive dispersal simulations for 70 different systems with 20 different values of M/C. For clarity, the system is labelled with four symbols, M/C, P ₀ (in units of bar), the thickness of the ring ($h = {R_{out,0}} - {R_{in,0}}$, in units of mm) and ϕ ₀.

3. Results

3.1. Macroscale dispersal behaviour

A space–time (R–t) diagram showing the spaciotemporal variations in the particle packing fraction ϕ_p(R, t) was constructed using the circumferentially averaged ϕ_p to understand the overall explosive dispersal behaviour of the particle ring. Figure 2(a–d) displays four typical R–t diagrams for systems with different M/C values. Upon the release of pressurized gases, an incident shock wave impinges on the inner surface of the particle ring (the dynamics of shocks are illustrated in figure 3a), resulting in a compaction front (CF) traversing the particle bed. Across the CF, ϕ ₀ immediately jumps to ϕ_comp, the packing fraction of the compacted particles. A rarefaction wave (RW) reflects off the external surface of the particle ring upon the arrival of the CF, resulting in a discernible decrease in ϕ_p in its wake. The dynamics of CF and the ensuing RW are sufficiently recognized in the R–t diagram of the volumetric strain rate ${\dot{\varepsilon }_v}$ in particles (the calculation of ${\dot{\varepsilon }_v}$ is presented in Appendix B), as plotted in figure 2(e).

Figure 2. R–t diagrams of the ϕ_p of particle rings in the dispersal system (a) 9.7-245-10-0.5, (b) 103.7-100-30-0.6, (c) 1024-20-50-0.6 and (d) 4875-20-140-0.6. (e) R–t diagrams of the volumetric strain rate ${\dot{\varepsilon }_v}$ in particles in the dispersal system 494-200-140-0.6. CF, compaction front; RW, rarefaction wave. The white and pink dashed lines indicate the boundaries delimiting the particle cloud, IPF and OPF, and dense core band (CB), IS_CB and OS_CB, respectively.

Figure 3. (a) Temporal variations in ϕ_p,ave for typical dispersal systems. (b) Scaled temporal variations in ϕ_p,ave for the same systems shown in panel (a) wherein the time is scaled by t_ring.

Despite similar early-time shock compaction phases, systems with different M/C ratios proceed with distinctly different late-time dispersal phases. The ring in system 9.7-245-10-0.5 undergoes persistent expansion, while particle shedding progressively disintegrates the ring from inwards to outwards (figure 2a). In contrast, as shown in figure 2(b), the ring in system 103.7-100-30-0.6 remains densely packed during expansion, which strikingly proceeds to undergo noticeable contraction accompanied by significant outbound particle jetting commencing from the outer surface (see the inset in figure 2b). The contraction observed in system 103.7-100-30-0.6 becomes more pronounced in system 1024-20-50-0.6, wherein a densely packed core band (CB) not only survives the entire contraction but also embarks on a second expansion (figure 2c). The internal and external surfaces of the densely packed CB, IS_CB and ES_CB are defined such that the packing fraction of bulk enclosed by IS_CB and ES_CB remains above 0.3. Only 5 % of the volume of particles resides inside the inner particle front (IPF) and another 5 % resides beyond the outer particle front (OPF). The trajectories of the IPF, OPF, IS_CB and OS_CB are superimposed in figures 2(c) and 2(d).

At the onset of the second expansion of the CB, inbound particle jetting initiated from the IS_CB occurs (figures 2c and 2d). In contrast to the internal particle shedding during the first ring expansion, as observed in system 9.7-245-10-0.5 (figure 2a), where the loose particles still move outwards, the particles entrained by the inbound particle jets travel inwards, eventually moving to the central area. The relatively short-lived CB in system 1024-20-50-0.6 becomes enduring in system 4875-20-140-0.6 and we observe multiple pulsations of CB with decreasing amplitudes (see figure 2d). Notably, inbound and outbound particle jetting occur at each contraction-to-expansion and expansion-to-contraction turning point, respectively. Particle shedding, and inbound and outwards particle jetting are elaborated in § 5.

3.2. Characterization of dispersal behaviour

The diverse dispersal behaviours of particles shown in figures 2(a)–2(d) lead to a fundamental question: Can we properly characterize and classify the diverse dispersal processes? This question is critical for adequately evaluating the dispersal process and tailoring the dispersal systems for a specific application. In this section, we aim to address this question by establishing a unified characterization framework.

Unarguably, efficiency is among the most important properties of the dispersal process from an engineering perspective. The time required for the packed particles to be dispersed into a dilute state, t_dis, is a useful indicator of the dispersal efficiency, although t_dis scales with the length scale of the system. Here, the dispersed ring is referred to as dilute when the average particle volume fraction in the annular region delimited by IPF and OPF, ϕ_p,ave, is less than 0.1. Figure 3(a) plots the variations in ϕ_p,ave in several typical dispersal systems. Notably, one or more kinks are commonly observed in ϕ_p,ave(t) curves, which are caused by the sudden deflections of IPF or OPF due to the burst of inbound or outbound particle jetting. As expected, systems with an increased M/C, due to heavier particle rings, smaller input energies or both, take more time to be dispersed to the nominal dilute state. As indicated in figure 4(a), t_dis increases logarithmically with the M/C and approaches infinity when the M/C is greater than 2000. For systems with an M/C of O(10³), ϕ_p,ave does not reach 0.1 during the computational time. Hence, t_dis was determined by extrapolating the ϕ_p,ave(t) curves, as described in Appendix C.

Figure 4. Variations in (a) t_dis, (b) t_ring, (c) t_dense with M/C. The dashed portions of t_dis(M/C) and t_dense(M/C) curves at the upper limit of the M/C range suggest t_dis and t_dense herein approach infinity rather than have definite values.

We instead use a dimensionless parameter ξ as the efficiency indicator by scaling t_dis with the characteristic time of the ring expansion, t_ring, ξ = t_dis/t_ring, to properly evaluate the efficiency of dispersal processes in systems with varying M/C values. Here, t_ring is the time it takes the ring to expand to twice the initial diameter, namely, t_ring = R_{out ,0}/V_out, with V_out representing the initial surface expansion velocity of the ring. The M/C dependence of t_ring is plotted in figure 4(b). Figure 3(b) depicts the variations in ϕ_p,ave with scaled time t/t_ring, with significantly reduced deviations observed between ϕ_p,ave(t) curves for different systems. Figure 5 shows the variation in ξ with M/C, and the error bars represent the variability between systems with different combinations of P ₀, R_{out ,0} and ϕ ₀. For values less than a threshold mass ratio (M/C)_ξ ~ 2500, ξ randomly fluctuates below 10. Afterwards, ξ increases above 10 and approaches infinity towards the upper limit of the M/C range. Therefore, the dispersion is at least one order slower than the ring expansion. For efficient dispersal, the dispersion should be completed during the time scale of ring expansion, leading to the criterion for the efficiency:

(3.1)\begin{equation}\xi < 10.\end{equation}

Another issue that is equally important is the spatial concentration of the dispersed particles, which is strongly affected by the emergence of densely packed CB, as shown in figures 2(c) and 2(d). If the CB survives for the entire dispersal time, t_dis, the bulk of particles will adhere inside the narrow CB, although the nominal dilute state has been achieved. This dispersal process, albeit efficient, is highly inhomogeneous. Figure 4(c) presents the variation in t_dense with M/C, which is fitted by an exponential function asymptotically approaching infinity. Rings in the systems with an M/C ratio larger than 1600 undergo slowly decaying pulsation and eventually dwell in an annular region located not far from its initial position, leading to an infinite t_dense. The derivation of t_dense for the systems with the enduring CB is described in detail in Appendix C. The ratio between the CB survival duration to the dispersal time, which is defined as κ = t_dense/t_dis, produces a homogeneity indicator. The M/C dependence of κ plotted in figure 5 shows a similar trend with ξ(M/C). Here, κ remains below unity until M/C exceeds a threshold, (M/C)_κ ~ 500, beyond which κ increases above unity, substantially increases with M/C and eventually approaches infinity at the upper limit of the M/C range. A homogeneous dispersal does not allow the presence of the CB at the end of the dispersal. Accordingly, the criterion for homogenous dispersal is expressed as follows:

(3.2)\begin{equation}\kappa < 1.\end{equation}

Figure 5. M/C dependences of dimensionless parameters ξ, κ and χ. The critical M/C values corresponding to the thresholds for the efficiency, homogeneity and completeness, (M/C)_ξ, (M/C)_κ and (M/C)_χ are indicated.

In addition to efficiency and homogeneity, a third issue related to dispersal completeness becomes prominent, as a non-trivial proportion of particles is observed to dwell in the central area, as shown in figures 2(c) and 2(d). The incompleteness of dispersal is conspicuously the undesired outcome for most applications. Therefore, the extent of dispersal completeness should be regarded as one of the essential dispersal properties. The volume fraction of particles eventually residing inside the central region, χ, is a suitable measurement of the completeness. Estimating χ in systems where inbound particle jetting has not yet ceased during the computational times received special attention (a detailed description of the estimation method is provided in Appendix C). Figure 5 shows the variation in χ with M/C. The χ(M/C) curve exhibits a plateau below 0.1 as M/C increases from O(10⁰) to O(10²) until M/C reaches a threshold, (M/C)_χ ~ 250. A substantial increase in χ(M/C) ensues, accompanied by strong fluctuations between 0.3 and 0.6. Complete dispersal is expected for a minimum fraction of particles residing in the central area at the end of dispersal. We thus set the criterion for complete dispersal as follows:

(3.3)\begin{equation}\chi < 0.1.\end{equation}

Notably, the M/C thresholds corresponding to the three criteria are in the sequence of (M/C)_χ < (M/C)_κ < (M/C)_ξ. Specifically, (M/C)_ξ is one order larger than (M/C)_χ and M/C)_κ.

3.3. Classification of the dispersal mode

With the aid of the criteria for the efficiency, homogeneity and completeness of explosive dispersal, diverse dispersal behaviours are classified into four distinct modes, namely: (i) an ideal dispersal that satisfies all three criteria for efficiency, homogeneity and completeness ((3.1)–(3.3)); (ii) partial dispersal that satisfies the criteria for efficiency and homogeneity but fails to meet the criterion for completeness; (iii) retarded dispersal that only meets the criterion for efficiency; and (iv) failed dispersal that meets none of the three criteria. The classification of various dispersal processes in the parameter space of ξ, κ and χ is illustrated in figure 6, where the symbols representing the dispersal systems are rendered according to the respective M/C value. As the M/C increases from O(10⁰) to O(10¹), to O(10²), and finally to O(10³) and higher, the dispersal mode transitions from ideal to partial, retarded and finally to failed, highlighting the dominant role played by the M/C in determining the dispersal mode. In explosive dispersal experiments using bursters as explosion sources, the retarded and failed dispersal modes are rarely observed because the maximum (M/C)_exp in the reported experiments rarely reaches the order of O(10²) (Loiseau et al. 2018; Pontalier et al. 2018), which is equivalent to (M/C)_gas (gases at ambient temperature) of the order of O(10⁰) (see Appendix A) and most likely corresponds to the ideal or partial dispersal modes.

Figure 6. Mode categorization of various dispersal systems in the parameter space of ξ, κ and χ. Dispersal systems classified into different dispersal modes are denoted by symbols with distinct shapes. Symbols are coloured according to the M/C of the corresponding system.

We should be fully aware that the order of magnitude rather than the actual values of the M/C thresholds, (M/C)_i (i = ξ, κ, χ), indicated in figure 5, serves as the preliminary guide for the design or evaluation of the dispersal system. The phase map constructed in the parameter space shown in figure 6 provides a more reliable and refined framework. However, the calculation of ξ, κ and χ for a specific dispersal system requires knowledge of the whole dispersal process, which is often inaccessible or only accessible in hindsight. Therefore, the correlation between the structure of the dispersal system and the resulting dispersal mode must be established, which requires an in-depth understanding of the physics governing the transitions between different dispersal modes.

4. Analysis

4.1. Macroscale particle-flow coupling

The dynamics of the particle ring driven by the central pressurized gases are primarily governed by the complex coupling between particles and central flows. Figures 7(a)–7(d) present the space–time (R–t) diagrams of the pressure fields for the systems shown in figures 2(a)–2(d), where the circumferentially averaged gas pressure is used to plot the R–t diagrams. Similar to the dispersal behaviour of particle rings, the central flow field undergoes a two-stage evolution. The first stage, coinciding with the shock compaction phase of the ring, is dominated by the reciprocation of shock fronts between the centre and the inner surface, as evident in figure 7(a). Each incident shock invokes a CF traversing the particles, which is discernible in figure 2(e). The reciprocation of shock fronts becomes invisible after the significant expansion of the inner surface of the ring. Thereafter, the rapid expansion of the ring induces the overexpansion of the central gases, leading to a substantial decrease in pressure and eventually the reversal of the pressure gradient direction. Subject to the inwards directed pressure gradient forces, the expanding ring decelerates and tends to implode, which subsequently promotes the recovery of the pressure in the central gas pocket. Once the central pressure is increased above the ambient pressure, the outwards directed pressure gradient forces are restored across the thickness of the ring. The pinching ring subsequently gains outwards momentum, with a second expansion likely ensuing. This macroscale particle–gas coupling is adequately manifested by the synchronization between the pulsation of the particle ring and the fluctuation of the pressure inside the gas pocket, as plotted in figure 8(a), which share the same frequency with a 1/2 phase difference.

Figure 7. R–t diagrams of circumferentially averaged pressure in the dispersal system (a) 9.7-245-10-0.5, (b) 103.7-100-30-0.6, (c) 1024-20-50-0.6 and (d) 4875-20-140-0.6. The white and pink dashed lines indicate the IFP (OPF) and IS_CB (ES_CB), respectively.

Figure 8. (a) Fluctuations of the centre of mass of CB in terms of R_CB(t) and the pressure in the centre, P_g(R = 0). (b) Temporal variation in mass fraction of gases retained inside the gas pocket, χ_g, in the system 1024-20-50-0.6 and 4875-20-140-0.6. (c) Radial profiles of pressure across the thickness of ring at three sequent contraction-to-expansion transitions in the system 4875-20-140-0.6. The dotted lines in panels (a) and (b) indicate the overall temporal variations in the averaged R_CB, peak P_g and peak χ_g.

The pulsation of the particle ring resembles the pulsation of a gas bubble in the scenario of an underwater explosion (Wang et al. 2021). In contrast to the enduring bubble pulsation in the underwater explosion, which slowly loses momentum mainly due to the weak viscous dissipation inside the water, the particle ring only sustains a limited number of pulsation cycles as a result of quickly weakening particle–gas coupling. Multiple causes should be considered. The inelastic collision and friction between particles modestly dissipate the kinetic energy of the ring, acting similar to the viscosity of fluids. Nevertheless, the net gas flow-out associated with the gas infiltration plays a major role in diminishing the particle–flow coupling. As indicated in figure 8(b), only 66 % and 38 % (mass fraction) of gases are retained in the gas pockets at the first expansion–contraction transitions in systems 4875-20-140-0.6 and 1024-20-50-0.6, respectively. Notably, the gases flow into the gas pocket when the central pressure becomes negative (i.e. below ambient pressure). Thus, the χ_gas fluctuates in sync with the central pressure but decreases in the long term, indicating net gas loss over time. Consequently, the overall central pressure decreases with increasing damping fluctuation. The pressure difference between the central gas pocket and the ambient gas outside the CB is accordingly quickly normalized, as indicated by the flattening of the radial profiles of pressure across the thickness of the ring, as shown in figure 8(c), substantially reducing the primary driving forces of particles, namely, the pressure gradient forces and drag forces.

Two other important mechanisms are responsible for dissolving the particle–flow coupling, as illustrated in figure 9. As shown in figure 9(a), the initial impulse imparted to the particle ring is sufficiently strong to enable the expanding ring to quickly break away from the influence of the central flows. Likewise, the central gas pocket barely experiences the presence of the ring. This scenario occurs in systems with an M/C of the order of O(10⁰) (see figures 2a and 7a). However, the particle–gas coupling is sustained only if the ring, more precisely the densely packed CB, continues to provide sufficient confinement for the gas pocket. Once the CB disintegrates, the central gases thrust out, terminating particle–gas coupling (figure 9b).

Figure 9. Illustrations of two particle–flow decoupling mechanisms. (a) Fast expanding DB escapes the central flows with negative P_g. (b) DB prematurely disintegrates into the dilute particle cloud which cannot confine the central gases.

Different mechanisms dominate particle–gas coupling in different dispersal systems, leading to markedly varying degrees of particle–gas coupling. Here, we classify particle–gas coupling into decoupling, weak coupling, medium coupling and strong coupling regimes. In the decoupling regime, rings are only subjected to the initial impetus imparted by the central pressurized gases and thereafter evolve on their own (see figures 2a and 7a), which normally occurs in systems with an M/C of the order of O(10⁰). In the weak coupling regime, the innermost layers of particles accumulate in the centre due to the negative overpressure (below the ambient pressure) inside the gas pocket, while the bulk of particles are dispersed into the outer space. This phenomenon is illustrated in figures 2(b) and 7(b). As more particles are influenced by the central negative pressure and accumulate in the centre, the particle–gas coupling is classified as medium. In the medium coupling regime, the contraction of the particle ring is a solitary event since the densely packed CB disintegrates during ring contraction, as shown in figures 2(c) and 7(c), which often leads to abnormally high proportions of undispersed particles. In contrast, the strong particle–gas coupling results in the persistent pulsation of the ring, more strictly speaking, the densely packed CB, as shown in figures 2(d) and 7(d).

Obviously, in the decoupling and weak particle–gas coupling regimes, only the minimum fraction of particles accumulates in the centre. However, in the medium and strong particle–gas coupling regimes, a non-trivial proportion of particles is hurled into the centre during each contraction-to-expansion transition. The dispersal accordingly becomes incomplete. The most disastrous scenario occurs when the ring disintegrates during the contraction, and hence, the majority of particles shift inwards and collide with each other in the centre, leading to singularly high values of χ. The emergence of the worst scenario accounts for the significant fluctuation of the χ(M/C) curve beyond (M/C)_χ and the large variability in data among different systems with the same M/C values (see figure 5). However, the strong particle–flow coupling entails an inhomogeneous dispersal since the bulk of particles is generally located in an annular band where the pulsating CB eventually stops. A high concentration annular region is formed in this region.

4.2. Characterization of macroscale gas–flow coupling

The degree of particle–flow coupling depends on the relative importance of the characteristic times associated with the flow evolution inside the central gas pocket and the dynamics of the particle ring. The central flow evolution is represented by the temporal variation in the central pressure, which is a function of the initial state, the equation of state (EOS) of gases and the evolution path. We employ the time required by the overpressure in the centre to transition from positive (above) to negative (below the ambient pressure), t_pr, as the characteristic time of flow evolution. In § 3.2, we introduce two characteristic times, t_ring and t_dense, for the dynamics of the dispersed particle ring. For rings experiencing medium or strong particle–gas coupling, ring pulsation is the third defining event. Thus, we introduce a third characteristic time, t_ring,expcon, at which the first expansion-to-contraction transition occurs. The derivation of t_ring,expcon from the trajectory of the centre of mass of the CB is presented in Appendix C. Figures 10(a) and 10(b) plot the variations in t_pr and t_ring,expcon with the M/C. In contrast to the moderate increase in t_pr in the range of 1–6 ms throughout, t_ring,expcon remains infinity due to the absence of ring contraction until the M/C exceeds a critical value, (M/C)_cr ~ 350. At values greater than (M/C)_cr, t_ring,expcon plunges to the order of O(10¹) as M/C increases to the order of O(10³), and varies between 10 and 16 ms thereafter.

Figure 10. Variations in (a) t_pr and (b) t_ring,expcon with the increasing M/C. The dashed portion of the t_ring,expcon(M/C) curve is drawn only as a guide for the eye rather than represent the actual tend of t_ring,expcon since t_ring,expcon remains infinity in this range of M/C.

The degree of the gas–flow coupling can be quantitively evaluated by determining the relative importance of three pivotal events of ring dynamics with respect to the flow evolution, specifically the ratios between t_ring, t_ring,dense, t_ring,expcon to t_pr, which are denoted by $\varPi$, $\varOmega$ and $\varPsi$, respectively:

(4.1)\begin{gather}{\varPi } = \frac{{{t_{ring}}}}{{{t_{pr}}}},\end{gather}

(4.2)\begin{gather}{\varOmega } = \frac{{{t_{ring,dense}}}}{{{t_{pr}}}},\end{gather}

(4.3)\begin{gather}{\varPsi } = \frac{{{t_{ring,exp\text{-}con}}}}{{{t_{pr}}}}.\end{gather}

Figure 11 shows the variations in $\varPi$, $\varOmega$ and $\varPsi$ with increasing M/C. Both $\varPi$ and $\varOmega$ show a semiexponential dependence on M/C when M/C increases from O(10⁰) to O(10³). In contrast, $\varPsi$ plummets from infinity when M/C is less than (M/C)_cr to the order of O(10⁰), and plateaus between 3 and 2 as M/C increases to the order of O(10³). The convergence of $\varPsi$ combined with the infinite value of $\varOmega$ indicates that a consistent particle–flow coupling pattern emerges towards the upper limit of M/C.

Figure 11. M/C dependences of dimensionless parameters $\varPi$, $\varOmega$ and $\varPsi$. The critical M/C corresponding to the transitions between the decoupling, weak, medium and strong coupling are indicated by (M/C)$_\varPi$, (M/C)$_\varPsi$ and (M/C)$_\varOmega$, respectively.

The particle–gas decoupling regime requires the ring expansion to be faster than the evolution of the central pressure, leading to

(4.4)\begin{equation}{\varPi } < 1.\end{equation}

In the weak particle–flow coupling regime, ring expansion begins to be affected by the negative pressure in the overexpanded central gas pocket, while the effect of the flow is not sufficient to induce ring contraction as a whole. In this regime, t_ring,expcon is at least one order higher than t_pr. The criterion for weak particle–flow coupling is

(4.5a,b)\begin{equation}{\varPi } > 1,\quad {\varPsi } > 10.\end{equation}

The medium particle–flow coupling requires a single round of ring contraction, suggesting that the dense CB only survives the first expansion–contraction cycle. Thus, t_ring,expcon and t_ring,dense should be of the same order as t_pr, which is expressed as follows:

(4.6a,b)\begin{equation}{\varPsi } \sim O({10^0}),\quad {\varOmega } \sim O({10^0}).\end{equation}

The strong particle–flow coupling regime is embodied by an enduring dense CB whose survival time far exceeds the time scale of the flow evolution. The corresponding criterion is

(4.7)\begin{equation}{\varOmega } > 10.\end{equation}

For the group of dispersal systems studied here, the M/C thresholds corresponding to the transitions of the decoupling and weak, weak and medium, and medium and strong coupling regimes are (M/C)$_\varPi$ ~ 100, (M/C)_ψ ~ 350 and (M/C)$_\varOmega$ ~ 800, respectively. As argued in § 3.3, the M/C thresholds delimiting different particle–flow regimes derived from one group of dispersal systems should be applied with caution to another group with different granular materials or different explosion sources.

Employing the criteria listed in (4.4)–(4.7), the dispersal systems studied here are classified into different particle–flow regimes in the parameter space of $\varPi$, ψ and $\varOmega$, as shown in figure 12. As the M/C increases from O(10⁰) to O(10³) over three orders of magnitude, the particle–flow regime experiences the crossover from decoupling to weak, then to medium and finally to strong coupling.

Figure 12. Categorization of the particle–flow coupling regimes for various dispersal systems in the parameter space of $\varPi$, $\varOmega$ and $\varPsi$. Dispersal systems classified into different coupling regimes are denoted by symbols with distinct shapes. Symbols are coloured according to the M/C of the corresponding system.

4.3. Correlation between the dispersal mode and particle–flow coupling

The explicit correspondences between the dispersal modes and the particle–flow coupling regimes are illustrated in figure 13. The particle–gas decoupling regime guarantees an ideal dispersal mode, while the strong particle–flow coupling regime is necessary for a failed dispersal mode. The dispersal mode tends to be ideal towards the lower limit of the weak coupling regime. The retarded dispersal mode occurs in the upper limit of the medium coupling regime. The partial dispersal mode likely occurs towards either the upper limit of the weak coupling regime or the lower limit of the medium coupling regime.

Figure 13. Phase diagram of the dispersal modes and the particle–flow coupling regimes along the M/C axis, illustrating the correspondence between them. The phase boundaries are determined via the criterion for the respective dispersal mode or particle–flow coupling regime and specifically intended for the group of dispersal systems investigated in the present work.

The correspondence between the dispersal modes and the particle–flow coupling regimes provides a plausible approach to either predicting the dispersal mode for an available system or engineering the system for the desired dispersal mode. More specifically, if the correlations between the structure of the dispersal system and the dimensionless parameters $\varPi$, $\varOmega$ and $\varPsi$ are established, we are able to predict the particle–flow coupling regime and proceed to estimate the resulting dispersal mode. Conversely, the particle–flow coupling regime corresponding to a specified dispersal mode sets the range limits for $\varPi$, $\varOmega$ and $\varPsi$, which in turn place constraints on the structure of the dispersal system.

Since $\varPi$, $\varOmega$ and $\varPsi$ are ratios between t_ring, t_ring,dense, t_ring,expcon and t_pr, respectively, modelling these characteristic times as a function of a variety of structural parameters is among our priorities. The first and third parameters are associated with shock compaction and the first contraction of the ring, during which the bulk of particles likely remains densely packed such that pressure diffusion alongside gas infiltration is the dominant particle–gas coupling mechanism. The second characteristic time t_ring,dense involves particle shedding and multiple inbound and outbound jetting, whose mechanisms are far from well understood, as addressed in § 5. Hence, our main task in the next section is to predict t_pr, t_ring and t_ring,expcon by modelling the particle ring as a coherent granular medium rather than the collection of discrete grains.

4.4. Prediction of the particle–gas coupling regime

Predicting the t_pr, t_ring and t_ring,expcon for different dispersal systems requires modelling of the shock compaction phase and the subsequent pulsation of the particle ring, which involve distinct mechanisms and are accounted for separately. The final state of the shock-compacted ring sets the initial condition for the ring pulsation model.

The shock compaction and ring pulsation models both regard the particle ring as a compressible porous medium whose total mass is retained throughout. The transmitted shock through the particles is minimal and negligible due to the relatively high packing fraction (ϕ ₀ > 0.5). The driving forces of particles are the pressure gradient forces and drag forces across the thickness of the particle ring, F _∇P and F_drag (units N kg⁻¹), which are established by the diffusional pressure field (Britan & Ben-Dor 2006). Since the flow velocity relative to the particles depends on the local pressure gradient dictated by the Darcy and Forchheimer laws (Britan & Ben-Dor 2006), F_drag is proportionate to F _∇P, as shown in (4.8):

(4.8)\begin{equation}{F_{drag}} = \frac{{1 - {\phi _p}}}{{{\phi _p}}}{F_{\boldsymbol{\nabla }P}} ={-} \frac{{1 - {\phi _p}}}{{{\phi _p}}}\frac{{\boldsymbol{\nabla }{P_g}}}{{{\rho _p}}}.\end{equation}

Note that the direction of F _∇P is opposite to that of the pressure gradient ∇P_g, which may be directed either inwards or outwards depending on whether P_g is greater than or less than the ambient pressure P_amb. The deduction of (4.8) is presented in Appendix D. Meanwhile, the gases flow out of or flow into the central gas pocket via gas infiltration through particles. By ignoring the nonlinear term, the Darcy law prescribes the gas flow-out or flow-in rate as follows:

(4.9)\begin{equation}{\dot{m}_g} = 2{\rm \pi}{\rho _g}{R_{in}}\frac{k}{\mu }\boldsymbol{\nabla }{P_g},\end{equation}

where ρ_g and μ are the instantaneous gaseous density and viscosity inside the central gas pocket, respectively, and the permeability k is a function of the packing fraction ϕ_p described by the Ergun equation (Felice 1994):

(4.10)\begin{equation}k = \frac{1}{{150}}\frac{{{{(1 - \phi _p^{})}^3}}}{{\phi {{_p^2}^{}}}}d_p^2.\end{equation}

The value of ρ_g evolves with the volumetric variation of the central gas pocket and the cumulative mass flux as follows:

(4.11)\begin{equation}{\rho _g} = {\rho _{g,0}}\frac{{R_{in,0}^2}}{{R_{in}^2}} + 2\frac{k}{\mu }\int_0^t {\frac{{{\rho _g}}}{{{R_{in}}}}\boldsymbol{\nabla }{P_g}\,\textrm{d}t.}\end{equation}

The gases in the central gas pocket are assumed to undergo isothermal expansion:

(4.12)\begin{equation}{P_g} = {\rho _g}R{T_0}.\end{equation}

Once ρ_g is determined using (4.11), the central gaseous pressure P_g is estimated. We then calculate the gaseous temperature using the ideal gas EOS and the corresponding viscosity with the Sutherland model.

4.4.1. Shock compaction model

The shock compaction model aims to predict the kinetic energy imparted to the particle rings at the end of shock compaction. The schematic representation of the geometry considered in the model is shown in figure 14, where a CF propagates at a velocity of V_comp. The packing fraction jumps from ϕ ₀ to ϕ_comp (ϕ_comp = 0.68) across the CF, while the particles about to be compacted by the CF gain the velocity of u_comp(R_comp), where R_comp is the radius of the CF. In a cylindrical geometry, the mass conservation in the annular compacted band requires the particle velocity and acceleration, u_comp(R) and ${\dot{u}_{comp}}(R)$, to satisfy (4.13) and (4.14), respectively:

(4.13)\begin{gather}{u_{p,comp}}(R) = \frac{{{V_{in}}{R_{in}}}}{R},\end{gather}

(4.14)\begin{gather}{\dot{u}_{p,comp}}(R) = \frac{{{{\dot{V}}_{in}}{R_{in}}}}{R}\textrm{ + }\frac{{V_{in}^2}}{R} - \frac{{{{({V_{in}}{R_{in}})}^2}}}{{{R^3}}}\textrm{,}\end{gather}

where V_in is the velocity of the inner surface, also the particle velocity here, V_in = u_comp(R_in). At the CF,

(4.15)\begin{gather}{u_{p,comp}}({R_{\textrm{co}mp}}) = \frac{{{V_{in}}{R_{in}}}}{{{R_{comp}}}},\end{gather}

(4.16)\begin{gather}{\dot{u}_{p,comp}}({R_{comp}}) = \frac{{{{\dot{V}}_{in}}{R_{in}}}}{{{R_{comp}}}} + \frac{{V_{in}^2}}{{{R_{comp}}}} - \frac{{V_{in}^2R_{in}^2}}{{R_{comp}^3}}.\end{gather}

The V_comp and V_in should meet the Rankine–Hugoniot condition:

(4.17)\begin{equation}{V_{in}} = \frac{{{V_{comp}}{R_{comp}}}}{{{R_{in}}}}\left( {1 - \frac{{{\phi_0}}}{{{\phi_{comp}}}}} \right).\end{equation}

The momentum balance of the annular compacted band shown in figure 14 is calculated using (4.18):

(4.18) \begin{align}& {\rho _p}{\phi _{comp}}\int_{{R_{in}}(t)}^{{R_{comp}}(t)} {{{\dot{u}}_{p,comp}}(R)R\,\textrm{d}R} ={-} {\rho _p}{\phi _0}{u_{p,comp}}({R_{comp}}){V_{comp}}{R_{comp}}\nonumber\\ & \quad + {\rho _p}{\phi _{comp}}\int_{{R_{in}}(t)}^{{R_{comp}}(t)} {{F_{\boldsymbol{\nabla }P}}(R) \cdot R\,\textrm{d}R} + {\rho _p}{\phi _{comp}}\int_{{R_{in}}(t)}^{{R_{comp}}(t)} {{F_{drag}}(R) \cdot R\,\textrm{d}R} . \end{align}

Figure 14. Schematic representation of the wedge volumetric element with unit cross-sectional area taken into consideration in the (a) shock compaction and the (b) ring pulsation model. Inset in panel (b) shows the numerically derived compression curve and the fitting curve for the particle packings investigated in the present work.

The first term on the right-hand side of (4.18) arises from the growing mass of the compacted band. The second and third terms on the right-hand side of (4.18) represent the total pressure gradient force and the total drag force exerted on the compacted band with a cross-section of unit area. As a first-order estimation, we assume a linear pressure gradient across the thickness of the compacted band:

(4.19)\begin{equation}{F_{\boldsymbol{\nabla }P}} ={-} \frac{{\boldsymbol{\nabla }P}}{{{\rho _p}}} = \frac{{{P_g} - {P_{amb}}}}{{{\rho _p}({R_{comp}} - {R_{in}})}}.\end{equation}

Substituting equations (4.13)–(4.17) and (4.19) into (4.18) yields

(4.20)\begin{align}{\dot{V}_{in}}({R_{comp}} - {R_{in}}) + V_{in}^2\left[ {\frac{{{R_{comp}}}}{{{R_{in}}}} + \frac{{{\phi_{comp}}}}{{{\phi_{comp}} - {\phi_0}}}\frac{{{R_{in}}}}{{{R_{comp}}}} - 2} \right] = \frac{{{P_g} - {P_0}}}{{{\phi _{comp}}{\rho _p}}}\frac{1}{2}\left( {\frac{{R_{comp}^{}}}{{{R_{in}}}} + 1} \right).\end{align}

Equation (4.20) describes the evolution of ${\dot{V}_{in}}$ with the initial condition of ${\dot{V}_{in}} = 0$ and R_comp = R_in at t = 0. The integration of ${\dot{V}_{in}}$ results in V_in and V_comp (4.18), whose integrations in turn yield R_in and R_comp. Equations (4.10)–(4.20) constitute the complete formulations of the shock compaction model, which can be solved numerically as elaborated in Appendix E.

The trajectories of the inner surface and CF in the dispersal system 4875-20-140-0.6 are plotted in figure 15(a) against their counterparts derived from simulations, showing quite good agreement. We also compared the predicted velocities of the ring outer surfaces at the end of the shock compaction phase, V_out,pre, and the simulation results, V_out,num, in a variety of dispersal systems in figure 15(b). A good consistency is evident, supporting the reliability and accuracy of the shock compaction model in terms of predicting the shock compaction dynamics of particles in the radial geometry.

Figure 15. (a) Trajectories of the ring inner surface and the compaction front in the system 4875-20-140-0.6 predicted from the shock compaction model (R_in,pre, R_CF,pre) and derived from the simulations (R_in,num, R_CF,num). (b) Comparison of the velocities of ring outer surface at the end of the shock compaction phase predicted from the shock compaction model, V_out,pre (empty symbols), and derived from the simulations, V_out,num (solid symbols).

4.4.2. Ring pulsation model

After shock compaction, the particle ring undergoes the pulsation phase and itself is allowed to dilate or shrink, which entails varying ϕ_p capped by ϕ_comp. Figure 14(b) illustrates the force balance of a representative wedge element during the ring pulsation phase. Inside the wedge, a representative volume element ($\varLambda$) with inner and outer radii of R and R + δR is subjected to two body forces, F _∇P and F_drag, and two surface forces exerted by the particles in contact with the inner and outer surfaces of $\varLambda$, F_inner and F_outer. Both F_inner and F_outer arise from the granular pressure, P_gra, namely, F_inner = P_gra ⋅ R and F_outer = P_gra ⋅ (R + δR). The force balance of $\varLambda$ is described in (4.21):

(4.21)\begin{equation}{\rho _p}{\phi _p}R\delta R{\dot{V}_\varLambda } = {\rho _p}{\phi _p}{F_{\boldsymbol{\nabla }P}}R\delta R + {\rho _p}{\phi _p}{F_{drag}}R\delta R + {P_{gra}}R - {P_{gra}}(R + \delta R),\end{equation}

where ${\dot{V}_\varLambda }$ is the acceleration of $\varLambda$. After substituting equation (4.8) into (4.21), (4.21) reduces to

(4.22)\begin{equation}{\dot{V}_\varLambda } ={-} \frac{{\boldsymbol{\nabla }{P_g}}}{{{\rho _p}{\phi _p}}} - \frac{{{P_{gra}}}}{{{\rho _p}{\phi _p}R}}.\end{equation}

The mass conversation of the particle ring requires

(4.23)\begin{equation}{\phi _p} = \frac{{(R_{out,0}^2 - R_{in,0}^2){\phi _0}}}{{(R_{out}^2 - R_{in}^2)}}.\end{equation}

A reasonable assumption is a steady diffusional pressure field achieved across the thickness of the ring, which yields the pressure gradients at the inner and outer surface of the ring (Morrison 1970):

(4.24)\begin{gather}\boldsymbol{\nabla }{P_g}({R_{in}}) = \frac{{p_{amb}^2 - P_g^2}}{{2{P_g}({R_{out}} - {R_{in}})}},\end{gather}

(4.25)\begin{gather}\boldsymbol{\nabla }{P_g}({R_{out}}) = \frac{{p_{amb}^2 - P_g^2}}{{2{P_{amb}}({R_{out}} - {R_{in}})}}.\end{gather}

By substituting equations (4.23)–(4.25) into (4.22), we obtain the formulae for the accelerations of the innermost and outermost volume elements ($\varLambda$_in and $\varLambda$_out), ${\dot{V}_{{\varLambda _{in}}}}$ and ${\dot{V}_{{\varLambda _{out}}}}$:

(4.26)\begin{gather}{\dot{V}_{{\varLambda _{in}}}} = \frac{{({R_{out}} + {R_{in}})}}{{{\rho _p}(R_{out,0}^2 - R_{in,0}^2){\phi _0}}}\left[ {\frac{{P_g^2 - P_{amb}^2}}{{2{P_g}}} - \left( {\frac{{{R_{out}}}}{{{R_{in}}}} - 1} \right){P_{gra}}} \right],\end{gather}

(4.27)\begin{gather}{\dot{V}_{{\varLambda _{out}}}} = \frac{{({R_{out}} + {R_{in}})}}{{{\rho _p}(R_{out,0}^2 - R_{in,0}^2){\phi _0}}}\left[ {\frac{{P_g^2 - P_{amb}^2}}{{2{P_{amb}}}} - \left( {1 - \frac{{{R_{in}}}}{{{R_{out}}}}} \right){P_{gra}}} \right].\end{gather}

The granular pressure, P_gra, which reflects the compression resistance of the granular medium, arises from the energy stored in elastic–plastic layers around each grain contact point. The value of P_gra as a function of the solid volume fraction α can be determined by differentiating the configuration energy B(α), as shown in (4.28) and (4.29) (Richard et al. 2010):

(4.28) \begin{gather}B(\alpha ) = \left\{ {\begin{array}{@{}ll} {a{{\left[\begin{array}{@{}l@{}} (1 - \alpha )\log (1 - \alpha ) + (1 + \log (1 - {\alpha_0}))(\alpha - {\alpha_0})\\ - (1 - {\alpha_0})\log (1 - {\alpha_0}) \end{array} \right]}^n}}&{\textrm{if}\;{\alpha_0} < \alpha < 1}\\ 0&{\textrm{otherwise}} \end{array}} \right.,\end{gather}

(4.29) \begin{gather}{P_{gra}} = \left\{ {\begin{array}{@{}ll} {\alpha {\rho_p}\dfrac{{\textrm{d}B(\alpha )}}{{\textrm{d}\alpha }}\textrm{ = } - an\alpha {\rho_p}(\log (1 - \alpha ) + 1){{\left( {\dfrac{{B(\alpha )}}{a}} \right)}^{(n - 1)/n}}}&{\textrm{if}\;{\alpha_0} < \alpha < 1}\\ 0&{\textrm{otherwise}} \end{array}} \right..\end{gather}

The positive parameter α ₀ in (4.28) and (4.29) corresponds to the solid volume fraction when the granular pressure is zero. Parameters a and n are also characteristic of a particular powder and, more precisely, of its response during quasistatic loading. By fitting the numerically derived compression curve of glass beads, as shown in the inset of figure 14(b), in the solid volume fraction range of 0.6–0.68, the fitting parameters are determined to be α ₀ = 0.61, a = 500 and n = 1.004. Notably, the variations in α consist of the contributions from both ϕ_p and ρ_p since α = (ϕ_p ρ_p)/ρ_{p ,0}, where ρ_{p ,0} is the initial material density of grains. At the pressure investigated in this paper, the material density of grains barely changes, ρ_p = ρ_{p ,0}. Hence, α in (4.28) and (4.29) is equivalent to ϕ_p.

Equations (4.10)–(4.12) and (4.26)–(4.29) constitute the complete formulations of the ring pulsation model, which is solved numerically using the method presented in Appendix E. Figures 16(a)–16(c) compare the predicted and simulation-derived trajectories of the centre of mass of CB (R_CB,pre and R_CB,num), central pressure (P_g,pre and P_g,num) and mass fraction retained in the gas pocket (χ_g,pre and χ_g,num) for system 4875-20-140-0.6, respectively. The predicted R_CB(t), P_g(t) and χ_g(t) curves all exhibit fluctuating characteristics and long-term convergence resembling those observed in the simulations (see figures 8a and 8b). Although the fluctuations in the amplitudes of the R_CB(t) and P_g(t) curves are overestimated, the predicted fluctuation periods are consistent with those derived from the simulations, which lends credence to the capacity of the ring pulsation model to predict t_pr and t_ring,expcon.

Figure 16. Comparisons of the temporal variations of (a) the centre of mass, R_CB, (b) the central pressure, P_g (R = 0), and (c) the mass ratio retained inside the gas pocket, χ_g, predicted from the ring pulsation model (dashed lines) and results derived from simulations (solid lines) in system 4875-20-140-0.6. Note that the final states of ring and central gases predicted by the shock compaction model serves as the initial state for the ring pulsation model.

Figures 17(a)–17(c) show the prediction errors of t_ring, t_pr and t_ring,expcon for dispersal systems with M/C spanning three orders of magnitude. Here, the prediction error ε_i (i = t_ring, t_pr and t_ring,expcon) is defined as ${\varepsilon _i} = |t_i^{pre} - t_i^{sim}|/t_i^{sim}$, where $t_i^{pre}$ and $t_i^{sim}$ represent the theoretical predictions and numerical results averaged over systems with the same M/C but varied combinations of structural parameters, respectively. As shown in figure 17(a), the prediction errors of t_ring vary within 30 %. The prediction errors of t_pr decrease from in excess of 40 % to less than 20 % as the M/C increases from O(10¹) to O(10³) (see figure 17b). The predicted t_ring,expcon also shows good agreement with the simulation results when the M/C is greater than (M/C)$_\varPsi$, as indicated in figure 17(c). However, due to the lack of a particle shedding mechanism, the ring pulsation model fails to account for t_ring,expcon approaching infinity when M/C is less than (M/C)$_\varPsi$, as observed in the simulations.

Figure 17. Prediction errors of (a) t_ring, (b) t_pr and (c) t_ring,expcon for systems with M/C spanning over four orders of magnitude. The variation in numerically derived t_ring,expcon with the M/C is also coplotted in panel (c).

Figures 18(a) and 18(b) plot the variations in the predicted $\varPi$ and $\varPsi$, respectively, with the M/C over three orders of magnitude. Although the threshold (M/C)$_\varPsi$ signifying the transition from the weak to medium particle–flow coupling regimes cannot be properly identified since the predicted $\varPsi$ (M/C) curve remains less than 10 throughout, we can readily determine the threshold (M/C)$_\varPi$ from the predicted $\varPi$ (M/C) curve in which the decoupling regime is differentiated from other regimes. The predicted $(M/C)_\varPi ^{pre}$ (~280) is quite close to that deduced from simulations, $(M/C)_\varPi ^{sim}$~100. Therefore, the combination of the shock compaction and ring pulsation models shows a high prediction accuracy for the decoupling regime. If a particular dispersal system is predicted to be within the particle–flow decoupling regime, an ideal dispersal mode is guaranteed.

Figure 18. Comparisons of variations in theoretically predicted and numerically derived (a) $\varPi$ and (b) $\varPsi$ with increasing M/C.

Although the shock compaction and ring pulsation models do not explicitly incorporate the M/C, the range of structural parameters, including P ₀, ρ_{g ,0}, R_{in ,0}, R_{out ,0}, ϕ ₀ and ρ_p, included in the models are integral parts of the M/C, as calculated using (2.1). Hence, the shock compaction and ring pulsation models reveal the underlying mechanisms by which the M/C influences the imparted energy during the shock compaction phase and the energy transfer between the centre flows and particle ring during the subsequent ring pulsation phase, respectively.

5. Discussion

The inability of the combined shock compaction and ring pulsation models to predict dispersal modes other than the ideal dispersal mode arises from their failure to account for particle stripping from the bulk. The stray particles that lose contact with neighbouring particles do not provide sufficient confinement to the central gases or experience pressure gradient forces. Macroscale particle–flow coupling is not valid for stray particles. The stray particles are thereafter referred to as loosened particles. The particle loosening mechanism and rate critically affect the extent and sustainability of the macroscopic particle–flow coupling, thereby altering the fundamental properties of the dispersal process. Here, we focus on addressing the origins and initiations of different particle loosening events. The respective modelling process should be a future pursuit.

As evidenced in the R–t diagrams of ϕ_p in different dispersal systems (figures 2a–2d), three primary particle loosening events are initiated at different phases of ring dynamics: (i) particle shedding initiated from the internal surface immediately after shock compaction; (ii) outbound particle jetting initiated from the external surface upon the onset of the first ring contraction; and (iii) inbound (outbound) particle jetting initiated from the internal (external) surface upon contraction-to-expansion (expansion-to-contraction) transitions after the first ring contraction. Figures 19(a)–19(c) show the morphological evolutions of the dispersed ring dominated by different particle loosening events. The first event generates a diffusive dilute particle layer textured by inwards protrusions (figure 19a). In contrast, the second event produces a number of well-defined particle jets with sharp cusps (figure 19b). Despite resembling shapes with the jets formed in the second event, the finger-like jets arising from the third event, either inbound or outbound, are much coarser, with blunt heads (figure 19c). These three particle-loosening events are referred to as internal particle shedding, external particle jetting and bidirectional particle fingering, respectively.

Figure 19. Morphological evolutions of the dispersed rings dominated by (a) the internal particle shedding, (b) the external particle jetting and (c) the bidirectional particle fingering. Panels (a), (b) and (c) are extracted from systems 103.7-200-50-0.6, 1024-20-50-0.6 and 4875-20-140-0.6, respectively.

The internal particle shedding illustrated in figure 19(a) dominates the dispersal process in the ideal dispersal mode. Figure 20(a) plots the growth of internal particle protrusions, l_shedding, defined by the average perturbation amplitude of the multimode internal surface, and the volume fraction of particles constituting the protrusions, ϕ_shedding, in system 103.7-200-50-0.6. The l_shedding and associated ϕ_shedding both remain minimal until the first shock impingement after the shock compaction, t_imp ~ 1 ms, highlighting the important role played by the shock impingement in the initiation of the internal particle shedding. Figure 20(b) depicts the pressure fields at the instant immediately after t_imp. The pressure inside the grooves of the perturbed inner surface is significantly increased due to the shock focusing at this site. If the ring is divided into a series of wedge elements, the wedge with the concave inner surface is pushed by the stronger pressure gradient forces and consequently moves faster than that with the bulged inner surface, as illustrated in figure 20(c). As a result, the minute cusps of the rough inner surface evolve into detectable protrusions. Afterwards, the protrusions consisting of loosened particles lag behind the bulk of the ring due to the lack of efficient pressure gradient forces. Meanwhile, more particles migrate from the bulk to the roots of the protrusion. In this manner, the internal particle shedding initiated from the inner surface progresses towards the outer surface. A complete account of shock-induced particle shedding is presented in our previous article (Li et al. 2022).

Figure 20. (a) Growth of internal particle protrusions, l_shedding, and the temporal variation in the volume fraction of particles constituting the protrusions, ϕ_shedding. Top left inset in panel (a) is the R–t diagram of P_g with the indication of t_imp; bottom right inset in panel (a) illustrates the definition of l_shedding. (b) Pressure contour at the instant just after t_imp superimposed by isobaric lines (white lines) and the disturbed inner surface (black line). (c) Illustration of the internal particle shedding mechanism initiated by the shock focusing upon the grooves of the rough inner surface. All plots in panels (a) and (b) are derived from the system 103.7-200-50-0.6.

The external particle jetting illustrated in figure 19(b) is widely observed in all dispersal modes except the ideal dispersal mode. Figure 21(a) plots the growth of the external jet length, l_jet,ex, defined by the average perturbation amplitude of the multimode external surface, and the temporal variation of the external surface velocity, V_out, in system 1024-20-50-0.6. The onset of destabilization of the external surface is signified by the substantial increase in l_jet,ex that coincides with the beginning of the deceleration of the external surface. This observation seems to suggest a Rayleigh–Taylor instability (RTI), which occurs as a heavy fluid decelerates into a light fluid (Kun et al. 2018). However, debate persists regarding whether the shock-loaded granular medium can be treated as a one-phase material on the bulk scale since the coupling between the solid grains and interstitial fluids is far from complete and in equilibrium (Han, Xue & Bai 2021). Previous studies have found that the generally accepted viscoplastic constitutive equations are inadequate to describe the behaviour of granular media in the wake of shocks on the bulk scale (Han et al. 2021). Hence, the RTI does not provide the definite answer to the origin of the external jetting.

Figure 21. (a) Growth of external jet length, l_jet,ex (left y axis), and the temporal variation of the external surface velocity, V_out (right y axis). Inset in panel (a) illustrates the derivation of l_jet,ex. (b) Granular velocity fields, u_p(R,θ), during the shock compaction (1.04 ms), the rarefaction propagation (1.28 ms) and the fully developed external jetting (14.76 ms). (c) Circumferential variations of u_p(θ) in the outmost layer at different times. Plots in panels (a)–(c) are derived from results in the system 1024-20-50-0.6.

Here, we attempt to provide insights into the particle-scale physics governing external jetting. Figure 21(b) shows the evolution of the granular velocity fields, u_p(R,θ), during shock compaction (1.04 ms), rarefaction propagation (1.28 ms) and the development of external jets (14.76 ms) in system 1024-20-50-0.6. The most striking finding is the presence of profuse localized high-velocity clusters. These clusters, which already emerge during shock compaction, albeit diffusive and inconspicuous, are significantly increased in strength and size during the rarefaction propagation phase. The circumferential variations of u_p in the outer layers with a thickness of 100d_p at different times are plotted in figure 21(c) to quantify this trend. During the propagation of the rarefaction wave, a prominent regularly fluctuating pattern of u_p(θ) emerges from the scattering distribution of u_p(θ) formed at the end of shock compaction. The fast-moving clusters in the outermost layers thrust through the external surface, forming jet-like protrusions, as suggested at the top of figure 21(b).

The heterogeneities exhibited in the velocity fields (figure 21b) may be attributed to the heterogenous momentum transmission channelled by the force chain network intrinsic to the granular media (Majmudar & Behringer 2005; Clark et al. 2015; Tadanaga et al. 2018). Extensive studies have revealed the highly heterogeneous strain/flow fields and intense granular temperature arising from complex force chain structures (Clark et al. 2015; Huang et al. 2016; Tadanaga et al. 2018; Xue et al. 2018; Mo et al. 2019). However, to our knowledge, no investigation has been reported concerning the augmentation and coarsening of heterogeneities with the aid of rarefaction relaxation, which closely correlates with the wavelength selection mechanism of the resulting jets and is worth further investigation.

Bidirectional fingering normally occurs in retarded and failed dispersal modes. Examining the R–t diagram of ϕ_p in system 4875-20-140-0.6 (figure 2d) and the morphological evolution of the corresponding ring (figure 19c), we ascertain two striking characteristics of the bidirectional fingering: (i) at each expansion-to-contraction (contraction-to-expansion) transition of the ring, a batch of finger-like protrusions are catapulted from the external (internal) surface; and (ii) a new generation of fingers always bud between the neighbouring fingers from the previous generations. In contrast to the filamentary external jets, the bidirectional fingers adopt a mushroom-like shape featuring a flared head and an elongated stem.

Figures 22(a) and 22(b) present the pressure field and particle velocity field in system 1024-20-50-0.6 before the second expansion-to-contraction transition of the ring, respectively. As shown in figure 22(a), fully developed external jets arising from the prior external jetting create a highly undulated external surface and grooves form between adjacent jets with bulged bottoms. For a particle ring with non-uniform thickness, the pressure gradient becomes non-uniform circumferentially. As illustrated in figure 23, a stronger pressure gradient develops across the necks of the ring aligned with the grooves. Consequently, the grooves gain higher velocities than the jet bottoms. As the inwards directed pressure gradient forces occur, the bulged bottoms of jets reverse the moving direction first while the grooves continue to expand outwards, as suggested in figure 22(b). A large chunk of particles eventually is catapulted from the groove, leaving a finger-like protrusion in the wake, as shown in figure 22(c). As the newly born fingers grow upwards, new grooves in between these fingers are generated, precipitating the formation of the next generation of fingers. The same mechanisms also account for the inbound fingering occurring at the contraction-to-expansion transitions. Figure 23 illustrates the mechanism governing the initiation of bidirectional fingering.

Figure 22. Pressure contour in the system 1024-20-50-0.6 at t = 18 ms superimposed by particles (white dots), the internal and external surfaces of CB (white dashed lines). (b,c) Mophologies of dispersal ring at (b) t = 22.4 and (c) 26 ms wherein the particles are coloured according to the instantaneous particle velocities.

Figure 23. Schematic of the initiation of bidirectional fingering due to the circumferential non-uniform pressure gradients.

An understanding of the physics governing the initiation and evolution of various particle loosening events provides an opportunity to properly model the particle loosening rate and the dynamics of loosened particles, which should consider particle–flow coupling on the microscale and the macroscale. For internal particle shedding, microscale particle–flow coupling refers to the interaction between reciprocating shocks and particles. The strength of the first incident shock after shock compaction determines the internal particle shedding rate. Due to the strongly decayed incident shocks after the shock compaction phase, the internal particle shedding becomes trivial for systems with an M/C of the order of O(10³). For external jetting, the number of jets is associated with the heterogeneity of the granular flows, which arises from particle-scale momentum transmission. Regarding bidirectional fingering, the proportion of ejected particles correlates with the circumferential variabilities of pressure gradients resulting from the undulating inner and outer surfaces. Therefore, the complex dispersal structures are the products of particle–flow coupling relations on both the macroscale and microscale.

6. Conclusions

In the present study, we reveal a remarkable variety of dispersal behaviours and hierarchical structures in explosive dispersed particle rings. Diverse dispersal processes are classified into four dispersal modes in terms of three fundamental attributes. Transitions between distinct dispersal modes are accounted for based on the prevailing macroscale particle–flow coupling. The phase diagrams of the dispersal mode and particle–flow coupling have been mapped in the space of the mass ratio, which is selected as one of the foremost important parameters. The combination of the continuum-based shock compaction and ring pulsation models is capable of predicting the dynamics of the ring and the central flow evolution at early time points, thereby successfully identifying the ideal dispersal mode. In addition, the theoretical models reflect the complex effects of the mass ratio, which incorporates the contributions of various structural parameters. By properly considering the particle loosening mechanisms, we expect that the modified models will allow the predictions of all dispersal modes, paving the way for the design and optimization of dispersal systems for specific engineering demands.