2.1. Governing equations and numerical algorithm

Numerical simulations were performed based on the 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, Ozel & Kolehmainen Reference Sundaresan, Ozel and Kolehmainen2018; Tian et al. Reference Tian, Zeng, Meng, Chen, Guo and Xue2020). 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 physical particles that a computational parcel represents is quantified using a scaling factor, namely 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 reported χ in previous literature ranges from O(10¹) to O(10³) (Osnes et al. Reference Osnes, Vartdal and Pettersson Reif2018; Koneru et al. Reference Koneru, Rollin, Durant, Ouellet and Balachandar2020). In the present work, χ is of the order of O(10¹).

For the gas phase, the volume-averaged governing equations [(2.1)–(2.3)] constructed in the Eulerian frame are based on a five-equation transport model, i.e. a simplified form of the Baer–Nunziato model, which has been modified to account for compressible multiphase flows ranging from dilute to dense gas–particle flows (Baer & Nunziato Reference Baer and Nunziato1986)

(2.1) \begin{gather} \dfrac{\partial(\phi_{f}\langle {\rho_f}\rangle)}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }({\phi _f}\langle {\rho_f}\rangle {{\tilde{\boldsymbol{u}}}_f}) = 0, \end{gather}

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

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

The gas volume fraction, velocity, density, pressure and the total energy of the gas are represented by ϕ_f, u_f, ρ_f, P_f and E_f, E_f = ρ_f e_f + 0.5 ρ_f u_f u_f, respectively. In (2.1)–(2.3), $\langle \;\rangle $ and $\,\,\widetilde {}\,\,$ denote phase-averaged and mass-averaged variables, respectively, ρ_p.i and u_{p ,i} are the density and velocity of parcel i, D_p,i is the drag force coefficient of parcel i and ϕ_p.i = w_i,f V_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, V_p,i and V_f are the volumes of parcel i and the fluid cell).

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

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

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

