2.1. Governing equations and numerical algorithm

Numerical simulations were performed based on CCFD-DPM, a Euler–Lagrange (E-L) numerical approach suitable for modelling of gas–particle flows in laboratory-scale systems (Sundaresan, Ozel & Kolehmainen Reference Sundaresan, Ozel and Kolehmainen2018; Koneru, Rollin & Durant Reference Koneru, Rollin and Durant2020; Tian et al. Reference Tian, Zeng, Meng, Chen, Guo and Xue2020). The CCFD-DPM approach tracks and accounts for particle–particle contact interactions.

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{equation} \frac{\partial \left(\phi _{f}\left\langle \rho _{f}\right\rangle \right)}{\partial t}+\nabla \cdot \left(\phi _{f}\left\langle \rho _{f}\right\rangle \tilde{\boldsymbol{u}}_{f}\right)=0,\end{equation}

(2.2) \begin{equation} \begin{array}{l} \frac{\partial \left(\phi _{f}\tilde{\boldsymbol{u}}_{f}\left\langle \rho _{f}\right\rangle \right)}{\partial t}+\nabla \cdot \left(\phi _{f}\left\langle \rho _{f}\right\rangle \tilde{\boldsymbol{u}}_{f}\tilde{\boldsymbol{u}}_{f}+\phi _{f}\left\langle P_{f}\right\rangle \right)\\[8pt] =\left\langle P_{f}\right\rangle \nabla \phi _{f}+\sum _{i}\left\{\phi _{p,i}\rho _{p,i}D_{p,i}\left(\tilde{\boldsymbol{u}}_{p,i}-\left(\tilde{\boldsymbol{u}}_{f}\right)_{i}\right)\right\}, \end{array}\end{equation}

(2.3) \begin{equation} \begin{array}{l} \frac{\partial \left(\phi _{f}\left\langle \rho _{f}\right\rangle \tilde{E}_{f}\right)}{\partial t}+\nabla \cdot \left(\phi _{f}\left\langle \rho _{f}\right\rangle \tilde{E}_{f}\tilde{\boldsymbol{u}}_{f}+\phi _{f}\left\langle p_{f}\right\rangle \tilde{\boldsymbol{u}}_{f}\right)\\[10pt] =\left\langle P_{f}\right\rangle \nabla \phi _{f}\cdot \tilde{\boldsymbol{u}}_{p}+\sum _{i}\left\{\phi _{p,i}\rho _{p,i}D_{p,i}\left(\boldsymbol{u}_{p,i}-\left(\tilde{\boldsymbol{u}}_{f}\right)_{i}\right)\cdot \tilde{\boldsymbol{u}}_{p,i}\right\}. \end{array}\end{equation}

The gas volume fraction, density, velocity, pressure and 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 , respectively. In (2.1)–(2.3), $\langle \rangle$ and $\tilde{}$ denote phase-averaged and mass-averaged variables, respectively. Here ρ _p.i. and u _p,i are the density and velocity of particle i, D _p,i is the drag coefficient of particle i, and φ _p.i = w _i,f V _p,i /V _f is the contribution of particle i to the weighted particle volume fraction (w _i,f is the distributed weight that particle i contributes to the particle volume fraction in fluid cell f, and V _p,i and V _f are the volumes of particle i and the fluid cell).

The corresponding gas–particle coupling drag force is calculated using (2.4), and D _p,i can be expressed in the form of (2.5). Here sg is the specific weight of an individual particle with fluid, and r _p is the particle radius. Here C _d is the drag force coefficient, as shown in (2.6). This nonlinear drag force model is widely used in the numerical calculation of the interaction between shock waves and particle curtains and can consider the effects of the particle Reynolds number, Re _p , φ _p and Mach number, M _a . Thus, C _d is a function of Re _p , φ _p and M _a :

(2.4) \begin{equation} F_{Drag}=\frac{1}{8}\pi d_{p}^{2}C_{d}\left(Re_{p},M_{a},\phi _{p}\right)\rho _{f}\left| \boldsymbol{u}_{f}-\boldsymbol{u}_{p,i}\right| \left(\boldsymbol{u}_{f}-\boldsymbol{u}_{p,i}\right) \end{equation}

(2.5) \begin{equation} D_{p,i}=\frac{F_{Drag}}{m_{p}\left(\boldsymbol{u}_{f}-\boldsymbol{u}_{p,i}\right)}=\frac{3}{8sg}C_{d}\left(Re_{p},M_{a},\phi _{p}\right)\frac{\left| \boldsymbol{u}_{f}-\boldsymbol{u}_{p,i}\right| }{r_{p}},\end{equation}

(2.6) \begin{equation} C_{d}\left(Re_{p},M_{a},\phi _{p}\right)=C_{d,std}\left(Re_{p}\right)\xi _{1}\left(Re_{p},M_{a}\right)\xi _{2}\left(\phi _{p}\right).\end{equation}

The particle Reynolds and Mach numbers are defined as

(2.7) \begin{equation} {Re}_{p}=\frac{\rho _{f}\left| \boldsymbol{u}_{f}-\boldsymbol{u}_{p}\right| {d}_{p}}{\mu },\quad M_{a}=\frac{\left| \boldsymbol{u}_{f}-\boldsymbol{u}_{p}\right| }{c},\end{equation}

where μ is the dynamic viscosity of the fluid and c is the local speed of sound of the fluid. The drag force model introduces the effects of the Mach number and particle volume fraction by correcting the standard drag correlation C _d, which was developed by Clift & Gauvin (Reference Clift and Gauvin1970) using the Mach number multiplication factor (Parmar, Haselbacher & Balachandar Reference Parmar, Haselbacher and Balachandar2010) ξ ₁ and the particle volume fraction multiplication factor (Sangani, Zhang & Prosperetti Reference Sangani, Zhang and Prosperetti1991) ξ ₂, as follows:

(2.8) \begin{equation} C_{d,std}=\frac{24}{Re_{p}}\left(1+0.15Re_{p}^{0.687}\right)+0.42\left(1+\frac{42500}{Re_{p}^{1.16}}\right)^{-1},\end{equation}

(2.9) \begin{equation} \xi _{1}\left({Re}_{p},M_{a}\right)=C_{d}\left({Re}_{p},M_{a}\right)/C_{d,std}\left({Re}_{p}\right),\end{equation}

(2.10) \begin{equation} \xi _{2}\left(\phi _{p}\right)=\frac{1+2\phi _{p}}{\left(1-\phi _{p}\right)^{3}}.\end{equation}

The particle phase is represented by discrete particles whose motion is governed by Newton’s second law,

(2.11) \begin{equation} \frac{{d}\boldsymbol{u}_{p,i}}{{d}t}=D_{p,i}\left(\boldsymbol{u}_{f}-\boldsymbol{u}_{p,i}\right)-\frac{1}{\rho _{p}}\nabla \left\langle P_{f}\right\rangle +\frac{1}{m_{p}}\sum _{j}\boldsymbol{F}_{C,ij},\end{equation}

(2.12) \begin{equation} \frac{{d}\boldsymbol{x}_{p,i}}{{d}t}=\boldsymbol{u}_{p,i},\end{equation}

where u _p,i and x _p,i are the velocity and displacement of particle i, respectively, and F _C,ij represents the collision force between particles 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 particles are driven by the pressure gradient force, drag force and collision force between particles (2.13). A soft sphere model, represented by a linear-spring–dashpot model, was employed to model the collision force between particles. Hence, F _C,ij consists of a repulsive force and a damping force