(2.6)\begin{gather}{f_{base}} = \left\{ {\begin{array}{*{20}{@{}cl}} {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}

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

In (2.4)–(2.5), C_d is the dimensionless coefficient of the drag force, sg is the specific weight of individual particles and r_p is the particle radius. For dense particle flows (ϕ_f < 0.8), (2.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 (2.6) and (2.7).

The particle phase is represented by discrete parcels whose motion is governed by Newton's second law ((2.8) and (2.9))

(2.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 {{F_{C,ij}}} ,\end{gather}

(2.9)\begin{gather}\frac{{\textrm{d}\kern 0.06em {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 (2.8). A soft sphere model, 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

(2.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, δ_n and u_n,ij are the overlap and normal velocity difference between parcels in contact; γ_n,p is a function of the parcel restitution coefficient ε_p

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

To solve the equations governing the gases, the weighted essentially non-oscillatory scheme was used to reconstruct the primary flow variables. 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 to the time integration. The equations describing the parcel velocity and position were discretized by 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 Reference Liu, Osher and Chan1994). Numerical details with regard to CCFD-DPM can be found in our previous works (Tian et al. Reference Tian, Zeng, Meng, Chen, Guo and Xue2020).

The CCFD-DPM framework has been validated against several benchmark experiments involving shock driven particle laden flows (Tian et al. Reference Tian, Zeng, Meng, Chen, Guo and Xue2020), such as Rogue's experiments (Rogue et al. Reference Rogue, Rodriguez, Haas and Saurel1998) which investigate the shock dissipation through particle curtains, the experiments carried by Britan & Ben-Dor (Reference Britan and Ben-Dor2006), which access both the gaseous and solid pressures inside particle columns impinged head on by shocks and the experiments of shock induced interfacial instability of granular media (Xue et al. Reference Xue, Du, Shi, Gan and Bai2018). Specifically the CCFD-DPM-based simulations reproduce the characteristic expanding velocity and the signature finger-like jetting structure of the explosive dispersal, further validating the applicability of the CCFD-DPM method in investigating the explosive dispersal phenomena (Tian et al. Reference Tian, Zeng, Meng, Chen, Guo and Xue2020).

2.2. Numerical set-up

The two-dimensional numerical configuration is shown in figure 1(a) which has a high pressure gas pocket with the initial pressure P ₀ and the radius of R_gas. The gas pocket is enclosed by a densely packed particle ring with inner and outer radii of R_{in ,0} and R_{out ,0}, respectively. The volume fraction of particle ring ϕ ₀ ranges from random loose packing ϕ _0,min = 0.5 to the random dense packing ϕ _0,max = 0.65. The annular particle ring domain is filled by computational parcels generated by the radius expansion algorithm (Yan, Hai-Sui & Glenn Reference Yan, Hai-Sui and Glenn2009). A population of parcels with artificially small radii that ensure no particle overlap is randomly created within the specified domain. Then, all parcels are expanded until the specified parcel size distribution and desired volume fraction are both satisfied. The real particle has a diameter of d_p while the diameter of the parcel uniformly ranges from 0.75d_p to 1.25d_p to avoid potential crystallization during shock compaction. A random but homogenous arrangement of parcels is achieved as shown in the zoomed-in inset of figure 1(a) wherein the parcels are coloured by the local Voronoi volume fraction ϕ_p,voro. The parcel has a density of ρ_p same as the real particles.

Figure 1. (a) Two-dimensional numerical configuration of the explosive dispersal system (only one quarter of the configuration is shown). Distribution of the simulated dispersal systems in the parameter space (M/C, E_gas) (b) and (d_p, ρ_p) (c). The symbol size in (b) and (c) scales with the recurrence frequency of the corresponding parameter combination in all numerical cases.

For the simplified configuration shown in figure 1(a), the macroscale set of structural parameters includes the geometries of the central gas pocket and the particle ring, R_gas, R_{in ,0} and R_{out ,0}, the gas pressure and temperature, P ₀ and T ₀, the total mass of particle ring, ${m_{ring}} = {\rm \pi}{\rho _p}{\phi _0}(R_{out,0}^2 - R_{in,0}^2)$, etc. Both R_gas and P ₀ contribute to the explosion energy of the gas pocket which can be approximated by the Brode equation (Crowl Reference Crowl2003)

(2.12)\begin{equation}{E_{gas}} = \frac{{({P_0} - {P_{amb}}){V_{gas}}}}{{\gamma - 1}} = \frac{{{\rm \pi} ({P_0} - {P_{amb}})R_{gas}^2}}{{\gamma - 1}},\end{equation}

where V_gas is the volume of the gas pocket, P_amb is the ambient pressure and γ is the ratio of specific heat. The mass ratio based on the masses of pressurized gases and the particle ring is given by

(2.13)\begin{equation}\textrm{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}

Since the particles in the CCFD-DPM method are modelled as smooth spheres, the surface roughness, the angularity and the material properties are difficult to account for. The microscale structure of the dispersal system is mainly embodied by the particle diameter and density, d_p and ρ_p, both contributing to the inertia of the particle. Figure 1(b,c) shows the distributions of more than 140 numerical cases in the macro- and micro-scale parametric spaces, namely (E_gas, M/C) and (d_p, ρ_p) spaces, respectively. The symbol size in figure 1(b,c) scales with the recurrence frequency of the corresponding parameter combination in all numerical cases. Via varying the macroscale structural parameters, the M/C of numerical cases ranges from O(10⁰) to O(10⁴), spanning over five orders of magnitude. On the other hand, d_p ranges from O(10¹) to O(10²) μm while ρ_p ranges from O(10¹) to O(10³). The vast variations in the macro- and microscale structural parameters fully expose the variability of the resulting dispersal behaviours which are normally prohibited by the experimental means. Exact values of the structural parameters in each numerical case are presented in table 1 in Appendix A. For clarity, the system is labelled by three symbols, M/C, −d_p (in units of μm), − ρ_p (in units of kg m⁻³).

3.1. Macroscale disposal behaviour

In order to gain knowledge of the overall explosive dispersal behaviour of the particle ring, a space–time (r–t) diagram manifesting the spatio-temporal variations in particle volume fraction ϕ_p(r, t) was constructed using the circumferentially averaged ϕ_p. Figure 2(a–f) displays six typical r–t diagrams for systems with different combinations of macro- and microscale parameters. The inner and outer boundaries of the dispersed particle cloud denoted by the dotted lines are defined in such way that 90 % of particles are enclosed by the boundaries. The corresponding snapshots showing the positions of the dispersed particles at very late times are presented in figure 3(a–f).

Figure 2. The r–t diagrams of ϕ_p(r, t) for systems 9.7-100-2500 (a), 104-100-1503 (b), 159-100-2500 (c), 288-60-1503 (d), 1035-100-2500 (e), 2043-100-2500 (f). The inner and outer boundaries of particle clouds are denoted by dotted lines.

Figure 3. Positions of the dispersed particles at very late times for systems 9.7-100-2500 (a), 104-100-1503 (b), 159-100-2500 (c), 288-60-1503 (d), 1035-100-2500 (e), 2043-100-2500 (f). Particles are rendered according to the instantaneous velocities. Note that the colour ranging between yellow and dark red represents outbound velocities while the colour ranging between cyan and navy blue represents inbound velocities.

Figure 2(a) shows the most ideal dispersal which occurs at the lower limit of the M/C range. The particle ring undergoes a persistent expansion accompanied by substantial particle shedding from the inner surface, forming an annular particle cloud laden with uniformly distributed particles (figure 3a). The particle cloud becomes increasingly dilute as the volume swells until the momentum of particles is depleted by the drag forces. If the ring expansion is not fast enough, the innermost layers are susceptible to the negative central pressure due to the overexpansion of the central gases and eventually sucked into the centre, as seen in figure 2(b,c) (also see figure 3b,c). Hereafter, a negative pressure refers to a pressure below the ambient one, otherwise the pressure is positive.

As seen in figure 2(d–f), with increasing M/C, the persistent expansion of the ring is gradually taken over by the expansion–implosion pulsation which is compromised or even terminated by the intense particle shedding from both inner and outer surfaces of the ring. A variety of dispersal behaviours contrasting with the ideal dispersal emerge. For the system 288-60-1503, as shown in figure 2(d), the particle ring disintegrates during the first implosion, hurling a majority of particles towards the centre. These particles collide with each other and lose their momentum due to the frequent inelastic collisions. Eventually, most particles reside in the central region instead of being propelled out, resulting in an incomplete dispersal (figure 3d). Another unwanted dispersal outcome occurs when the ring disintegrates at the very moment the ring is about to implode (figure 2e) or the pulsating ring comes to rest in an equilibrium location (figure 2f). In these scenarios, a significant proportion of particles end up in an annular region near the initial location of the ring (see figure 3e,f) with finger-like jets bursting from the inner and outer surfaces, carrying away a fraction of particles. Hence, figure 2(e,f) demonstrates a highly incomplete and non-uniform dispersal.

Most engineering applications desire a uniformly distributed particle cloud via a complete explosive dispersal. The completeness of the explosive dispersal can be quantified by the proportion of particles which are eventually dispersed out, χ. The value of χ is derived by running the numerical experiments until the ring disintegrates or ceases to pulsate and finally tallying the outbound particles. The uniformity of the dispersed particle cloud is assessed by measuring the deviation of the mass centre radius of the actual particle cloud, R_mass, from that of the hypothetically homodispersed particle cloud, R_{mass ,homo}, $\kappa = |{{R_{mass}} - {R_{mass,homo}}} |/{R_{mass,homo}}$. The values of R_mass and R_{mass ,homo} are calculated when the average particle volume fraction inside the region delimited by the inner and outer boundaries falls to 0.1. The mathematical definitions of χ, R_mass, R_{mass ,homo} and corresponding numerical procedures of equations are elaborated in Appendix B.

Figure 4(a,b) shows the variations in χ and κ with increasing M/C, respectively. The size of symbols scales with 1/d_p, while the colour of the symbols is rendered according to ρ_p. The detectable dependences of χ and κ on M/C emerge from figure 4(a,b). As M/C increases from O(10²) to O(10⁴), χ decreases significantly while κ undergoes a substantial and non-convergent increase, indicating that the dispersal becomes increasingly incomplete and non-uniform beyond a M/C threshold of around O(10²). It is worth noting that there exist significant variabilities in χ and κ among dispersal systems in the M/C range of O(10²) to O(10³) which have identical M/C but different ρ_p and d_p. Specifically, the system 55-100-3506 yields an ideal dispersal with a close-to-unity χ and a minimum κ, χ = 0.99 and κ = 0.003. By contrast, the systems with larger but much lighter particles, systems 55-450-417 and 55-250-210, result in the incomplete and non-uniform dispersals with χ = 0.8 (system 55-450-417) and 0.3 (system 55-250-210), and κ = 0.06 (system 55-450-417) and 0.16 (system 55-250-210). The non-trivial influences on the dispersal behaviour brought by the microscale parameters will be accounted for in §§ 3.3 and 3.4.

Figure 4. Variations in χ and $\kappa$ with increasing M/C. The symbol size scales with 1/d_p. The colour of symbols is rendered according to ρ_p. The master curves with 80 % confidence intervals are superimposed in the plots. Symbols denoted as 1, 2, 3 correspond to systems 55-100-3506, 55-450-417, 55-250-210, respectively.

3.2. Macroscale gas–particle shell coupling

The overall explosive dispersal behaviour characterized by χ and κ is predominately dictated by the macroscale gas–particle coupling, as revealed in the previous study (Kun et al. Reference Kun, Lvlan, Jiarui, Chunhua and Baolin2023). Figure 5(a–f) presents the space–time (r–t) diagrams of the gaseous pressure fields, P_gas(r, t), for the systems shown in figure 2(a–f), wherein the circumferentially averaged P_gas is used to plot the r–t diagrams. We superimpose the inner and outer boundaries of the particle cloud (denoted by white dotted lines) as well as those of the dense core band (denoted by light yellow dashed lines) in figure 5. The dense core band remains densely packed, ϕ_dense ≥ 0.3, so that it suffices to effectively confine the central gases and maintain the pressure differential between the inner and outer boundaries. The disappearance of the dense core band signifies the complete disintegration of the particle ring.

Figure 5. The r–t diagrams of P_gas(r, t) for systems 9.7-100-2500 (a), 104-100-1503 (b), 159-100-2500 (c), 288-60-1503 (d), 1035-100-2500 (e), 2043-100-2500 (f). The inner and outer boundaries of particle clouds and dense core band are denoted by while dotted lines and light yellow dashed lines.

There are two characteristic events defining the evolution of the central pressure field enclosed by the inner surface of ring. The first is the multiple reflections of the shock waves and rarefaction fans between the centre and the inner surface, which is most prominent in early times, as shown in figure 5(a). All waves quickly attenuate as the ring rapidly expands. The overexpansion of the central gases takes over to dominate the gas dynamics, causing a drastic decline of the central pressure which eventually becomes negative. The resulting adverse pressure gradient, which is directed outwards, exerts inward pressure gradient forces on the particles. Thereby, the expanding ring decelerates. If the expanding ring completely disintegrates, particles in the inner layers may well reverse their motion and be sucked into the centre. The central gases gush out and quickly mix with the outside gases. In this scenario, a uniform pressure field is established (figure 5b,c). Otherwise, the dense core band eventually begins to implode, which in turn gives rise to the recovery of the central pressure and the rebuilding of a favourable pressure gradient field. In due time, the imploding dense core band will again reverse its motion and commence a second expansion. The synchronization between the central pressure fluctuation and the ring pulsation is termed macroscale gas–particle coupling which comes to an end when either the ring disintegrates (figure 5d,e) or the kinetic energy of the ring is damped away (figure 5f).

The sustainability of the macroscale gas–particle coupling in terms of the duration of the synchronization between the central gases and the enclosing ring, depends on three pivotal processes that occur on the macro-, meso- and microscale, as elaborated below. The first involves the initial momentum/energy transfer from the explosive source (central gases) to the enclosing ring. If the ring expansion is excessively fast, which occurs with small M/C, the central negative pressure field barely influences the dynamics of the ring, as shown in figure 5(a). The macroscale gas–particle coupling is short lived, only serving to provide the initial impetus to the particle ring. In this scenario, the time scale of the ring dynamics, which is defined as the time the ring takes to expand to twice the initial diameter, t_ring, is one order smaller than that of the pressure field, which is characterized by the time the central pressure takes to decrease to the ambient pressure, t_gas, namely τ_macro = t_ring/t_gas ~ O(10⁻¹). From figure 5(a–d), τ_macro = 0.6, 1.2, 4.9, 12.5, indicating an improved compatibility between the ring dynamics and the gas evolution with increasing M/C. Note that τ_macro does not exist for the systems 1034-100-2500 and 2043-100-2500 since the ring cannot expand to twice its initial diameter prior to the implosion.

3.3. Mesoscale gas–particle shell coupling

The initial momentum/energy transfer from the central gases to the particle ring occurs during the shock compaction phase which begins when the incident shock impinges on the inner surface of the ring and ends when the outer surface gains velocity. Two fundamental processes dominate the shock compaction, as schematized in figure 6(a,b). One is the transmission of blast loading exerted on the inner surface of the ring via the solid stresses arising from inter-particle contacts. The other is the inter-pore gas filtration driven by the pressure differential between the inner and outer surfaces of the ring. Alongside the inter-pore gas flows, a diffusional pressure field develops. Particles hence experience the pressure gradient forces and drag forces besides the solid stresses. All these forces together accelerate the particles and compact them to a compacted packing. Since the momentum/energy transportation is caused by the collective particle movements and inter-pore gas flows, the shock compaction is governed by the mesoscale gas–particle coupling.

Figure 6. Schematics of compaction wave propagation via inter-particle stresses (a) and gas filtration through pores in the granular packing (b) induced by the impingement of the incident shock.

While the particle ring gains the momentum, the interstitial gas flows cause substantial mass and energy losses of the central gases, thereby reducing the energy transferred to the ring. The resulting mass/energy reduction varies from system to system due to the distinctively different gas filtration. Figure 7 shows the early-time space–time (r–t) diagrams of ϕ_p (figure 7a,b), P_gas (figure 7c,d) and u_gas (figure 7e,f) for systems 1035-50-1828 and 1035-800-1828. The trajectories of the compaction front (CF) and the diffusional pressure front (DPF) are superimposed on each panel in figure 7. The CF delineates the compacted particle packing. As the CF reaches the outer surface of the ring, the shock compaction phase is completed with the compacted ring acquiring an initial expansion velocity V_ring. The distance of the DPF from the inner surface of the ring, L_DPF = R_DPF − R_in, characterizes the propagation depth of the transient diffusional pressure field whose determination is presented in Appendix C. For the system 1035-50-1828, the DPF propagates far slower than the CF, thereby there is no interstitial gas flowing out of the outer surface of the ring during the shock compaction phase (figure 7e). By contrast, the DPF distinctly precedes the CF in the system 1035-800-1828, inducing the intensified and persistent gas flows (figure 7f). Accordingly, the mass and energy losses mediated by the gas flowing out vary significantly. The mass ratios of gases retained in the central gas pocket after the shock compaction phase, χ_gas, are 96 % and 73 % for the systems 1035-50-1828 and 1035-800-1828, respectively. As expected, the ring in the system 1035-800-1828 gains a faster initial expanding velocity than that in the system 1035-50-1828, V_{ring, 1035-50-1828} = 8 m s⁻¹, V_{ring, 1035-800-1828} = 5 m s⁻¹.

Figure 7. Early-time space–time (r–t) diagrams of ϕ_p (a,b), P_gas (c,d) and u_gas (e,f) for systems 1035-50-1828 (a,c,e) and 1035-800-1828 (b,d,f). The trajectories of the CF and the DPF are superimposed on each panel.

Since the mass and energy losses of the central gases are associated with the propagation velocity ratio between the DPF and the CF, τ_meso = V_DPF/V_CF, if τ_meso can be expressed as a function of structural parameters, the energy diminishing effect arising from the mesoscale gas–particle coupling can be estimated for a given dispersal system. The advance of the DPF depends on the velocity of the local interstitial gas flows relative to the solid skeleton, V_DPF = u_gas(R_DPF) − u_solid(R_DPF). Ignoring the nonlinear terms, the interstitial gas flows can be accounted for by the Darcy equation, whose expression in polar coordinates is given by (3.1)

(3.1)\begin{equation}\varepsilon ({u_{gas}} - {u_{solid}}) =- \frac{k}{\mu }\frac{{\textrm{d}{P_{gas}}}}{{\textrm{d}r}},\end{equation}

where ε is the local porosity, μ is the dynamic viscosity of the interstitial gases and the permeability k is a function of the volume fraction ϕ_p described by the Ergun expression

(3.2)\begin{equation}k = \frac{1}{{150}}\frac{{{\varepsilon ^3}}}{{{{(1 - \varepsilon )}^2}}}d_p^2.\end{equation}

Thus, V_DPF is proportionate to the local gradient of the transient diffusional pressure field at the DPF, $(\textrm{d}{P_{gas}}/\textrm{d}r){|_{{R_{DPF}}}}$, which varies as the DPF advances. As a first order of approximation, we assume a linearly declining pressure field across the thickness of the ring, leading to an unvaried pressure gradient

(3.3)\begin{equation}\frac{{\textrm{d}{P_{gas}}}}{{\textrm{d}r}} =- \frac{{{P_5} - {P_{amb}}}}{{{h_{ring}}}},\end{equation}

where the initial pressure exerted on the inner surfaces of the ring is represented by the reflected pressure P ₅ which is estimated by the normal shock wave relation with the initial pressure P ₀. The ambient pressure P_amb represents the pressure exerted on the inner and outer surfaces of the ring. As a matter of fact, the central pressure decreases due to the expansion of the inner surface and the filtration induced mass flow out. Thus, the constant pressure assumption overestimates $\textrm{d}{P_{gas}}/\textrm{d}r$, which is somewhat offset by replacing the actual propagation depth of the diffusional pressure field with the thickness of the ring. Substituting (3.2)–(3.3) into (3.1) and replacing ϕ_p with ϕ ₀ to have

(3.4)\begin{equation}{V_{DPF}} = \frac{{d_p^2}}{{150}}\frac{{{{(1 - {\phi _0})}^2}}}{{\mu \phi _0^2}}\frac{{{P_5} - {P_{amb}}}}{{{h_{ring}}}}.\end{equation}

Equation (3.4) specifies the structural parameters affecting V_DFP, including the macroscale parameters such as P ₀, ϕ ₀ and h_ring, as well as the microscale parameter d_p. The time derivative of R_DPF yields an absolute propagation velocity of the DPF in contrast with the relative velocity implied by the V_DPF given in (3.4). Thus, we subtract the average particle velocity upon the DPF from the average R_DPF to derive V_DPF for numerical cases (details can be found in Appendix C). For the systems 1035-50-1828 and 1035-800-1828, the numerically derived V_DPF are 2.3 and 606 m s⁻¹, respectively, which agree well with the respective theoretical predictions given in (3.4), 2.6 and 658.2 m s⁻¹.

As to the V_CF, our previous work (Kun et al. Reference Kun, Lvlan, Jiarui, Chunhua and Baolin2023) proposed a shock compaction model accounting for the propagation of the CF in the particle packings driven by the diffusional pressure field. This model yields an adequate prediction of the V_CF given in (3.5) for the densely packed granular column with the upstream and downstream pressure kept constant

(3.5)\begin{equation}{V_{CF}} = \sqrt {\frac{{{P_5} - {P_{amb}}}}{{{\rho _p}}}\frac{{{\phi _{comp}}}}{{({\phi _{comp}} - {\phi _0}){\phi _0}}}} ,\end{equation}

where ϕ_comp is the maximum volume fraction reached by an assembly of particles with a certain size distribution, here ϕ_comp = 0.7. Although the studied configuration has a divergent geometry and the pressure on the inner surface of the ring is unsteady, (3.5) still provides a reasonable estimation for the V_CF. With (3.5), the systems 1035-50-1828 and 1035-800-1828 share the same V_CF of 63.4 m s⁻¹, compared with the respective numerical results, 64.8 and 54.2 m s⁻¹.

Combining (3.4) and (3.5) we obtain the expression of τ_meso as a function of a variety of structural parameters

(3.6)\begin{equation}\begin{array}{l} {\tau _{meso}} = \dfrac{{{V_{DPF}}}}{{{V_{CF}}}} = \dfrac{1}{{150}}\dfrac{{{{(1 - {\phi _0})}^2}}}{{\phi _0^2}} \cdot \sqrt {\dfrac{{({\phi _{comp}} - {\phi _0}){\phi _0}}}{{{\phi _{comp}}}}} \cdot \dfrac{{d_p^2\sqrt {({P_5} - {P_{amb}}){\rho _p}} }}{{\mu {h_{ring}}}}.\\ \end{array}\end{equation}

Figure 8 plots χ_gas for all systems with varying τ_meso. The symbols are rendered according to the corresponding M/C. An explicit τ_meso dependence of χ_gas is discernible as τ_meso ranges from O(10⁻²) to O(10¹). At the lower limit of the τ_meso range, τ_meso ~ O(10⁻²), the filtration is far slower than the compaction. The value of χ_gas is close to unity, indicating a minimum gaseous mass flow out. A slightly downward slope of χ_gas emerges when the τ_meso approaches O(10⁰), the gas filtration beginning to drain the central gases. In the mid- and upper range, τ_meso ~ O(10⁰) − O(10¹), χ_gas rapidly decreases with increasing τ_meso. The correlation between χ_gas and τ_meso can be described by a fitting function

(3.7) \begin{align} {\chi _{gas}} & = f(\textrm{lg}\,{\tau _{meso}}) = 0.0227{(\textrm{lg}\,{\tau _{meso}})^3} - 0.172{(\textrm{lg}\,{\tau _{meso}})^2}\notag\\ & \quad - 0.111\,\textrm{lg}\,{\tau _{meso}} + 0.933,{\tau _{meso}}\sim O({10^{ - 2}}) - O({10^1}). \end{align}

Figure 8. Variation in χ_gas with τ_meso. symbols are rendered according to the corresponding M/C. A master curve fitted by a polynomial function with an 80 % confidence interval is superimposed in the plot.

Notably, χ_gas for systems at the lower limit of the investigated M/C range, M/C ~ O(10¹), exhibits non-trivial deviations from the fitting master curve. As revealed in the last section, the expansion of the ring in this M/C range is exceedingly fast so that the central pressure undergoes a drastic drop. The constant surface pressure assumption underpinning (3.4) and (3.5) does not hold. Thus, the correlation between χ_gas and τ_meso can no longer be adequately described in (3.7). Since χ_gas is negligible (χ_gas > 92 %) in these cases, mainly thanks to the quite short duration of the shock compaction phase, we assume χ_gas = 1, ignoring the gas filtration altogether.

The gas filtration associated with the mesoscale gas–particle coupling reduces the effective mass in the central gas pocket, equivalent to an increased mass ratio, $\textrm{M}/(\textrm{C} \cdot {\chi _{gas}}) = \textrm{M}/(\textrm{C} \cdot f(\textrm{lg}\,{\tau _{meso}}))$ whereby the influence of the mesoscale coupling on the dispersal is properly accounted for. Although the derivation of (3.7) is based on the simulation data, the τ_meso dependence of the χ_gas can also be constructed via a theoretical model introduced in § 4.

3.4. Microscale gas–particle shell coupling

In contrast with the solid or liquid shells/rings, the particle shedding and the variations of volume fraction are inevitable throughout the dispersal of the particle shell/ring. Since the shed-off particles and particles in regions with very low ϕ_p (ϕ_p < 0.3) neither sustain persistent inter-particle contacts nor effectively confine the central gases, they cannot be regarded as part of the densely packed particle shell/ring. The actual particle ring whose pulsation synchronizes with the central pressure fluctuation hence has smaller mass than the initial one, equivalent to having a reduced M/C. Therefore, the particle shedding and localized pack loosening may well fundamentally change the dispersal behaviour. Since these events involve the individual particle dynamics, the resulting thinning of the effective particle ring embodies the microscale gas–particle coupling. In this section, we first explore the mechanism governing the particle shedding and pack loosening, then proceed to quantify the influence of these microscale processes on the overall dispersal and finally shed light on the outliers in figure 4(a).

The integrity of particle ring is maintained during the shock compaction phase and the incipient expansion phase when a favourable pressure gradient is imposed on the particles. The direction of the pressure gradient is soon reversed as the central pressure drops to a negative value. Figure 9(a,b) shows the radial profiles of the pressure and pressure gradient, P_gas(r) and ${\nabla _r}{P_{gas}}(r)$, during the overexpansion phase of the central gases for the system 288-60-1503 whose r–t diagram of ϕ_p is shown in figure 2(d). The corresponding radial profiles of the particle velocity and volume fraction, u_p(r) and ϕ_p(r), are presented in figure 9(c,d), respectively. The substantial central pressure decline emanates expansion waves into the particle ring, unloading the compacted particles in the wake. The particles adjacent to the inner surface of the ring become so loosely packed (figure 9d) that no appreciable pressure gradient is built therein (figure 9b), inducing a dilute inner band. Meanwhile an outermost layer of the ring pells off when the CF reflects off the outer surface and transitions to an inward travelling rarefaction wave. Particles swept by the rarefaction wave are prone to breaking away from the outer surface, forming a shedding outer band. In between the dilute inner band and the shedding outer band there exists a densely packed core band whose ϕ_p(r) peaks at the band centre and this is denoted as $\textrm{B}{\textrm{C}_{{\phi _{max}}}}$ in figure 9(d) with the red dotted line. Meanwhile, an adverse pressure gradient field with the radial profile resembling that of ϕ_p(r) is established across the thickness of the dense core band. The amplitude of the inward pointing pressure gradient forces, $|{{F_{{\nabla_P}}}} |$, increases on approaching $\textrm{B}{\textrm{C}_{{\phi _{max}}}}$ and decreases on moving away from $\textrm{B}{\textrm{C}_{{\phi _{max}}}}$ as expressed in (3.8)

(3.8)\begin{equation}\left. {\begin{array}{*{20}{c@{}}} {|{{F_{\nabla P}}} |\sim 0,\quad r < {R_{dense,in}}\ \textrm{or}\ r > {R_{dense,out}},}\\ {\dfrac{{\partial |{{F_{\boldsymbol{\nabla }P}}} |}}{{\partial r}} > 0,\quad {R_{dense,in}} \le r < {R_{{\phi_{max}}}},}\\ {\dfrac{{\partial |{{F_{\boldsymbol{\nabla }P}}} |}}{{\partial r}} < 0,\quad {R_{{\phi_{max}}}} \le r \le {R_{dense,out}}.} \end{array}} \right\}\end{equation}

Where R_dense,in, ${R_{{\phi _{max}}}}$ and R_dense,out represent the radii of the inner boundary, the centre ($\textrm{B}{\textrm{C}_{{\phi _{max}}}}$) and the outer boundary of the densely packed core band. The drag forces, F_drag, which are proportionate to the pressure gradient forces in the dense packings, have identical radial profiles (Kun et al. Reference Kun, Lvlan, Jiarui, Chunhua and Baolin2023).

Figure 9. Variations in radial profiles of P_gas(r) (a), ${\nabla _r}{P_{gas}}(r)$ (b), u_p(r) (c) and ϕ_p(r) (d) for the system 288-60-1503 during the overexpansion phase of central gases, 5.3 ms < t < 8.9 ms. The dotted lines in each plot from left to right represent the inner boundary of the particle cloud, the inner boundary of the dense core band, the band centre $\textrm{B}{\textrm{C}_{{\phi _{max}}}}$, the outer boundary of the dense core band and the outer boundary of the particle cloud, respectively. Typical radial profiles of ${\nabla _r}{P_{gas}}(r)$, u_p(r) and ϕ_p(r) are shown in (e).

Such distinctive profiles of ${F_{{\nabla _P}}}$ and F_drag result in a non-monotonic radial variation in u_p(r) (figure 9c,e). During the ring implosion phase, in the inner half of the dense core band (${R_{dense,in}} \le r < {R_{{\phi _{max}}}}$), the absolute value of the u_p(r) increases with r. Particles in the outer layers tightly compress against those in the inner layers. Thereby, the inner half of the dense core band remains densely packed. By contrast, the absolute value of u_p(r) decreases with r in the outer half of the dense core band (${R_{{\phi _{max}}}} \le r \le {R_{dense,out}}$), leading to a persistent pack loosening and eventually particle shedding. Notably, the imploding dense core band arrests the drifting particles in the dilute inner band, countering the particle shedding from the outer surface.

Figure 10 shows the temporal mass percentage of particles collected by the inner boundary of the dense core band, χ_add(t), and those breaking loose from the outer boundary, χ_shed(t), during the ring implosion. The variation in the thickness of the dense core band, h_dense(t), is also superimposed in figure 10. For the system 288-60-1503, the particle shedding effect always prevails over the particle collecting effect. The ring completely disintegrates when the former depletes the entire dense core band at ${R_{{\phi _{max}}}} = 0.134$. In this scenario, a disastrous dispersal, where a majority of particles collapse into the centre, is observed.

Figure 10. Temporal variations in χ_add(t), χ_shed(t) and h_dense(t) during the ring implosion phase for the system 288-60-1503.

The ring thinning effect originating from the particle shedding depends on the dynamics of the ring as well as the individual particle. Specifically, particles of small inertia are more prone to reverse their motion and less likely to break loose compared with those of large inertia. On the other hand, a short-lived implosion of ring does not allow particles enough time to complete the motion reversal so that the ring sheds more severely than the long-lasting one. Thus, it is necessary to quantify the time scales of the ring implosion as well as the particle motion reversal.

Assuming an outbound particle with density of ρ_p and diameter of d_p is about to reverse its motion driven by an inward directing pressure gradient force, the pressure gradient can be estimated by assuming a linear decline of pressure across the thickness of the ring, namely $\boldsymbol{\nabla }P = {P_{amb}}/{h_{ring}}$. The characteristic time such a particle takes to reverse its motion is

(3.9)\begin{equation}{t_p} = \sqrt {\frac{{2{d_p}{h_{ring}}{\rho _p}}}{{{P_{amb}}}}} .\end{equation}

The derivation of t_p is presented in Appendix D. Smaller t_p is, the faster the outbound particle reverses to inbound and the less likely it is to break loose.

As seen from figure 2(e,f), the expansion–implosion cycle of the ring is semi-symmetric in terms of the duration and the amplitude as well. Hence, the time scale of the ring implosion can be approximated by that of the ring expansion introduced in § 3.2, namely the t_ring which is expressed in (3.10) under the assumption of a constant expansion velocity

(3.10)\begin{equation}{t_{ring}} = \frac{{{R_{out,0}}}}{{{V_{ring,Gurney}}}}.\end{equation}

The denominator in (3.10), V_{ring ,Gurney}, is the estimation of V_ring by a modified Gurney equation (Mo et al. Reference Mo, Lien, Zhang and Cronin2019) which accounts for the effect of the volume fraction ϕ ₀ in addition to the structural parameters normally considered in a conventional Gurney equation

(3.11)\begin{gather}{V_{ring,Gurney}} = \sqrt {\frac{{2{e_{gas}}}}{{\dfrac{{\textrm{M/C}}}{{0.2\rho _p^{0.18}}} + 0.5}}} \cdot F({\phi _0},\textrm{M/C}),\end{gather}

(3.12)\begin{gather}F({\phi _0},\textrm{M/C}) = 1 + [0.162\,\textrm{exp}(1.127{\phi _0}) - 0.5] \cdot \textrm{lg}(\textrm{M/C}).\end{gather}

The e_gas in (3.11) is the explosion energy per unit mass of the central gases

(3.13)\begin{equation}{e_{gas}} = \frac{{{E_{gas}}}}{{{m_{gas}}}} = \frac{{({P_0} - {P_{amb}})R{T_0}}}{{(\gamma - 1){P_0}}} \simeq \frac{{R{T_0}}}{{(\gamma - 1)}}.\end{equation}

The ratio between t_ring and t_p, τ_micro = t_ring/t_p, measures the importance of the microscale particle shedding relative to the macroscale ring implosion. Combining (3.9)–(3.13) leads to the expression of τ_micro which is a function of both macroscale and microscale parameters

(3.14)\begin{gather}{\tau _{micro}} = {I_p}({d_p},{\rho _p}) \cdot {F_{macro}}(\textrm{M/C},{\rho _{gas}},{R_{out,0}},{h_{ring}},{\phi _0}),\end{gather}

(3.15)\begin{gather}{I_p} = \frac{1}{{d_p^{0.5}\rho _p^{0.5}}},\end{gather}

(3.16)\begin{gather}{F_{macro}} = \frac{{{R_{out,0}}}}{{2F({\phi _0},\textrm{M/C})}}\sqrt {\left( {\frac{{\textrm{M/C}}}{{0.2\rho_p^{0.18}}} + 0.5} \right)\frac{{{P_{amb}}}}{{{e_{gas}}{h_{ring}}}}} .\end{gather}

A smaller τ_micro, or equivalently a larger 1/τ_micro, indicates a longer particle motion reversal time relative to the entire implosion duration. Particles are more likely to remain outbound, breaking loose from the imploding ring. This argument is corroborated by the τ_micro dependence of the cumulative volume fraction of shed particles at the end of the first ring implosion, χ_shed, shown in figure 11. Indeed, χ_shed significantly increases as τ_micro decreases from O(10²) to O(10¹). Once τ_micro falls below 20, χ_shed converges to unity, indicating a complete falling apart of the ring during the first implosion. Notably, the symbols in figure 11 are rendered according to the respective M/C. A M/C dependence of τ_micro is discernible in the mid- and upper range of the M/C, M/C ~ O(10³) − O(10⁴). Specifically, a larger M/C results in a higher τ_micro, suggesting the ring thinning effect caused by the particle shedding becomes increasingly insignificant with increasing M/C.

Figure 11. The τ_micro dependence of χ_shed. Symbols are rendered according to M/C. The master curve with an 80 % confidence interval is superimposed in the plot.

As seen in (3.14), τ_micro is the product of a particle inertia related factor I_p (3.15) and a macroscale complex F_macro (3.16) which is a function of the M/C and the ring geometry. Decreasing d_p and/or ρ_p leads to an elevated I_p whereby τ_micro increases. The particle mass loss caused by the particle shedding will be suppressed. Equivalently, the effective ring in motion is thicker than its counterpart with a smaller τ_micro. The resulting dispersal behaviour resembles that with an amplified M/C. Figure 12 exhibits the r–t diagrams of ϕ_p for systems 55-100-3506 (figure 12a), 55-450-417 (figure 12b), 55-250-210 (figure 12c) which show worsening completeness and uniformity of dispersal, as first alluded to in § 3.1. These three systems have the same M/C, similar ring geometry but different d_p and ρ_p, and thereby different I_p. As I_p decreases from 1.7 (d_p = 100 μm, ρ_p = 3506 kg m⁻³) to 2.3 (d_p = 450 μm, ρ_p = 417 kg m⁻³), then to 4.4 (d_p = 250 μm, ρ_p = 210 kg m⁻³), the disintegration of the ring is significantly retarded and leads to a lower proportion of particles, χ = 0.99, 0.8, 0.3 for these three systems, respectively. This complete-to-incomplete dispersal transition caused by different microscale couplings is responsible for the bulk of the outliers in figure 4(a).

Figure 12. The r–t diagrams of ϕp for systems 55-100-3506 (a), 55-450-417 (b), 55-250-210 (c). The inner and outer boundaries of the particle clouds are denoted by white dashed lines.

3.5. Structure-dispersal correlation

The previous sections shed light on the multiscale gas–particle coupling relations underlying the varying explosive dispersal behaviours. The influences of the macro-, meso- and microscale coupling on the dispersal behaviour can be characterized via the parameter complex, M/C (2.13), τ_meso (3.6) and τ_micro (3.14), respectively. These non-dimensional complexes incorporate a range of macro- and microscale structural variables. The M/C plays a primary role in determining the overall dispersal performance in terms of χ and κ. The parameters τ_meso and τ_micro bring about the complementary but indispensable influences related to the mass reduction of gases and particles, respectively. In this section, we aim to construct a defining non-dimensional complex which adequately manifests the correlative effects arising from the M/C, τ_meso and τ_micro, thereby bridging the dispersal systems to the resulting dispersal performances.

An explicit τ_meso dependence of χ_gas, χ_gas(τ_meso) is established in § 3.3, as presented by (3.7), which is validated throughout the investigated M/C range except at the lower limit of the range, M/C ~ O(10⁰ − 10¹). Beyond that range, the original M/C ought to be replaced by M/(C⋅χ_gas(τ_meso)) to account for the influence of τ_meso. As to τ_micro, although it is clear that the diminishing mass of the particle ring in motion becomes stronger with a smaller τ_micro, equivalent to having a reduced M/C, the exact correlation between them is beyond the scope of the present study. Thus, to incorporate the microscale coupling effect we tentatively modify the M/C by multiplying it by lg(τ_micro). Finally, we construct a modified mass ratio, (M/C)*, which is formulized as

(3.17)\begin{equation}{(\textrm{M/C})^\ast } = \left\{ {\begin{array}{*{20}{@{}ll@{}}} {(\textrm{M/C}) \cdot \textrm{lg}({\tau_{micro}})}&{\textrm{(M/C)}\sim O({{10}^0} - {{10}^1})}\\ {\textrm{(M/C}) \cdot \dfrac{1}{{{\chi_{gas}}({\tau_{meso}})}} \cdot \textrm{lg}({\tau_{micro}})}&{\textrm{(M/C}) \ge O({{10}^1})} \end{array}} \right..\end{equation}

The correspondence between the M/C and (M/C)* for all systems is presented in figure 13. The size of the symbols scales with τ_meso, while the colour of the symbols is rendered according to τ_micro. Where the M/C is of the order of O(10³), the (M/C)* is proportionate to the M/C as the effects of τ_micro and τ_meso are minimum. Otherwise in the lower- and mid-range of M/C, M/C ~ O(10¹) − O(10²), the (M/C)* of systems with the upper limit of τ_meso and/or the lower limit of τ_micro significantly deviate from the scaling law.

Figure 13. Correlation between M/C and (M/C)*. The size of symbols scales with τ_meso, while the colour of symbols is rendered according to τ_micro.

Replacing the M/C with the (M/C)*, we replot figure 4(a,b) in figure 14(a,b). In contrast with the wide data variabilities shown in figure 4(a,b), the data converging into the master curves of χ((M/C)*) and κ((M/C)*) is substantially improved without apparent outliers, which corroborates the decisive role of (M/C)* in determining the dispersal characteristics. As seen in figure 14(a,b), the dispersal remains complete (χ ~ 1) and uniform (κ remains minimum), namely an ideal dispersal, until the (M/C)* is of the order of O(10²). Thereafter, dispersal incompleteness quickly exacerbates as χ begins a sudden and substantial decline, where χ drops to below 0.5 over one decade of (M/C)*. Meanwhile, κ increases moderately. Accordingly, the explosive dispersal transitions from a complete and uniform one into an incomplete one with a non-trivial proportion of particles failing to disperse out. Regardless, particles in the dispersed particle cloud remain mostly uniformly distributed as suggested by the low value of κ (κ < 0.1). Thus, this dispersal mode is referred to as a partial dispersal. When the (M/C)* reaches the order of O(10³) and beyond, the decreasing rate of χ significantly slows down, gradually approaching zero as (M/C)* approaches O(10⁴). By contrast the growth rate of κ expedites, indicating a quickly worsening uniformity of dispersal. In this scenario the majority of particles dwell in a highly concentrated region, giving rise to a failed dispersal.

Figure 14. Variations in χ (a) and κ (b) with increasing (M/C)*. The size of symbols scales with the 1/d_p, while the colour of symbols is rendered according to ρ_p.

4.1. Criteria for the dispersal mode transition

As revealed in § 3.5, the dispersal mode transitions from ideal to partial then to failed as the (M/C)* increases from O(10⁻¹) to O(10⁴). The (M/C)* thresholds signifying the ideal-to-partial and partial-to-failed mode transitions are denoted as $(\textrm{M/C})_{th,I}^\ast $ and $(\textrm{M/C})_{th,II}^\ast $, as marked in figure 14(a,b), $(\textrm{M/C})_{th,I}^\ast{\sim} 1 \times {10^2}$, $(\textrm{M/C})_{th,II}^\ast{\sim} 1 \times {10^3}$, respectively. Since the exact values of the $(\textrm{M/C})_{th,X}^\ast $ (X = I or II) are sensitive to the numerical models and parameter set-ups, the $(\textrm{M/C})_{th,X}^\ast $ (X = I or II) marked in figure 14(a,b) only serves as the reference rather than a guaranteed guide to predicting the dispersal mode. Instead of resorting to a single defining parameter, in this section we attempt to establish a theoretical model accounting for the entire ring dynamics resulting from a specific dispersal system. Combed with reasonable criteria for the dispersal mode transitions, the model allows for the prediction of the dispersal mode.

Firstly, the criterion for the ideal-to-partial dispersal mode transition is introduced. As discussed in § 3.2, in the ideal dispersal the time scale of the ring expansion is far smaller than that of the gas evolution, τ_macro ~ O(10⁻¹). Figure 15(a) shows the variation in τ_macro with the (M/C)*. Indeed all systems with (M/C)* falling below $(\textrm{M/C})_{th,I}^\ast $ satisfy this requirement. Thereafter τ_macro increases to O(10⁰), the time scale of the ring dynamics becoming comparable to that of the pressure evolution, leading to a vast number of particles being sucked into the centre by the negative central pressure. Therefore τ_macro ~ O(10⁻¹) serves as a requirement for an ideal dispersal.

Figure 15. Variations in τ_macro (a) and ξ_imp (b) with the (M/C)*. Dark red symbols in (b) represent the semi-failed dispersal while light red symbols represent validated failed dispersal.

At the other end of the dispersal mode spectrum, a failed dispersal featuring sustained macroscale coupling requires a complete ring to implode at least once. The likelihood of the particle ring surviving a complete implosion depends on the thickness of the ring as well as the implosion distance. A thinner ring with a larger radius is less likely to survive a complete implosion (see figure 2c–f). Figure 15(b) plots the ratio between the thickness and the outer radius of the ring at the very moment of motion reversal (see the inset of figure 15b), ξ_imp = h_dense,imp/R_dense,imp, for systems of the failed dispersal mode. The ring here refers to the dense core band introduced in § 3.4. The systems wherein the rings do not survive the first implosion are denoted by the dark red circles. Otherwise, the systems are denoted by the light red circles. Thereafter, these two groups of systems are referred to as the semi- and validated failed dispersal. The former clusters just below the $(\textrm{M/C})_{th,II}^\ast $ while the latter is above the $(\textrm{M/C})_{th,II}^\ast $. Clearly, the ξ_imp of all the light red circles exceeds 0.1, thereby systems with ξ_imp > 0.1 suffice to predict a validated failed dispersal.

4.2. Continuum-based explosive dispersal model

The last section sets the criteria for the ideal and validated failed modes, namely τ_macro ~ O(10⁻¹) and ξ_imp > 0.1, respectively. To derive τ_macro and ξ_imp, it is necessary to adequately reproduce two primary ring dynamics phases, namely the shock compaction and the ring pulsation phases, which are the focus of two sub-models. The shock compaction sub-model aims to account for CF which propagates through the thickness of the ring and couples with the transient interstitial gas flows, thereby estimating the kinetic energy imparted to the particle rings and the mass/energy loss of the central gases at the end of shock compaction. The final state of the shock compaction sub-model serves as the initial state of the following ring pulsation sub-model. The ring pulsation sub-model intends to yield the key dynamic parameters pertinent to the ring and central gases, such as t_ring, t_gas, h_dense,imp and R_dense,imp, whereby τ_macro and ξ_imp are determined. A value of τ_macro of the order of O(10⁻¹) results in the ideal dispersal mode. Otherwise, ξ_imp is assessed sequentially. If ξ_imp > 0.1 holds, a validated failed dispersal is expected. Otherwise the resulting dispersal is either partial or semi-failed.

4.2.1. Evolution of central gases and gas filtration

Considering the configuration of the archetype shown in figure 1(a), wherein a circular central gas pocket with the initial radius R_{in ,0} is filled with homogenous static gases with an initial pressure P ₀ and density ρ_{centre ,0}

(4.1)\begin{equation}{\rho _{center,0}} = \frac{{{P_0}}}{{\bar{R}{T_0}}},\end{equation}

where $\bar{R}$ is the specific gas constant, $\bar{R} = 288.7\ \textrm{J}\ \textrm{k}{\textrm{g}^{ - 1}}\ {\textrm{k}^{ - 1}}$. The enclosing ring with the initial inner and outer radii of R_{in ,0} and R_{out ,0} comprises the compressible porous medium with the initial porosity ε ₀ and the solid phase density ρ_solid.

Throughout the shock compaction and the following ring pulsation phases, gases in the central gas pocket are allowed to flow out (the central pressure is positive) or flow in (the central pressure is negative) due to the inter-pore gas filtration which is governed by Darcy's law (3.1). The mass loss or gain rate of the central gases depends on the gas velocity on the inner surface of the ring relative to the moving inner surface

(4.2)\begin{equation}{\dot{m}_{gas}}(t) = 2{\rm \pi} {R_{in}} \cdot {\rho _{gas}}({R_{in}}) \cdot [{u_{gas}}({R_{in}}) - {V_{in}}] \cdot \varepsilon ({R_{in}}),\end{equation}

where ρ_gas(R_in), u_gas(R_in) and ε(R_in) are the instantaneous interstitial gas density, gas velocity and the porosity on the inner surface of the ring, respectively, and V_in is the expanding velocity of the inner surface. Since the positive u_gas directs outward, a positive the mass loss rate induces a negative mass gain rate. The instantaneous mass retained inside the central gas pocket is

(4.3)\begin{equation}{m_{gas}}(t) = {m_{gas,0}} - \int_0^t {2{\rm \pi} {R_{in}} \cdot {\rho _{gas}}({R_{in}}) \cdot [{u_{gas}}({R_{in}}) - {V_{in}}] \cdot \varepsilon ({R_{in}})\,\textrm{d}t} ,\end{equation}

where m_{gas ,0} = ${\rm \pi}$R²_{in ,0}⋅ρ_{centre ,0} is the initial mass of the central gases. The cumulative mass flux combined with the volumetric variation of the central gas pocket caused by the expansion or implosion of the ring collectively leads to the density evolution of the central gases

(4.4)\begin{equation}\begin{array}{ccccc} {\rho _{center}} & = \dfrac{{R_{in,0}^2}}{{R_{in}^2}}{\rho _{center,0}} - \int_0^t {\dfrac{{2{\rho _{gas}}({R_{in}})}}{{{R_{in}}}} \cdot [{u_{gas}}({R_{in}}) - {V_{in}}] \cdot \varepsilon ({R_{in}})\,\textrm{d}t} \\ & = \dfrac{{R_{in,0}^2}}{{R_{in}^2}}{\rho _{center,0}} + \int_0^t {\dfrac{{2{\rho _{gas}}({R_{in}})}}{{{R_{in}}}} \cdot \dfrac{k}{\mu }\boldsymbol{\nabla }{P_{{R_{in}}}}\,\textrm{d}t} . \end{array}\end{equation}

Assuming an isentropic expansion for the central gases gives

(4.5)\begin{equation}\frac{{{P_{center}}}}{{\rho _{center}^\gamma }} = \frac{{{P_0}}}{{\rho _{center,0}^\gamma }},\end{equation}

where γ is the specific heat ratio, γ = 1.4. Substituting (4.1) and (4.4) into (4.5) to derive the instantaneous central pressure, which is also the pressure exerted on the inner surface of the ring, P_gas(R_in) = P_centre

(4.6)\begin{equation}{P_{gas}}({R_{in}}) = {P_{center}} = \frac{{{{\bar{R}}^\gamma }T_0^\gamma }}{{P_0^{\gamma - 1}}} \cdot {\left[ {\frac{{R_{in,0}^2}}{{R_{in}^2}}\frac{{{P_0}}}{{\bar{R}{T_0}}} + \int_0^t {\frac{{2{\rho_{gas}}({R_{in}})}}{{{R_{in}}}} \cdot \frac{k}{\mu }\boldsymbol{\nabla }{P_{{R_{in}}}}\,\textrm{d}t} } \right]^\gamma }.\end{equation}

As seen in (4.6), P_center depends on the ring dynamics and the gas filtration as well. The former dictates the volumetric variation of the central gas pocket via the variable R_in. The latter involves the dynamic flow parameters, u_gas(R_in) and ρ_gas(R_in), which are incorporated into the Morrison equation in polar coordinates

(4.7)\begin{equation}\frac{{\partial (\varepsilon {P_{gas}})}}{{\partial t}} = \frac{{\partial \left( {\dfrac{k}{\mu }\dfrac{{\partial {P_{gas}}}}{{\partial r}}{P_{gas}}r - \varepsilon {u_{solid}}{P_{gas}}r} \right)}}{{r\partial r}},\end{equation}

where u_solid is the velocity of the solid phase, V_in = u_solid(R_in). The boundary conditions of (4.7) are set by the pressure on the inner and outer surfaces of the ring, P_gas(R_in) = P_center and P_gas(R_out) = P_amb. Since the evolution of the central gases (4.6), the gas filtration (4.7) and the ring dynamics are coupled together, their governing equations ought to be solved together.

4.2.2. Shock compaction sub-model

Upon the release of the pressurized central gases, a diffusional pressure field develops across the thickness of the ring and progresses outwards, as demonstrated in figure 16(a). Meanwhile, the interstitial gases start to flow, driven by the diffusional pressure field. The solid phase of the ring immersed in the advancing diffusional pressure field is subjected to the outward pressure gradient force ${F_{{\nabla _P}}}$ and the associated drag force F_drag, both of which scale with the local pressure gradient

(4.8)\begin{gather}{F_{\boldsymbol{\nabla }P}} =- \frac{{\textrm{d}{P_{gas}}/\textrm{d}r}}{{{\rho _{solid}}}},\end{gather}

(4.9)\begin{gather}{F_{drag}} = \frac{\varepsilon }{{1 - \varepsilon }}{F_{\boldsymbol{\nabla }P}} =- \frac{\varepsilon }{{1 - \varepsilon }}\frac{{\textrm{d}{P_{gas}}/\textrm{d}r}}{{{\rho _{solid}}}}.\end{gather}

Figure 16. Schematics of configuration considered in the theoretical model. (a) A CF has a finite thickness of w_CF across which the porosity ε and solid velocity u_solid gradually decrease. (b) Propagation of the CF with a velocity of V_CF driven by ${F_{{\nabla _P}}}$ and F_drag. The radial profiles of ε(r), u_solid(r), ${F_{{\nabla _P}}}(r)$ and F_drag(r) are superimposed in the plots.

These forces compact the solid skeleton to a compacted band with a minimum porosity ε_min, ε_min = 1 − ϕ_comp = 0.3 in the present study. The compacted band extends from the inner surface of the ring to a CF, as shown in figure 16(a). Ahead of the CF there exists a transitional zone with a finite thickness of w_CF (w_CF ~ 10d_p) across which ε and u_solid gradually decrease from ε_min and u_solid(R_CF) to ε ₀ and zero, respectively. To approximate the variations in ε and u_solid inside the transitional zone, a Gauss error function is employed

(4.10)\begin{gather}\varepsilon (r) = 1 - \frac{{{\varepsilon _{min}} + {\varepsilon _0}}}{2} + \frac{{{\varepsilon _{min}} - {\varepsilon _0}}}{2} \cdot \textrm{erf}\left( {4\frac{{r - {R_{CF}}}}{{{w_{CF}}}}} \right),\quad {R_{CF}} \le r \le \textrm{ }{R_{CF}} + {w_{CF}},\end{gather}

(4.11)\begin{gather}{u_{solid}}(r) = \frac{{{u_{solid,CF}}}}{2} - \frac{{{u_{solid,CF}}}}{2} \cdot \textrm{erf}\left( {4\frac{{r - {R_{CF}}}}{{{w_{CF}}}}} \right),\quad {R_{CF}} \le r \le \textrm{ }{R_{CF}} + {w_{CF}},\end{gather}

(4.12)\begin{gather}\textrm{erf}(x) = \frac{2}{{\sqrt {\rm \pi}}}\int_0^x {{\textrm{e}^{ - {\eta ^2}}}\,\textrm{d}\eta } ,\end{gather}

where u_solid,CF = u_solid(R_CF). The values of ε(r) and u_solid(r) described in (4.10) and (4.11) agree well with the numerically derived profiles of ε and u_solid as seen in figure 16(a). Thanks to the significant variation of the porosity, the flow dynamics of the interstitial gases undergoes fundamental changes when flowing through the transitional zone ahead of the CF. Thus, prior to resolving (4.7) it is necessary to determine the location of the CF, namely R_CF(t), which is to be addressed below.

The momentum balance of the annular compacted band demonstrated in figure 16(b) is given by (4.13)

(4.13) \begin{align} {\rho _{solid}}(1 - {\varepsilon _{min}})\int_{{R_{in}}}^{{R_{CF}}} {{{\dot{u}}_{solid}}(r)r\,\textrm{d}r} & =- {\rho _{solid}}(1 - {\varepsilon _0}){u_{solid,CF}}{V_{CF}}{R_{CF}}\notag\\ & \quad + {\rho _{solid}}(1 - {\varepsilon _{min}})\int_{{R_{in}}}^{{R_{CF}}} {{F_{\boldsymbol{\nabla }P}}r\,\textrm{d}r} \notag\\ & \quad + {\rho _{solid}}(1 - {\varepsilon _{min}})\int_{{R_{in}}}^{{R_{CF}}} {{F_{drag}}r\,\textrm{d}r} , \end{align}

where ${\dot{u}_{solid}}$ is the acceleration of the solid phase. The first term on the right-hand side of (4.13) accounts for the growing mass of the compacted band. The second and third terms on the right-hand side of (4.13) represent the total pressure gradient force and the total drag force exerted on the compacted band with a cross-section of unit area. The mass conservation of the incompressible compacted band requires that u_solid is proportional to 1/r

(4.14)\begin{equation}{u_{solid}}(r) = \frac{{{V_{in}}{R_{in}}}}{r},\quad {R_{in}} \le r \le {R_{CF}},\end{equation}

where V_in is equal to the particle velocity herein, V_in = u_solid(R_in). Differentiating (4.14) leads to the expression of ${\dot{u}_{solid}}$

(4.15)\begin{equation}{\dot{u}_{solid}}(r) = \frac{{{{\dot{V}}_{in}}{R_{in}}}}{r} + \frac{{V_{in}^2}}{r} - \frac{{{{({V_{in}}{R_{in}})}^2}}}{{{r^3}}},\quad {R_{in}} \le r \le {R_{CF}}.\end{equation}

At the CF

(4.16)\begin{gather}{u_{solid,CF}} = \frac{{{V_{in}}{R_{in}}}}{{{R_{CF}}}},\end{gather}

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

The V_CF and V_in ought to meet the Rankine–Huguenot condition

(4.18)\begin{equation}{V_{in}} = \frac{{{V_{CF}}{R_{CF}}}}{{{R_{in}}}}\left( {\frac{{{\varepsilon_0} - {\varepsilon_{min}}}}{{1 - {\varepsilon_{min}}}}} \right).\end{equation}

Equations (4.13)–(4.18) constitute the complete set of equations governing the shock compaction of the solid phase which is coupled with the evolution of the interstitial gas flows. Appendix E elaborates on the algorithm adopted to numerically solve those equations with the reference frame fixed on the inner surface of the ring.

4.2.3. Ring pulsation sub-model

Once the CF reaches the outer surface of the ring, the ring dynamics transitions to the ring pulsation synchronizing with the fluctuation of the central pressure. To predict t_ring, t_gas, h_dense,imp and R_dense,imp, we only need to reproduce the first expansion-to-implosion transition (EIT) of the ring. Since significant particle shedding occurs during the implosion, the ring before the first EIT is regarded as incompressible with an unvaried porosity, ε = ε_min. The momentum balance of the ring as an entirety is expressed as

(4.19) \begin{align} {\rho _{solid}}(1 - {\varepsilon _{min}})\int_{{R_{in}}}^{{R_{out}}} {r{{\dot{u}}_{solid}}\,\textrm{d}r} & = {\rho _{solid}}(1 - {\varepsilon _{min}})\int_{{R_{in}}}^{{R_{out}}} {{F_{\boldsymbol{\nabla }P}} \cdot r\,\textrm{d}r} \notag\\ & \quad + {\rho _{solid}}(1 - {\varepsilon _{min}})\int_{{R_{in}}}^{{R_{out}}} {{F_{drag}} \cdot r\,\textrm{d}r} . \end{align}

As in (4.13), the first and second terms on the right-hand side of (4.19) represent the total pressure gradient force and the total drag force exerted on the ring with a cross-section of unit area. The values of ${F_{{\nabla _P}}}$ and F_drag are described by (4.8) and (4.9), respectively.

Since the ring is incompressible, we have

(4.20)\begin{gather}\textrm{ }{u_{solid}}(r) = \frac{{{V_{in}}{R_{in}}}}{r},\quad {R_{in}} \le r \le {R_{out}},\end{gather}

(4.21)\begin{gather}{\dot{u}_{solid}}(r) = \frac{{{{\dot{V}}_{in}}{R_{in}}}}{r} + \frac{{V_{in}^2}}{r} - \frac{{{{({V_{in}}{R_{in}})}^2}}}{{{r^3}}},\quad {R_{in}} \le r \le {R_{out}}.\end{gather}

Integrating ${\dot{u}_{solid}}$ from the inner to the outer surface of the ring leads to

(4.22) \begin{align}\int_{{R_{in}}}^{{R_{out}}} {{{\dot{u}}_{solid}}r\,\textrm{d}r} &= \int_{{R_{in}}}^{{R_{out}}} {\left( {{{\dot{V}}_{in}}{R_{in}} + V_{in}^2 - V_{in}^2\frac{{R_{in}^2}}{{{r^2}}}} \right)\textrm{d}r}\notag\\ &= {\dot{V}_{in}}{R_{in}}({R_{out}} - {R_{in}}) + V_{in}^2\frac{{{{({R_{out}} - {R_{in}})}^2}}}{{{R_{out}}}}.\end{align}

Substituting (4.8), (4.9) and (4.22) into (4.19) reduces to

(4.23)\begin{align}({R_{out}} - {R_{in}}){\rho _{solid}}\left[ {{{\dot{V}}_{in}}{R_{in}} + V_{in}^2\left( {1 - \frac{{{R_{in}}}}{{{R_{out}}}}} \right)} \right] =- \frac{1}{{1 - {\varepsilon _{min}}}}\int_{{R_{in}}}^{{R_{out}}} {\textrm{d}{P_{gas}}/\textrm{d}r \cdot r\,\textrm{d}r} .\end{align}

Equations (4.20), (4.21) and (4.23) constitute the key governing equations of the ring dynamics in the ring pulsation sub-model. The evolution of the central gases and the gas filtration through the ring take place concurrently. The corresponding numerical algorithm is also presented in Appendix E.

4.3. Validation and applicability of the explosive dispersal model

Figure 17(a) compares the simulation derived and theoretically predicted trajectories of the inner and outer surfaces of the ring as well as the CF for the system 1035-50-1828. The consistencies between the numerical results and the theoretical predictions validate the capacity of the explosive dispersal model in terms of accounting for the ring dynamics. From figure 17(a) the theoretically predicted V_ring, h_dense,imp and R_dense,imp are identified, $V_{ring}^{theo} = 9.2\ \textrm{m}\ {\textrm{s}^{ - 1}}$, $h_{dense,imp}^{theo} = 0.045\ \textrm{m}$ and $R_{dense,imp}^{theo} = 0.22\ \textrm{m}$, which agree well with their numerically derived counterparts since the derivation rates are 15 %, 24 %, 14 %, respectively. As to the gaseous evolution, the theoretical model also provides validated predictions of the central pressure evolution, $P_{center}^{num}(t)$, as well as the cumulative gaseous mass loss ratio, $\chi _{gas}^{num}(t)$, as shown in figure 17(b,c). The predicted time scale of central gases, $t_{gas}^{theo} = 5.68\ \textrm{ms}$, has a deviation rate of 0.14 %. Since the ring cannot expand to twice its initial diameter, $\tau _{macro}^{theo}$ approaches infinity, while $\xi _{imp}^{theo} = 0.2$, indicating a failed dispersal, which is consistent with the simulation results where the system 1035-50-1828 indeed produces a failed dispersal.

Figure 17. Simulation derived and theoretically predicted results for the system 1035-50-1828. (a) Trajectories of the inner and outer surfaces of the ring as well as the CF. (b) Central pressure evolution. (c) Cumulative gaseous mass loss ratio.

Besides the macroscale parameters, the model properly accounts for the mesoscale momentum/energy transportation, especially during the shock compaction phase. As shown in figure 18(a,b), the model reproduced well the diffusional pressure field P_gas(r, t) and the transient velocity field of solid phase u_solid(r, t) that is observed in the simulation.

Figure 18. Comparisons of the simulation derived and theoretically predicted radial profiles of P_gas(r, t) (a) and u_solid(r, t) (b) during the shock compaction phase for the system 1035-50-1828.

Figure 19 plots the numerically derived and theoretically predicted time scales of the ring (figure 19a), $t_{ring}^{theo}$ and $t_{ring}^{num}$, as well as the time scales of the central gases (figure 19b), $t_{gas}^{theo}$ and $t_{gas}^{num}$, for all investigated systems. An overall consistency between $t_{ring}^{theo}$ and $t_{ring}^{num}$ is discernible, as is the case for $t_{gas}^{theo}$ and $t_{gas}^{num}$. Notably, the estimations of the time scale of the ring via the Gurney equation (3.10), $t_{ring}^{Gurney}$, are also supplemented in figure 19(a). Clearly, the consistency between $t_{ring}^{theo}$ and $t_{ring}^{num}$ is substantially improved compared with $t_{ring}^{Gurney}$, especially in the mid- and upper range of the (M/C)* wherein $t_{ring}^{theo}$ and $t_{ring}^{num}$ do not exist since the ring cannot expand to twice its initial radius prior to the EIT. The resulting $\tau _{macro}^{theo}$ for all systems are show in figure 19(c) in comparison with $\tau _{macro}^{num}$. All systems with $\tau _{macro}^{theo}\sim O({10^{ - 1}})$ indeed result in ideal dispersal since their actual $\tau _{macro}^{num}$ are of the order of O(10⁻¹). Thereby the theoretical model can successfully identify the ideal dispersal.

Figure 19. Comparisons of the simulation derived and theoretically predicted t_ring (a), t_gas (b), τ_macro (c) and ξ_imp (d). Inset in (d) shows the results from the validated failed dispersal.

For systems with the order of $\tau _{macro}^{theo}$ larger than O(10⁻¹), ξ_imp is used to decide whether they are the validated failed mode. The theoretically predicted $\xi _{imp}^{theo}$ for all systems are plotted in figure 19(d) in comparison with $\xi _{imp}^{num}$. For the ideal and partial dispersal modes, $\xi _{imp}^{num}$ does not exist since the rings disintegrate prior to the first implosion. By contrast, the ring implosion always occurs in the theoretical model since the ring disintegration is not accounted for, yielding a finite $\xi _{imp}^{theo}$ for all systems. The overall (M/C)* dependence of $\xi _{imp}^{theo}$ is similar to that of $\xi _{imp}^{num}$. It is most noteworthy that the systems denoted by the light red circular symbols, which are guaranteed to result in a validated failed dispersal, all dwell in the domain of $\xi _{imp}^{theo}\; > 0.1$ (see the inset of figure 19d), indicating the theoretical model suffices to identify the validated failed mode.