(2.13) \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 the particles and δ _n and u _n,ij are the overlap and normal velocity differences between the particles in contact, respectively. Here γ _n,p is a function of the particle restitution coefficient ε _p ,

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

To solve the equations governing the gases, the weighted essentially non-oscillatory scheme (Liu Reference Liu1994) 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 (Toro Reference Toro2009). The third-order Runge–Kutta method was applied for time integration. The equations describing the particle 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 particles. Numerical details about CCFD-DPM can be found in our previous works (Tian et al. Reference Tian, Zeng, Meng, Chen, Guo and Xue2020).

2.2. Verification of the numerical method

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 investigated shock dissipation through particle curtains; experiments carried out by Britan et al. (Reference Britan, Ben-Dor, Igra and Shapiro2006), which assessed both the gaseous and solid pressures inside particle columns impinged head on by shocks; and experiments involving the shock-induced interfacial instability of granular media (Xue et al. Reference Xue, Du, Shi, Gan and Bai2018, Reference Xue, Shi, Zeng, Tian, Han, Li, Liu, Meng, Guo and Bai2020).

Since the 2-D configurations are used to investigate the shock attenuation effect in the present work, it is necessary to validate the 2-D CCFD-DPM modelling with the experiments and 3-D particle-resolved simulations involving shock propagation through particle layers. Specifically, Grinten et al. measured the pressure traces upstream of and inside a densely packed particle column placed against the close end of a shock tube (van der Grinten, van Dongen & van der Kogel Reference van der Grinten, van Dongen and van der Kogel1985). The schematic of the experimental set-up is shown figure 2(a). Figure 2(b) compares the pressure histories at different locations measured in the experiment (dashed line) and those derived from our 2-D simulations using a 2-D configuration presented in figure 2(a). A good agreement is evident in terms of the overall profiles of pressures.

Figure 2.

(a) Schematic of the experimental set-up in which a normal shock with the Mach number M _s = 1.155 head-on collides with a densely packed long particle column (φ _p = 70 %, L ₀ = 1895 mm, d _p = 375 $\unicode{x03BC}$ m) placed against the close end of the shock tube. Four pressure transducers (P_A –P_D ) were used to measure pressure traces upstream of and inside the particle column. (b) Comparison of pressure histories at P_A –P_D obtained in the experiment and our 2-D simulations.

As to the shock interaction between the mobile particle layers, Theofanous et al. (Reference Theofanous, Mitkin and Chang2016) and Wagner et al. ( Reference Wagner, Beresh, Kearney, Trott, Castaneda, Pruett and Baer2012) carried out a series of experiments of shock impingements of thin particle curtains, which serve as a definitive set of validation cases for the pertinent numerical simulations. We performed corresponding 2-D simulations and obtained the unsteady flows as well as the particle dispersion behaviours. Figures 3(a,b) and 3(d,e) compare the numerically derived pressure histories upstream and downstream of particle curtain with the experimental results reported by Theofanous et al. and Wagner et al., respectively. Here the pressure is scaled via the post incident shock pressure,

Figure 3.

(a,b,d,e) The numerically and experimentally derived pressure histories upstream and downstream of shock impinged thin particle curtains. (c,f) The numerically and experimentally derived trajectories of UI and DI of shock dispersed particle curtains. (a–c) Experimental results are extracted from the reference (Theofanous et al. Reference Theofanous, Mitkin and Chang2016) in which the main parameters of particle curtain are φ _p = 32 %, L ₀ = 23.2 mm, d _p ranging from 860 $\unicode{x03BC}$ m to 960 $\unicode{x03BC}$ m; (d–f) experimental results are extracted from the reference (Wagner et al. Reference Wagner, Beresh, Kearney, Trott, Castaneda, Pruett and Baer2012) in which the main parameters of particle curtain are φ _p = 21 %, L ₀ = 1.5 mm, d _p ranging from 106 $\unicode{x03BC}$ m to 125 $\unicode{x03BC}$ m. Here TS and RS represent the transmitted shock in the particle layers and the shock reflected off the upstream interface.

(2.15) \begin{equation} P^{*}=\frac{P-P_{1}}{P_{2}-P_{1}},\end{equation}

and the time is scaled via the time the incident shock takes to travel the length of particle bed, L ₀:

(2.16) \begin{equation} \tau =\frac{t}{L_{0}/u_{s}}.\end{equation}

As seen in figure 3, our simulations are able to accurately capture the arrival times of the reflected shock (RS) and transmitted shock (TS) and provide pressure profiles consistent with the experimental measurements. Besides, our simulations also reproduce the same dispersion rates of particle curtain in terms of the trajectories of upstream and downstream interfaces (UI and DI).

Particle resolved direct numerical simulations are able to depict the overall shock-induced flows with particle scale details rather than pressure histories at several locations. Figure 4 shows the evolutions of streamwise mean pressure (figure 4 a,c) and velocity (figures 4 b and 4 d) profiles attained from 3-D particle-resolved viscous (Osnes et al. Reference Osnes, Vartdal, Omang and Reif2020) (figure 4(a,b) and inviscid simulations (Mehta et al. Reference Mehta, Neal, Salari, Jackson, Balachandar and Thakur2018) (figure 4 c,d). Here the velocity is scaled via post incident shock velocity, u ₂:

Figure 4.

Comparisons of the streamwise mean pressure and flow velocity profiles, P ^*(x/L ₀) (a), U ^*(x/ L ₀) (c), P ^*(x/d _p ) (b) and U ^*(x/d _p ) (d), at different times derived from the present 2-D (solid line), 3-D (dotted line) CCFD-DPM simulations and 3-D viscous, (a)–(b) (Osnes et al. Reference Osnes, Vartdal, Omang and Reif2020) and inviscid, (c)–(d) (Mehta et al. Reference Mehta, Neal, Salari, Jackson, Balachandar and Thakur2018) particle-resolved simulations (dashed line).

(2.17) \begin{equation} U^{*}=\frac{u}{u_{2}}.\end{equation}

In the viscous simulation (Osnes et al. Reference Osnes, Vartdal, Omang and Reif2020), a particle curtain with 30 layers of particles (φ _p = 10 %, L ₀/ d _p = 30) consisting of immobile particles (d _p = 34 $\unicode{x03BC}$ m) are impinged by a normal incident shock with M _s = 2.6. In the inviscid simulation, a stronger shock (M _s = 3) head-on impinges a thinner particle curtain (L ₀/ d _p = 19) consisting of larger mobile particles (d _p = 100 μm) while the initial volume fraction is kept consistent, φ _p = 10 %.

The particle resolved simulations have the advantages of capturing the particle scale flow fluctuations, especially during the initial shock induced transient. A previous study (Ling et al. Reference Ling, Wagner, Beresh, Kearney and Balachandar2012) found that a sustained variation in flow characteristics (velocity, pressure, density) across the particle layers persists after the passage of transmitted shock. Also the length scale of this variation is of the order of the layer's thickness. Therefore, particle unresolved simulations using the particle-scale volume averaged homogeneous formulations of drag force suffice to adequately account for the streamwise variation of mean flows in the wake of transmitted shock. The mean pressure and velocity profiles derived from the 2-D (solid curves) and 3-D (dotted curves) CCFD-DPM simulations are superimposed in figure 4 in which the structural parameters of numerical configuration are identical to those used in the particle resolved simulations. Evidently, the advancing of the RS and TS as well as the profiles of pressure and velocity derived from either 2-D or 3-D CCFD-DPM simulations are consistent with those derived from the particle resolved simulations, further corroborating the validity of the current 2-D modelling in describing the shock induced unsteady flows upstream of, inside and downstream of dense particle layers.

Finally it is necessary to justify the capacity of the particle unresolved simulations in accounting for the shock attenuation effect by elucidating the partition of energy influx caused by the TS in particle resolved and unresolved methods. The primary difference in this regard arises from the origin and resolution of the attained particle-scale fluctuations, which are an integral part of the total kinetic energy (Mehta, Jackson & Balachandar Reference Mehta, Jackson and Balachandar2019a ). The particle-resolved method accounts for the particle-scale fluctuations arising from multiple factors, including the shock phenomenon at individual particles, interparticle interactions via particle wakes and viscous effects (Regele et al. Reference Regele, Rabinovitch, Colonius and Blanquart2014; Osnes et al. Reference Osnes, Vartdal, Omang and Reif2019, Reference Osnes, Vartdal, Omang and Reif2020; Mehta et al. Reference Mehta, Jackson and Balachandar2019a ,Reference Mehta, Neal, Salari, Jackson, Balachandar and Thakur b , Reference Mehta, Salari, Jackson and Balachandar2018). By contrast, the particle-unresolved method using Ergun-type drag models considers the particle-scale fluctuations through the variations in the local flow parameters (Re), particle size (d _p ) and local volume fraction (φ _p ) (Ling et al. Reference Ling, Wagner, Beresh, Kearney and Balachandar2012; Houim & Oran Reference Houim and Oran2016; Theofanous & Chang Reference Theofanous and Chang2017; Meng et al. Reference Meng, Zeng, Tian, Li, He and Guo2019; Sugiyama et al. Reference Sugiyama, Ando, Shimura and Matsuo2019). The resolution of the particle-scale fluctuations observed in unresolved simulations depends on the ratio between the grid length and the particle diameter as well as the algorithm related to the information exchange between the fluid and particle phases (Regele et al. Reference Regele, Rabinovitch, Colonius and Blanquart2014; Mehta et al. Reference Mehta, Neal, Salari, Jackson, Balachandar and Thakur2018). Here Re, d _p and φ _p are the dominating parameters of the unsteady flows after the initial shock-induced transient because the steady forces on the particles as a function of these variables become predominant shortly at this time (Sun et al. Reference Sun, Saito, Takayama and Tanno2005; Parmar et al. Reference Parmar, Haselbacher and Balachandar2011, Reference Parmar, Haselbacher and Balachandar2012). Although the particle-unresolved methods may underestimate the particle-scale fluctuations during the initial shock-induced transient due to the lack of unsteady forces, the shock attenuation can be adequately accounted for because it is associated primarily with the long-term unsteady flows after the passage of the transmitted shock. The study by Osnes et al. (Reference Osnes, Vartdal, Omang and Reif2019, Reference Osnes, Vartdal, Omang and Reif2020) revealed that the initial shock-induced transient is negligible compared with the time scale of shock attenuation.

Figure 5.

Computational configuration of a planar shock impinging head-on the UI of a fixed particle bed with a length of L ₀. The ambient condition is imposed upon the DI of the particle bed.

Figure 6.

Distributions of numerical cases in the parameter spaces (M _s , φ _p ) (a), ( $\mathit{L}_{0}^{\mathrm{*}}$ , d _p ) (b) and (Re _p , θ) (c).

2.3. Numerical set-up

The objective of this work is to obtain a fundamental understanding of the shock attenuation effect through granular media by proposing scaling laws for the upstream and downstream pressures rather than by precisely resolving these pressures. Hence, we consider the 2-D configuration shown in figure 5. A fixed particle bed consisting of randomly positioned/packed cylindrical particles is placed downstream of a planar incident shock (IS). Regele et al. noted that a proper 2-D numerical configuration can qualitatively capture the multidimensional behaviour of a real 3-D particle cloud impinged by an incident shock (Regele et al. Reference Regele, Rabinovitch, Colonius and Blanquart2014). The manageable cost of 2-D simulations allows us to investigate shock attenuation in a high-dimensional parameter space with variables spanning several decades. The computational domain upstream of the incident shock wave is long enough to ensure that no upstream disturbances influence the flows inside the particle bed during the computation duration. To minimize the variability of the streamwise mean flow parameters due to spanwise particle-scale fluctuations, the cross-sectional length of the fixed particle bed H is large enough, and periodic boundary conditions are used on the spanwise boundaries (y = $\pm$ H/2). Specifically, H ranges from 1100–110 d _p as φ _p increases from 0.01–0.5.

As indicated in the experimental and numerical studies (Britan et al. Reference Britan, Ben-Dor, Igra and Shapiro2001; Bakken et al. Reference Bakken, Slungaard, Engebretsen and Christensen2003; Osnes et al. Reference Osnes, Vartdal, Omang and Reif2019; Mehta et al. Reference Mehta, Salari, Jackson and Balachandar2019b ), the shock strength M _s , the length of the particle bed $\mathit{L}_{0}^{\mathrm{*}}$ = L ₀/d _p , the particle diameter d _p , and the volume fraction φ _p are the primary factors controlling shock attenuation. Hence, we devised 126 numerical cases with vastly varying M _s , $\mathit{L}_{0}^{\mathrm{*}}$ , d _p and φ _p whose locations in the parameter spaces (M _s , φ _p ) and ( $\mathit{L}_{0}^{\mathrm{*}}$ , φ _p ) are shown in figures 6(a) and 6(b), respectively. Specifically, M _s varies from 1.2–2.58; accordingly, the postincident shock flow varies from subsonic (M ₂ = 0.29) to supersonic (M ₂ = 1.23). The state of the particle layers ranges from moderately concentrated (φ _p = 0.01) to dense (0.01< φ _p <0.5) particle clouds and then to randomly close packed beds (φ _p ≥ 0.5). The mean values of the particle diameter, which obey a normal distribution (Δd _p /d _p = ±5.5 %), span three orders of magnitude and range from 1–910 μm. The length of the particle beds varies from dozens of particle layers ( $\mathit{L}_{0}^{\mathrm{*}}$ = 25) to semi-infinite layers ( $\mathit{L}_{0}^{\mathrm{*}}$ = 2200). The exact parameter values of each numerical case are listed in table 2 in Appendix A. The distributions of the numerical cases in the non-dimensional parameter space (Re _p , θ) are plotted in figure 6 c). The θ ranges from O(10⁻¹) to O(10³), which is consistent with the range of θ variation observed in the previous experiments (figure 1). Thereafter, the numerical cases are denoted by the combination of four critical variables, i.e. C-M _s -φ _p - $\mathit{L}_{0}^{\mathrm{*}}$ -d _p (in μm).

3.1. Transmission-dominated shock attenuation

Figure 7(a) shows a typical x–t diagram of the mean pressures P ^* averaged over the transverse direction in a moderately concentrated particle bed consisting of a handful of particle layers (1.67-0.01-110-910). The corresponding unsteady flows are dominated by various waves discerned from the numerical schlieren (figure 7 b), which is consistent with other investigations (Levy et al. Reference Levy, Ben-Dor, Skews and Sorek1993; Rogue et al. Reference Rogue, Rodriguez, Haas and Saurel1998; Ling et al. Reference Ling, Wagner, Beresh, Kearney and Balachandar2012; Regele et al. Reference Regele, Rabinovitch, Colonius and Blanquart2014). As the IS impinges on the UI of the particle bed, a part of it is transmitted into the bed as a TS, and another part is reflected, becoming the primary reflected shock (PRS). As the TS penetrates layers of particles, it bounces back and forth inside the interparticle pores, some of which find a way to emerge from the UI as weak compression waves (CWs). The CWs eventually catch up with the PRS, which hence becomes increasingly stronger. Moreover, the TS becomes increasingly dispersed with decreasing strength. When the TS emerges from the downstream interface (DI) of the particle bed at t = t _TS , a downstream travelling contact surface (CS) is created, trailed by a later emerging contact surface (CS₁) previously formed at the UI. Simultaneously, a reflected rarefaction wave (RRW) travels upstream. A steady pressure gradient is set up upon the RRW reaching the UI at t = t _RRW.

Figure 7.

(a) The x–t diagram of mean pressures P ^*, (b) numerical schlieren, evolution of the P ^* profiles in (c) the early and (d) late times, (e) temporal variations in pressures upstream ( $P_{{UI}}^{\mathrm{*}}$ , X = 0) of, inside (P ^*(X = 0.5)) and downstream ( $P_{{DI}}^{\mathrm{*}}$ , X = 1) of particle bed for the case C-1.67-0.01-110–910. Note that the induced density ρ is scaled by ρ ₂. The dashed curve in (d) represents the steady P ^* profile predicted by the Morrison model with Δ $P_{{s}}^{\mathrm{*}}$ = $P_{{s},{UI}}^{\mathrm{*}}$ - $P_{{s}, {DI}}^{{*}}$ . Here PRS, primary reflected shock; CWs, compression waves; RRW, reflected rarefaction wave; CS₁, contact surface formed at the UI; CS₂, contact surface formed at the DI.

The evolution of the P ^* profiles shown in figure 7(c) manifests the decay of the TS as well as the gradual increase in the upstream pressure due to the arrival of CWs. Specifically, the strength of the TS, $P_{TS}^{\mathrm{*}}$ , is defined as the pressure jump across its front, which decreases by 12 % after the TS penetrates 110 layers of particles while the propagation velocity remains largely constant. After t_TS , in the wake of RRW, the profile of P ^*(X) converges to a steady $P_{{steady}}^{\mathrm{*}}$ (X), corresponding to a steady diffusional pressure field dictated by the pressure differential between the upstream and downstream directions, Δ $P_{s}^{{*}}$ = $P_{{s},{UI}}^{\mathrm{*}}$ - $P_{{s}, {DI}}^{{*}}$ (see figure 7d ). Note that $\mathit{P}_{{steady}}^{\mathrm{*}}$ (X) agrees well with the profile predicted by the Morrison model (Morrison Reference Morrison1972) with Δ $\mathit{P}_{{s}}^{\mathrm{*}}$ = $\mathit{P}_{{s},{UI}}^{{*}}$ - $\mathit{P}_{{s}, {DI}}^{{*}}$ . As shown in figure 7(e), $\mathit{P}_{{DI}}^{\mathrm{*}}$ jumps to $\mathit{P}_{{s}, {DI}}^{{*}}$ immediately after the TS emerges from the DI, $\mathit{P}_{{s}, {DI}}^{{*}}$ = 0.88, while $\mathit{P}_{{UI}}^{\mathrm{*}}$ gradually increases and converges to $\mathit{P}_{{s},{UI}}^{\mathrm{*}}$ at t = t_RRS , $\mathit{P}_{{s},{UI}}^{\mathrm{*}}$ = 1.72. In this scenario, the pressure drop in terms of Δ $\mathit{P}_{{s}}^{\mathrm{*}}$ clearly depends on the propagation of the TS, which is referred to as the transmission-dominated mode hereafter.

3.2. Transmission-filtration-dominated shock attenuation

For a denser and longer particle bed, the TS may not emerge from the UI of the bed, as indicated in figure 8. Although the x–t diagram of the streamwise mean pressure P ^* (figure 8 a) shows a pattern similar to that in figure 7(a), no RRW can be identified from the numerical schlieren (figure 8 b), indicating complete dissipation of the TS midway through the bed. The evolution of the profiles of P ^*(X) (figure 8c) demonstrates a rapid decay of the TS, which degrades into a compression wave when it travels to X = X _TS = 0.67; therefore, no pressure jump across the TS front is discernible. Hereafter, the distance from the UI upon which the TS degrades into a transmitted compression wave (TCW) is referred to as the maximum transmission length, $\mathit{L}_{{TS}}^{\mathrm{*}}$ , which is a function of M _s , φ _p and d _p, as revealed below. The determination procedure of the $\mathit{L}_{{TS}}^{\mathrm{*}}$ is presented in Appendix B.

Figure 8.

(a) The x–t diagram of mean pressures P ^*, (b) numerical schlieren, evolution of the P ^* profiles in in (c) the early and (d) times, (e) temporal variations in pressures upstream ( $\mathit{P}_{{UI}}^{\mathrm{*}}$ , X = 0) of, inside (P ^*(X = 0.5)) and downstream ( $\mathit{P}_{{DI}}^{\mathrm{*}}$ , X = 1) of particle bed for the case C-1.67-0.01-1100-1. The dashed and dotted curves in (d) represent the steady pressure profiles predicted by the Morrison model with Δ $\mathit{P}_{{s}}^{\mathrm{*}}$ = $\mathit{P}_{5}^{\mathrm{*}}$ – $\mathit{P}_{1}^{\mathrm{*}}$ , Δ $\mathit{P}_{{s}}^{\mathrm{*}}$ = $\mathit{P}_{{s},{UI}}^{\mathrm{*}}$ – $\mathit{P}_{{s}, {DI}}^{{*}}$ , respectively.

As the TCW emerges from the DI, the downstream pressure $\mathit{P}_{{DI}}^{\mathrm{*}}$ slowly increases until the arrival of CS₁ at τ = $\tau _{{\textit{CS}_{1}}}$ (figure 8 e). The $\mathit{P}_{{DI}}^{\mathrm{*}}$ eventually stabilizes at a value appreciably higher than the ambient pressure, $\mathit{P}_{{s}, {DI}}^{{*}}$ = 0.35, contradicting the previously held assumption that $\mathit{P}_{{s}, {DI}}^{{*}}$ remains ambient if no TS emerges. As described in Appendix C, the noticeable increase in $\mathit{P}_{{DI}}^{\mathrm{*}}$ is associated with the thermodynamic evolution of heated gases between CS₁ and CS₂ flowing through the DI. Moreover, the upstream pressure $\mathit{P}_{{UI}}^{\mathrm{*}}$ (τ) (figure 8 e) undergoes a sharp jump followed by a gradual increase, converging to $\mathit{P}_{{s},{UI}}^{\mathrm{*}}$ = 2.96 at late times. Notably, the increase in $\mathit{P}_{{UI}}^{\mathrm{*}}$ does not completely cease when the last CW (CW _T ) starting from X _TS reaches the UI at τ _CWT = 1.95. Rather, we observe a persistent increase in $\mathit{P}_{{DI}}^{\mathrm{*}}$ , which lasts until CS₁ emerges from the DI at τ _CS1 = 6.6. The fast-flowing interstitial gases swept by the CW _T continue to be slowed down by friction between gases and the solid skeleton, thereby generating compression waves that can also emerge from the UI and further increase $\mathit{P}_{{DI}}^{\mathrm{*}}$ . The overall pressure profile becomes steady upon τ _CS1, $\mathit{P}_{{s}}^{\mathrm{*}}(\mathit{X})$ , which is consistent with Morrison’s prediction for $\mathit{P}_{{s},{UI}}^{\mathrm{*}}$ = 2.96 and $\mathit{P}_{{s}, {DI}}^{{*}}$ = 0.35, as shown in figure 8(d). Since the pressure drop in this scenario depends on both the TS and the gas filtration, this type of shock attenuation is referred to as the transmission-filtration-dominated mode.

3.3. Filtration-dominated shock attenuation

The TS can barely survive a few layers of particles when the IS penetrates a densely packed particle bed, as observed in the case of C-1.67-0.5-110-910. Hence, gas filtration takes over and becomes the dominant process. A typical diffusional pressure field develops alongside the gas filtration through the bed, as shown in figure 9(a), and a steady flow field forms after CS₁ emerges from the DI (figure 9 b). The friction between gases and the rigid skeleton slows down the induced flows (figure 9c), resulting in a gradual build-up of gaseous pressure (figure 9 d). Moreover, the upstream travelling compression waves caused by the slowing down of induced flows give rise to a slow but persistent increase in $\mathit{P}_{{UI}}^{\mathrm{*}}$ after the sharp jump at the PRS (figure 9 g). Eventually, $\mathit{P}_{{UI}}^{\mathrm{*}}$ converges to a $\mathit{P}_{{s},{UI}}^{\mathrm{*}}$ slightly lower than the pressure after the shock wave is reflected off a solid surface known as $\mathit{P}_{5}^{\mathrm{*}}$ . In the case of C-1.67-0.5-110-910, $\mathit{P}_{{s},{UI}}^{\mathrm{*}}$ = 3.29 while $\mathit{P}_{5}^{\mathrm{*}}$ = 3.38.

Figure 9.

(a) The x–t diagram of mean pressures P ^*, (b) numerical schlieren, (c–f) profiles of (c) V ^*(X), (d) P ^*(X), (e) T ^*(X), (f) ρ ^*(X), (g) temporal variations in pressures upstream ( $\mathit{P}_{{UI}}^{\mathrm{*}}$ , X = 0), inside (P ^*(X = 0.5)) and downstream ( $\mathit{P}_{{DI}}^{\mathrm{*}}$ , X = 1) of the particle bed for the case C-1.67-0.5-110-910. The dashed and dotted curves in (d) represent the steady pressure profiles predicted by the Morrison model with Δ $\mathit{P}_{{s}}^{\mathrm{*}}$ = $\mathit{P}_{5}^{\mathrm{*}}$ – $\mathit{P}_{1}^{\mathrm{*}}$ and Δ $\mathit{P}_{{s}}^{\mathrm{*}}$ = $\mathit{P}_{{s},{UI}}^{\mathrm{*}}$ – $\mathit{P}_{{s}, {DI}}^{{*}}$ , respectively. The temperature and velocity are scaled by T ₂ and u ₂, respectively.

In contrast with the isothermal filtration assumption of the conventional filtration model, shock-induced filtration may exhibit pronounced non-isothermal effects, as explored in our previous work (Li, Zeng & Xue Reference Li, Zeng and Xue2023). The scaled temperature profiles, T ^*(X), shown in figure 9(e) indicate an intensely heated gas band confined in CS₁ and CS₂ traversing the bed, resulting in a non-monotonic variation in the gaseous density, ρ ^*(X) (figure 9 f). The downstream gas density, $\rho _{{DI}}^{\mathrm{*}}$ , is also affected by the nozzling effect at the DI characterized by sudden acceleration as it approaches the DI (figure 9 f). As the heated gas band flows through the DI, the downstream pressure, $\mathit{P}_{{DI}}^{\mathrm{*}}$ , slowly increases to a steady value, $\mathit{P}_{{s}, {DI}}^{{*}}$ = 0.11 (figure 9 g). The evolution of $\mathit{P}_{{DI}}^{\mathrm{*}}$ due to the combination of the non-isothermal effect and the nozzling effect is elucidated in Appendix C. The diffusional pressure profile predicted by the Morrison model with $\mathit{P}_{{s},{UI}}^{\mathrm{*}}$ = 3.29 and $\mathit{P}_{{s}, {DI}}^{{*}}$ = 0.11 provides better agreement with the simulated result than the theoretical estimation with $\mathit{P}_{{s},{UI}}^{\mathrm{*}} = \mathit{P}_{5}^{\mathrm{*}}$ and $\mathit{P}_{{s}, {DI}}^{{*}} = \mathit{P}_{1}^{\mathrm{*}}$ = 0 (figure 9 d). Since $\mathit{P}_{{s},{UI}}^{\mathrm{*}}$ and $\mathit{P}_{{s}, {DI}}^{{*}}$ are associated with non-isothermal gas filtration, the shock wave attenuation in this scenario is referred to as the filtration-dominated mode.

3.4. Mode transition of shock wave attenuation

As mentioned in the previous sections, the mechanism governing shock attenuation undergoes a fundamental transition when the filtration process takes over the TS propagation. Accordingly, the shock attenuation effects in terms of the pressure drop, ΔP ^* = $\mathit{P}_{{s},{UI}}^{\mathrm{*}}$ - $\mathit{P}_{{s}, {DI}}^{{*}}$ , ought to experience the mode transition. The most pronounced characteristic of the mode transition is the ratio between the maximum transmission length, $\mathit{L}_{{TS}}^{\mathrm{*}}$ , and the length of the particle bed, $\mathit{L}_{0}^{\mathrm{*}}$ , where χ = $\mathit{L}_{{TS}}^{\mathrm{*}}/\mathit{L}_{0}^{\mathrm{*}}$ . A transmission-dominated mode (Mode I) requires χ to be larger than one unit, χ >1, while the χ of the order of O(10⁻²) and below guarantees a filtration-dominated mode (Mode III). When χ is in between these values, a transmission-filtration-dominated mode (Mode II) prevails.

To accurately determine χ and the resulting shock attenuation mode for a specific combination of an IS and a particle bed, it is necessary to derive the corresponding $\mathit{L}_{{TS}}^{\mathrm{*}}$ . Figures 10(a) and 10(b) show the dependences of $\mathit{L}_{{TS}}^{\mathrm{*}}$ on φ _p , where d _p spans four orders of magnitude (figure 10 a) and the strength of the IS increases from weak (M _s = 1.2) to strong (M _s = 2.58) (figure 10 b). Here $\mathit{L}_{{TS}}^{\mathrm{*}}$ decreases dramatically with increasing φ _p towards the lower limit of the φ _p range. We observe a one-order-of-magnitude decrease in $\mathit{L}_{{TS}}^{\mathrm{*}}$ when φ _p only increases by several per cent from the lower limit of the dense regime, φ _p = 0.01. The decrease of $\mathit{L}_{{TS}}^{\mathrm{*}}$ is significantly slowed down when φ _p exceeds 0.1. The $\mathit{L}_{{TS}}^{\mathrm{*}}$ converges to several particle layers when the particle bed approaches densely packed packing (φ _p > 0.4), which is set by the resolution limit of the particle unresolved method. Although similar decaying dependences of $\mathit{L}_{{TS}}^{\mathrm{*}}$ on φ _p evident despite the distinctively varying M _s and d _p , the decay rate of $\mathit{L}_{{TS}}^{\mathrm{*}}$ is most significant for particle beds consisting of small beads (d _p = 1 μm) subject to a weak IS (M _s = 1.2). We note that the d _p effect becomes trivial as d _p exceeds 100 μm.

Figure 10.

Variations in the $\mathit{L}_{{TS}}^{\mathrm{*}}$ with the increasing φ _p for particle beds with different d _p (a) and subjected to the IS with different M _s (b). The shaded bands in (a) encompass the results attained from cases with the same d _p but different M _s (M _s ∼ 1.2–2.58) while the shaded bands in (b) envelop the results attained from cases with the same M _s but different d _p (d _p ∼ 1–910 $ \unicode{x03BC}$ m).

In previous studies on the transmission-dominated mode of shock attenuation (Osnes et al. Reference Osnes, Vartdal, Omang and Reif2019, Reference Osnes, Vartdal, Omang and Reif2020; Mehta et al. Reference Mehta, Salari, Jackson and Balachandar2019b ), the collective effects of the M _s and d _p shown in figure 9 were incorporated into the particle Reynolds number Re _p . These studies reported a faster decay of TS with an increased Re _p (Osnes et al. Reference Osnes, Vartdal, Omang and Reif2020). The product of the first two terms in the expression of θ (1.1), which is denoted as θ ₁ as follows, explicitly accounts for the attenuation effects from φ _p as well as Re _p:

(3.1) \begin{equation} \theta _{1}=f_{m}\, \frac{\phi _{p}}{1-\phi _{p}}=\left(\frac{\alpha \phi _{p}}{{Re}_{p}}+\beta \right)\, \frac{\phi _{p}}{1-\phi _{p}}.\end{equation}

Figure 10(a) plots the numerically derived $\mathit{L}_{{TS}}^{\mathrm{*}}$ against θ ₁ with α = 150 and β = 1 in double logarithmic coordinates. Clearly, θ ₁ fails to scale the $\mathit{L}_{{TS}}^{\mathrm{*}}$ attained from a wide range of combinations M _s , φ _p and d _p, as indicated by considerable variabilities.

Another dimensionless variable that has the potential to scale $\mathit{L}_{{TS}}^{\mathrm{*}}$ is related to the filtration Reynolds number first introduced by Morrison (Reference Morrison1972),

(3.2) \begin{equation} {Re}_{f}=\frac{b}{a^{2}}\frac{P_{5}-P{1}^{}{}{}}{P_{1}}\, \frac{P_{1}\rho _{1}}{L_{0}\mu _{1}^{2}},\end{equation}

where subscripts 1 and 5 refer to the parameters of the downstream ambient gases and gases behind the reflected shock, respectively, and the Darcy, a, and Forchheimer, b, coefficients are functions of φ _p and d _p . For particle beds composed of spherical particles, a and b are expressed as

(3.3) \begin{equation} a = 180\frac{\phi _{p}^{2}}{\left(1-\phi _{p}\right)^{2}}\frac{1}{d_{p}^{2}},\end{equation}

(3.4) \begin{equation} b=1.8\frac{\phi _{p}}{1-\phi _{p}}\frac{1}{{d}_{p}}.\end{equation}

Substituting (3.3) and (3.4) into (3.2), we have

(3.5) \begin{equation} {Re}_{f}=\frac{1}{18000L_{0}^{*}}\, \frac{\left(P_{5}-P_{1}\right)}{P_{1}}\, \frac{\rho _{1}P_{1}d_{p}^{2}}{\mu _{1}^{2}}\, \frac{\left(1-\phi _{p}\right)^{3}}{\phi _{p}^{3}}.\end{equation}

Here Re _f is found to underlie the transition of filtration-induced gas flow from viscous loss-dominated Darcy flow (Re _f →0) to inertial loss-dominated Forchheimer flow (Re _f →∞) (Britan, Shapiro & Ben-Dor Reference Britan, Shapiro and Ben-Dor2007). Hence, Re _f is indicative of the nature of interphase friction, which is the predominant mechanism for the filtration-dominated mode of shock attenuation. Leaving out the first term involving $\mathit{L}_{0}^{\mathrm{*}}$ , the remainder of Re _f constitutes a new dimensionless parameter, Re _f,1, which manifests the effects arising from M _s (correlated with P ₅), d _p and φ _p:

(3.6) \begin{equation} {Re}_{f,1}=\frac{\left(P_{5}-P_{1}\right)}{P_{1}}\, \frac{\rho _{1}P_{1}d_{p}^{2}}{\mu _{1}^{2}}\, \frac{\left(1-\phi _{p}\right)^{3}}{\phi _{p}^{3}}.\end{equation}

The plot of $\mathit{L}_{{TS}}^{\mathrm{*}}$ against Re _f,1 in the double logarithmic coordinates is shown in figure 11(b), revealing a significant improvement in data convergence. The relative root mean square error (RRMSE) of all the numerically derived $\mathit{L}_{{TS}}^{\mathrm{*}}$ values from the linearly fitted master curve denoted by the black line in figure 11(b) is 5.6 %.

Figure 11.

Variations in the $\mathit{L}_{{TS}}^{\mathrm{*}}$ with the increasing θ ₁ (a) and Re _f,1 (b).

Inspired by θ ₁ and Re _f,1, which correlate with the attenuation mechanisms in the transmission- and filtration-dominated modes, respectively, we construct a dimensionless parameter composed of Re _p and a power function of φ _p as

(3.7) \begin{equation} \eta ={Re}_{p}\, \frac{\left(1-\phi _{p}\right)^{3}}{\phi _{p}^{4}}.\end{equation}

The correlations between η and θ ₁ as well as Re _f,1 are shown in figures 12(a) and 12(b), respectively. Clearly, there exists an explicit semilinear Re _f,1 dependence of η in contrast with the scattering η corresponding to a given θ ₁.

Figure 12.

Variations in the η with the increasing θ ₁ (a) and Re _f,1 (b).

The variation in $\mathit{L}_{{TS}}^{\mathrm{*}}$ with η is shown in figure 13 in the double logarithmic coordinates. The majority of the data collapse into a narrow band centred on a linearly fitted master curve, yielding a satisfactorily small RRMSE of 4.1 %. Therefore, the master curve in figure 13 represents the scaling law for $\mathit{L}_{{TS}}^{\mathrm{*}}$ :

Figure 13.

Variations in the $\mathit{L}_{{TS}}^{\mathrm{*}}$ with the increasing η.

(3.8) \begin{equation} \lg L_{{TS}}^{*}=0.29\lg \eta +0.22.\end{equation}

When η approaches zero, which occurs with extremely weak IS (Re _p →0) and/or exceedingly small d _p (Re _p →0) and/or the most densely packed bed (φ _p →1), $\mathit{L}_{{TS}}^{\mathrm{*}}$ converges to its lower limit, which is of the order of O(10⁰). Note that φ _p in (3.7) is for the 2-D system. Thus, for the 3-D bed, we should change φ_p,3-D into its 2-D counterpart, φ_p,2-D , using φ_p,2-D ≈ 1.22φ_p,3-D .

Using the scaling law of $\mathit{L}_{{TS}}^{\mathrm{*}}$ (3.8), we can calculate the ratio between $\mathit{L}_{{TS}}^{\mathrm{*}}$ and $\mathit{L}_{0}^{\mathrm{*}}$ , i.e. χ, for all numerically and experimentally investigated shock attenuation systems, including the cases plotted in figure 1. An explicit correlation between χ and θ is shown in figure 14. As θ ranges from O(10⁻¹) to O(10³), χ varies from O(10³) to O(10⁻²), whereby the shock attenuation crosses over all three modes. Notably, the χ values of systems whose θ values are very close vary enough so that they become quite distinguishable. Specifically, the χ of systems with θ ∼ 3 highlighted in figure 1 varies from 3.5 to 6.5. It is worth noting that whenever the results from 3-D configurations are compared with their 2-D counterparts, or vice versa, the 3-D volume fraction, φ_p,3-D , ought to be converted into 2-D equivalent, φ_p,2-D . Following the conversion principle proposed by Borchard-Ott (Reference Borchard-Ott2011), we adopt a the linear interpolation between the porosity at the maximum (hexagonal) and minimum (none particle) packing density states for monosized packing with disks (2-D) and spheres (3-D):

Figure 14.

Variations in the θ with the increasing χ.

(3.9) \begin{equation} \phi _{{p},3D}=\phi _{{p},2D}\frac{\phi _{{p},3D}^{\textit{HCP}}-0}{\phi _{{p},2D}^{\textit{HCP}}-0}=\frac{0.74-0}{0.097-0}\phi _{{p},2D}=0.82\phi _{{p},2D}.\end{equation}

Hereafter, whenever the experimental data are plotted, their volume fraction would be converted to their 2-D equivalent according to (3.9).

3.5. Scaling laws for shock attenuation

To assess the key role played by the χ in the upstream and downstream pressures of the particle bed impinged by the IS of varying intensity, we plot the variations in $\mathit{P}_{{s},{UI}}^{\mathrm{**}}$ and $\mathit{P}_{{s}, {DI}}^{{*}}$ with χ in figures 15(a) and 15(d), respectively, wherein the symbols are rendered with the flow Mach number behind the IS, M ₂. Here, $\mathit{P}_{{s},{UI}}^{\mathrm{**}}$ is scaled via

Figure 15.

Variations in the $\mathit{P}_{{s},{UI}}^{\mathrm{**}}$ with the increasing χ (a), θ (c) and χ ^* (e). Variations in the $\mathit{P}_{{s}, {DI}}^{{*}}$ with the increasing χ (b), θ (d) and χ ^* (f).

(3.10) \begin{equation} P_{{s},{UI}}^{**}=\frac{P_{{s},{UI}}-P_{2}}{P_{5}-P_{2}},\end{equation}

where P ₂ is the pressure behind the IS. Compared with the wide variations observed in the θ dependences of $\mathit{P}_{{s},{UI}}^{\mathrm{**}}$ (figure 15 b) and $\mathit{P}_{{s}, {DI}}^{{*}}$ (figure 15 e), the deviations of $\mathit{P}_{{s},{UI}}^{\mathrm{**}}$ (χ) and $\mathit{P}_{{s}, {DI}}^{{*}}$ (χ) from the fitted master curves are significantly reduced, with relative standard deviations of RRMSE _UI,χ = 11 %, RRMSE _DI,χ = 35 %.

Closer examination of the χ dependences of $\mathit{P}_{{s},{UI}}^{\mathrm{**}}$ in figure 15(a) reveals that $\mathit{P}_{{s},{UI}}^{\mathrm{**}}$ for systems with similar χ varies with M ₂. Specifically, a system with a smaller M ₂ yields a lower $\mathit{P}_{{s},{UI}}^{\mathrm{**}}$ than those with the same χ but a larger M ₂. In fact, χ only indicates the nature of the shock-induced unsteady flows inside the granular media in terms of its location in the mode spectrum ranging from transmission-dominated to filtration-dominated. The intensity of flows, which is closely related to the long-term friction interaction between gases and particles, also plays an important role (Mehta et al. Reference Mehta, Salari, Jackson and Balachandar2019b ) because the gradual build-up of the upstream and downstream pressures depends on the frictional process. Specifically, the stronger interstitial flows lead to a relatively enhanced friction process, which in turn strengthens $\mathit{P}_{{s},{UI}}^{\mathrm{**}}$ , equivalent to a system with a larger χ. Thus, a modified χ ^* is proposed via multiplication of χ with 1/M ₂, χ ^* = χ/M ₂, whereby the effects of the flow intensity are accounted for. Using χ ^*, all $\mathit{P}_{{s},{UI}}^{\mathrm{**}}$ satisfyingly collapse onto the master curve (see figure 15 c), which can be fitted by

(3.11) \begin{equation} P_{{s},{UI}}^{**}=\frac{1}{1+\left(\chi ^{*}\right)^{0.63}/5.2},\end{equation}

with a relatively low relative standard deviation of RRMSE _UI,χ* of 12 %.

The χ ^* can also improve the scaling of $\mathit{P}_{{s}, {DI}}^{{*}}$ , as shown in figure 15(f), with a reduced relative standard deviation of RRMSE _DI,χ* of 22 %. A consistent χ ^* dependence of $\mathit{P}_{{s}, {DI}}^{{*}}$ emerges from figure 15(f) and can be fitted by

(3.12) \begin{equation} P_{{s},\mathrm{DI}}^{*}=\frac{\chi ^{*}}{3.4+\chi ^{*}}.\end{equation}

The considerable improvement in the scaling of $\mathit{P}_{{s}, {DI}}^{{*}}$ via χ ^* can be understood in terms of the thermal effects of shock-induced unsteady flows on $\mathit{P}_{{s}, {DI}}^{{*}}$ . As discussed in Appendix C, the increase in $\mathit{P}_{{s}, {DI}}^{{*}}$ , particularly for the filtration-dominated mode (Mode III), is caused by the heated gas layers between CS₁ and CS₂. A stronger IS results in hotter gases, which in turn results in an elevated $\mathit{P}_{{s}, {DI}}^{{*}}$ , which is equivalent to a larger χ. The correlation between the intensity of the IS and $\mathit{P}_{{s}, {DI}}^{{*}}$ can be effectively accounted for by χ ^*. Notably, a major portion of the experimentally attained results fall just to the right of the master curve, which may be due to the overestimation of the corresponding $\mathit{L}_{{TS}}^{\mathrm{*}}$ and the resulting χ. The scaling law of $\mathit{L}_{{TS}}^{\mathrm{*}}$ expressed in (3.8) is based on 2-D configurations whereby the particles are spherical cylinders. It is known that wave scattering and viscous losses in 3-D particle beds consisting of spheres are more significant than those in 2-D beds. The $\mathit{L}_{{TS}}^{\mathrm{*}}$ for the 3-D bed should be smaller than its 2-D counterpart with the equivalent φ _p . The angularity and roughness of the particles further increase the degree of shock attenuation. Therefore, the $\mathit{L}_{{TS}}^{\mathrm{*}}$ for the experimentally investigated particle beds, which are mostly composed of crushed particles, may be much smaller than those predicted by (3.8). A similar phenomenon is observed for χ ^*. The actual locations of these data in the $\mathit{P}_{{s}, {DI}}^{{*}}$ - χ ^* parameter space should be somewhat moved left towards the lower limit of χ ^*. A better coincidence between these data and the master curve is expected.

Equations (3.7), (3.8), (3.11) and (3.12) constitute a complete procedure to theoretically predict $\mathit{P}_{{s},{UI}}^{\mathrm{**}}$ and $\mathit{P}_{{s}, {DI}}^{{*}}$ . For a given combination of a particle bed defined by φ _p , d _p and L ₀ and an IS with a M _s , we first calculate the dimensionless parameter η using (3.7) and then substitute it into (3.8) to obtain the maximum transmission $\mathit{L}_{{TS}}^{\mathrm{*}}$ , which in turn determines the χ ^*, χ ^* = $L_{TS}^{*}/(L_{0}^{*}\cdot M_{2})$ . Equations (3.11) and (3.12) provide the estimations of $\mathit{P}_{{s},{UI}}^{\mathrm{**}}$ and $\mathit{P}_{{s}, {DI}}^{{*}}$ , respectively, on the basis of χ ^*. Since χ ^* incorporates the effects from all potential structural parameters that contribute to the shock-induced unsteady flows and the resulting $\mathit{P}_{{s},{UI}}^{\mathrm{**}}$ and $\mathit{P}_{{s}, {DI}}^{{*}}$ , it allows us to readily discern the dependence of $\mathit{P}_{{s},{UI}}^{\mathrm{**}}$ and $\mathit{P}_{{s}, {DI}}^{{*}}$ on these parameters. Figure 16 presents the variations of $\mathit{P}_{{s},{UI}}^{\mathrm{**}}$ (figures 16 a and 16 b) and $\mathit{P}_{{s}, {DI}}^{{*}}$ (figures 16 c and 16 d) in the parameter space (φ _p , d _p ) for short ( $\mathit{L}_{{TS}}^{\mathrm{*}}$ = 100, figures 16 a and 16 c) and long ( $\mathit{L}_{{TS}}^{\mathrm{*}}$ = 5000, figures 16 b and 16 d) columns subjected to an IS with M _s = 1.67.

Figure 16.

Variation diagrams of the predicted $\mathit{P}_{{s},{UI}}^{\mathrm{**}}$ (a) and (b), and $\mathit{P}_{{s}, {DI}}^{{*}}$ (c) and (d) in the parameter space (φ _p , d _p ) for short ( $\mathit{L}_{0}^{\mathrm{*}}$ = 100, (a) and (c)) and long ( $\mathit{L}_{0}^{\mathrm{*}}$ = 5000, (b) and (d)) columns subjected to an IS with a M _s = 1.67.

3.6. Comparison of different pressure correlations

For the theoretical prediction of P _s,UI , Daniel and Wagner (Daniel & Wagner Reference Daniel and Wagner2022) proposed a pressure correlation given in (3.13) by modifying the solid wall reflected shock Mach number, $M_{r_{wall}}$ , using an empirically determined function of the volume fraction, $\mathit{M}_{r_{wall}}^{2}\phi _{{p}}^{0.1}$ . Here, the power of 0.1 was determined using a least-squares fit of the experimentally measured $\mathit{P}_{{s},{UI}}^{\mathrm{**}}$ :

(3.13) \begin{equation} \frac{P_{s,UI}}{P_{2}}=\frac{2\gamma M_{r_{wall}}^{2}\phi _{p}^{0.1}-\gamma +1}{\gamma +1}.\end{equation}

The scaled $\mathit{P}_{{s},{UI}}^{\mathrm{**}}$ is hence expressed as follows:

(3.14) \begin{equation} P_{s,UI}^{**}=\frac{P_{s,UI}-P_{2}}{P_{5}-P_{2}}=\frac{M_{r_{wall}}^{2}\phi _{p}^{0.1}-1}{M_{r_{wall}}^{2}-1}.\end{equation}

Clearly, the $\mathit{P}_{{s},{UI}}^{\mathrm{**}}$ given in (3.14) takes into account the decreased strength of the RS caused by different φ _p values but fails to incorporate the d _p and $\mathit{L}_{0}^{*}$ dependences of the RS strength, which significantly decreases as φ _p approaches the lower limit of the dense granular regime, as shown in figures 16(a) and 16(b). Plotting the empirical estimations of $\mathit{P}_{{s},{UI}}^{\mathrm{**}}$ (denoted as $\mathit{P}_{\textit{UI},\textit{emp}}^{\mathrm{**}}$ in figure 17 a) against the simulation-derived $\mathit{P}_{{s},{UI}}^{\mathrm{**}}$ (see figure 17 a) reveals the significant limitation of (3.13). In particular, (3.14) yields a negative value for cases with quite low φ _p .

Figure 17.

(a) Comparison between the empirical estimations of $\mathit{P}_{{s},{UI}}^{\mathrm{**}}$ (calculated from (3.13) and denoted as $\mathit{P}_{\textit{UI},\textit{emp}}^{\mathrm{**}}$ in (a)) and the actual $\mathit{P}_{{s},{UI}}^{\mathrm{**}}$ . (b) Comparison between the empirical estimations of $\mathit{P}_{{s}, {DI}}^{{*}}$ (calculated from (3.16) and denoted as $\mathit{P}_{\textit{DI},\textit{emp}}^{\mathrm{*}}$ in (b)) and the actual $\mathit{P}_{{s}, {DI}}^{{*}}$ .

Daniel and Wagner (Daniel & Wagner Reference Daniel and Wagner2022) also established a semiempirical pressure correlation for P _s,DI through an observed relationship between the pressure drop of the particle curtain and a term that approximates the induced dynamic pressure,

(3.15) \begin{equation} P_{s,UI}-P_{s,DI}\propto u_{2}^{2}\phi _{p}^{0.5}\rho _{1}.\end{equation}

The modified dynamic pressure term, $\mathit{u}_{2}^{2}\phi _{{p}}^{0.5}\rho _{1}$ , neglects the relative velocity between the particles and the induced flow and replaces the density of interstitial gases with the initial density ρ ₁. The plot of the pressure drop, P _{s,UI ^–} P _s,DI, versus the modified dynamic pressure reveals a constant linear relationship given by the best-fit line as

(3.16) \begin{equation} P_{s,UI}-P_{s,DI}=9.6u_{2}^{2}\phi _{p}^{0.5}\rho _{1}.\end{equation}

The scaled downstream pressure $\mathit{P}_{{s}, {DI}}^{{*}}$ can be derived from (3.16) as follows:

(3.17) \begin{equation} P_{s,DI}^{*}=P_{s,UI}^{**}\frac{P_{5}-P_{2}}{P_{2}-P_{1}}-\frac{9.6u_{2}^{2}\phi _{p}^{0.5}\rho _{1}}{P_{2}-P_{1}}+1.\end{equation}

Similar to the empirical estimation of $\mathit{P}_{{s},{UI}}^{\mathrm{**}}$ (3.14), the $\mathit{P}_{{s}, {DI}}^{{*}}$ expressed in (3.17) does not depend on d _p and $\mathit{L}_{0}^{*}$ which becomes non-trivial as the φ _p approaches the lower limit of the dense granular regime, as shown in figures 16(c) and 16(d). In contrast with the good agreement between our theoretical predictions of $\mathit{P}_{{s}, {DI}}^{{*}}$ and the actual $\mathit{P}_{{s}, {DI}}^{{*}}$ derived from the experiments and simulations (figure 15 f), the empirical estimations (denoted as $\mathit{P}_{\textit{DI},\textit{emp}}^{\mathrm{*}}$ in figure 17 b) exhibit large deviations from the actual $\mathit{P}_{{s}, {DI}}^{{*}}$ , as indicated in figure 17(b). Specifically, we observe negative values of $\mathit{P}_{\textit{DI},\textit{emp}}^{\mathrm{*}}$ for cases with low φ _p, which is unphysical and again demonstrates the limitations of (3.17).