Influence of particle non-dilute effects on its dispersion in particle-laden blast wave systems

Multiphase systems of solid granular particles interacting with gaseous blast waves generated by the sudden release of large amounts of energy are widely found in nature and industrial applications. Such systems are responsible for the creation of complex turbulent structures of cosmic scales in supernova remnants.^1,2 Volcanic eruptions with volcanic ash as the dispersed phase created from the explosion itself involves the interaction of these granular particulates with the multiple blast waves generated by the de-pressurization of the high temperature magma.^3–6 In industrial processes where reactive and highly combustible powders generate complex multiphase systems, the prediction and control of dust explosions is of great concern from a safety standpoint.^7–9 In defense applications, metallic explosives are loaded with solid particulates to control as well as increase blast pressure.^10–12 To control and mitigate these high explosion pressures, usually non-reactive particles are used, while to enhance the same, reactive particles are preferred. At laboratory scales, instabilities in shock-driven multiphase flows play a key role in laser-driven inertial confinement fusion experiments.¹³ Because of these many applications, understanding the fundamental physics involved in these multiphase systems and being able to predict the motion of gas and particles is critical.

Although the above systems are complex, one can construct a representative but simplified multiphase canonical system based on particle-laden blast waves (PLBWs). Figure 1 shows the gas features such as the Main Shock (MS), the Second Shock (SS) discontinuities formed by the rapid expansion of a spherical high energy density charge of air as well as the Particle Front (PF) produced by the dispersion of the embedded particles at later stages of system evolution. The gas phase portion of the PLBW problem has been studied extensively by researchers using numerical as well as experimental methods. The problem also called the “spherical shock-tube blast”¹⁴ or the “bursting-sphere blast”¹⁵ was studied experimentally by Boyer et al.¹⁶ Numerical studies have been carried out by Brode et al.^17,18 and later by Liu et al.¹⁹ to capture and study the evolution of gas phase discontinuities in the pressure field. The spherical shock-tube blast problem has also been investigated analytically by Friedman²⁰ and McFadden²¹ to derive simplified expressions for the early stage evolution of MS and SS in spherical blasts.

Introduction of solid particles of finite inertia makes the analysis of the system difficult due to the additional levels of unsteadiness and non-linearity introduced by the coupling between the gas and particle phases. On the numerical front, these multiphase systems were investigated by researchers such as Theofanous and Chang²² using the Eulerian–Eulerian method by assuming both the gas and particle phases represented by continuum fields. However, the method that is inexpensive and tractable is applicable only at very high volume fractions or in a dense particle phase and inherently assumes mechanical equilibrium between the gas and particle phases. To tackle low to moderate particle volume fractions, researchers use the Eulerian–Lagrangian method, which considers the particles to be discrete point representations while assuming the underlying gas phase to be in the continuum regime. The method also provides the flexibility to use polydisperse particle distributions to emulate more realistic systems.^23,24 The EL method has gained renewed interest in recent times in modeling the interaction of inertial particles with gas shocks. The studies in this regard by researchers such as Ling et al.²⁵ and Shallcross et al.²⁶ studied the dispersion of particles from a dust curtain due to its interaction with a shock wave. There have also been studies by Kosinski et al.,²⁷ Gosteev et al.,²⁸ and Marayikkottu et al.²⁹ in understanding the dust entrainment of particles behind a shock wave from a substrate using the EL method. Researchers such as Balachandar et al.^30–32 have used the Eularian–Lagrangian method to study the evolution of canonical PLBW systems, but all previous studies of PLBW systems have focused on particle volume fractions in the dilute regime where the effect of inter-particle interactions is insignificant to the overall dynamics of the multiphase system evolution. The objective of this work is to study the effect of non-dilute aspects of the particle phase on the unsteady PLBW system such as inter-particle collisions and near-particle aerodynamic interactions, as shown in Fig. 1) using a Eulerian–Lagrangian numerical approach.

With respect to the first effect, Brilliantov et al.,³³ Vinkovic et al.,³⁴ and Chinnappan et al.³⁵ have demonstrated that inter-particle collisions may be modeled using a modified version of the Direct Simulation Monte Carlo (DSMC) method. In addition, soft-sphere-based deterministic model by Cundall and Stack³⁶ is used by Shallcross et al.,²⁶ while Harris and Crighton’s stress-based model³⁷ is used by Ling et al.²⁵ to model inter-particle collisions. The DSMC method, however, is less costly computationally through its use of coarse grained simulated particles instead of physical particles, thus making it more tractable for simulations with a large number of particulates. However, the second factor related to aerodynamic interactions between neighboring particles due to high particle volume fraction and their effect on drag is less well known. Analytical models are difficult to develop because they must account for the surface of every particle, creating an N-body problem, and experiments are hindered by the difficulties of measuring the force and the local flow field on an individual particle in a particle cloud. In the last 20 years, the Immersed Boundary Method (IBM) based Navier–Stokes solvers have been used by researchers such as Hill et al.,³⁸ Van Der Hoef et al.,³⁹ Beetstra et al.,⁴⁰ and Tenneti et al.⁴¹ to deduce quasi-steady drag laws for dense monodisperse particle distributions. In recent times, researchers such as Mehrabadi et al.,⁴² Hosseinzadeh-Nik et al.,⁴³ Osnes et al.,⁴⁴ and Mehta et al.⁴⁵ used particle-resolved simulations to understand the complicated, unsteady effects associated with the interaction of a shock with a dense distribution of particles. In this study, we will employ our DSMC code CHAOS,⁴⁶ which is based on highly efficient unstructured/octree automatic mesh refinement techniques⁴⁶ to model the flow inter-leaving through complex structures of monodisperse dense distributions of spherical particles to obtain drag laws that account for aerodynamic interactions.

The remainder of the paper is organized as follows. Section II discusses the Eulerian–Lagrangian numerical methodology used in this work to model the unsteady two-phase PLBW flow as well as validate our use of the inelastic-DSMC algorithm to model energy transfer between particles. In Sec. III, we obtain parameters for a particle cloud drag model using the traditional DSMC gas-flow algorithm for particle distributions of 5% and 10% volume fraction. In the  Appendix, we provide additional information regarding the statistical variation of the aerodynamic interaction drag force due to local particle volume fraction such as clustering. Finally, in Sec. IV, we use the aforementioned models to the PLBW for cases spanning a range of volume fractions and charge energy. The main conclusions along with recommendations for future work are given in Sec. V.

The Eulerian–Lagrangian solver is implemented in the original FLASH⁴⁷ research code framework developed at the University of Chicago. The Eulerian gas field is solved using the finite volume (FVM) solver provided with the original code, while an in-house developed Lagrangian particle solver is integrated into the FVM solver to numerically study the evolution of the multiphase PLBW system under consideration.

The continuum gas phase is governed by the (generalized) compressible Navier–Stokes equation given as

where $U(x→,t)$ is the conservative state vector as a function of space $x→$ and time $t$⁠, $Fd(U)$ is the inviscid flux, and $Fv(U∇U)$ is the viscous flux. The independent state and flux vectors are expanded as follows:

where $ρg$ is the fluid density, $u→g$ is the fluid velocity vector, $Pg$ is the gas pressure, $Eg$ is the total energy of the gas per unit mass, and $I$ is the identity matrix. The fluid viscosity stress tensor and the viscous heating terms are represented by $τ→g$ and $qg$⁠, respectively, and an ideal gas equation of state is used to enforce closure of the set of governing equations given above.

The dispersed or the particle phase is modeled as discrete point representations. Since the particle distribution of the system under consideration quickly expands causing the volume fraction to rapidly decrease, the effect of particle field on the gas is not considered in this study.⁴⁸ The evolution of these Lagrangian representative particles is governed by the numerical integration of Newton’s equations of motion given as

where $x→p$ represents the position vector of a particle in the domain, $u→p$ is the velocity of the particle, the mass of the particle is indicated by $mp$⁠, $F→p$ represents the force established on the particle by momentum transfer with the underlying gas field, and $F→coll$ represents the forces experienced by particles due to inter-particle collisions. Note that the force, discussed below, is zero for the dilute particle cases. Unlike the work of Balachandar et al.,³² compressibility-rarefaction and unsteadiness effects were not considered in Eq. (3) due to the large size and high inertia of particles in this study. For example, Ling et al.³⁰ suggests that for the density ratio between the particulates and the gas on the order of $O(103)$ or higher, the effect of unsteady forces on the particle dynamics is not significant. The second-order accurate Leapfrog method⁴⁹ is used for this time integration. The fluid–particle interaction force $F→p$ for the study is based on the fluid drag given by

where $CD$ is the drag coefficient formulated as a function of non-dimensional parameter and particle Reynolds number (⁠$Rep$⁠). The subscripts $g$ and $p$ on velocity $u→$⁠, and density $ρ$ indicate continuous (gas) phase and particle (disperse) phase, respectively. $Dp$ represents the particulate diameter. For this study, $CD$ formulation by Sternin et al.⁵⁰ is used as a function of particle Reynolds number $Rep$ given as

where $ρg$ and $μg$ are the gas density and viscosity in the vicinity of the particle. The characteristic length scale of the particle is the diameter of the spherical particle represented by $Dp$⁠. Note that the gas viscosity varies as a function of gas temperature according to Sutherland’s viscosity formulation⁵¹ in these simulations. A separate comparison of the currently used drag model with that of Loth et al.⁵² and a recent DNS study by Nagata et al.,⁵³ which does include compressibility and rarefaction effects, showed a maximum deviation of 3.9% and 3.5%, respectively, for the range of flow parameters of interest, i.e., Mach number $<0.8$ and Knudsen number $∼7×10−4$⁠.

Up to this point, the particle phase implementations pertains to the dilute particle regime (⁠$<5%$ volume fraction). When the local volume fraction of particles is higher than this volume fraction, dense particle effects such as inter-particle collisions and aerodynamic interactions between particles become important. Inter-particle collisions may be modeled through an inelastic version of the Direct Simulation Monte Carlo (DSMC) method as was shown by Brilliantov et al.,³³ Vinkovic et al.,³⁴ and Chinnappan et al.,³⁵ where particulates are assumed to be inelastic and binary in nature. As in the original version of DSMC, the number of instantaneous collisions is decoupled from particle motion and is computed for each computational cell of volume $V$ containing solid particulates for a time interval $Δt$⁠, where the latter approximates the inter-particle collision time step as given by

In this expression, $Dp$ is the particle radius, $N$ is the number of simulated particles in the system, the ratio of the number of real solid particulates to simulated or computational particles is represented by $pFNUM$⁠, the maximum relative velocity between particles in a cell is represented by $u→pr,max$⁠, and the volume of the computational cell with edge length on the order of the particle mean-free path length is represented by $V$⁠. Note that, through the use of $pNUM$ in the above equation, the actual number of inter-particle collisions is computed, even though, the number of simulated particles used in the simulation are much lower than the actual number of physical particles as shown in Table I.

The acceptance–rejection algorithm is used to select particle pairs for collisions, based on a probability of collision between two randomly selected pairs of granular particles $i$ and $j$ that is given by

where the terms $|u→pi−u→pj|$ and $|u→pr,max|$ are the magnitude of the relative velocity of the selected particles and the maximum relative velocity of the particles in a computational cell, respectively. According to the algorithm, a pairwise collision occurs if a randomly drawn number (⁠$∈[0,1)$⁠) is less than the computed probability $Pcoll$⁠. If a collision does not occur, the particle pair retains its original velocity.

For a collision pair that undergoes an inelastic collision,

where $k→$ is a unit vector in a random direction between the particles, $mpi$ and $mpj$ are the mass of the particles, $u→pi$ and $u→pj$ are the pre-collision velocities of particles $i$ and $j$ with their post-collision velocities identified by primes, and $u→ij$ represents the relative velocity between the particles. The visco-elastic nature of the particles affects the two-phase flow system evolution since it affects the rate at which kinetic energy is dissipated. In this work, since the particles considered are rigid, a constant coefficient of restitution $en$ is used. Our approach, however, can easily be extended to include the dependence of $en$ on particle relative velocity.

In order to verify our inelastic-DSMC approach, we performed a zero-dimensional (0D) simulation on a domain of size $1×1×1mm3$ with periodic boundaries and no gaseous species. A total of 2500 spherical particles of uniform diameter 50 $μ$m and unit mass were initialized randomly in the domain with velocities sampled from the interval [$−50$⁠,50] m/s. The temporal evolution of the granular temperature, $Tgr=mp(u¯p−upi)2¯$⁠, was compared with an analytic result by Haff et al.⁵⁴ given as

where $Tgr(0)=8.267J$ is the granular temperature at time $t=0$ and $n¯$ is the number density of particles in the system. The relaxation time $τ$ is defined as the time taken by the system of monodisperse particles of constant coefficient of restitution $en=0.8$ to reach a quarter of its initial granular temperature. The temporal evolution of the system granular temperature obtained from our inelastic-DSMC simulation (shown in Fig. 2) is in good agreement with Eq. (10). The system evolves by dissipating particle kinetic energy through inelastic binary collisions and with increasing time we can see that more particles approach zero velocities in the simulation indicating homogeneous cooling as shown in Fig. 3.

In dense particle regimes, as the particles get closer to each other, the gas field in the vicinity of these particle clusters is also modified. This effect causes the drag model for an individual particle given by Eq. (4) to become insufficient in modeling the particle force terms in such a modified flow field. In order to select a viable mathematical formulation for the effect of aerodynamic interaction in dense particle distributions from the plethora of available models, a separate meso-scale study was performed as discussed in the Sec. III.

Using our DSMC code CHAOS,⁴⁶ which has been optimized to study gaseous flows through porous media such as a monodisperse dense distribution of spherical particles, we start with a three-dimensional cubical simulation domain of edge length 30 $μ$m for the simulations performed in this section. Since we are studying near-particle aerodynamic interaction effects in atmospheric gas pressure, half-Maxwellian subsonic boundary conditions are used on all the boundaries to emulate the atmospheric back pressure. An inert argon gas of free-stream gas number density $2.45×1024/m3$ (corresponding to 1 atm gas pressure) and a temperature of 300 K is used as the working fluid. Gas flow velocities corresponding to free-stream Reynolds numbers 10, 50, 100, and 200 based on the individual particle diameter as the characteristic length are used in the study and directed in the positive Z-direction. The corresponding free-stream Mach numbers are in the range 0.2 to 4.3. The particle distribution is initialized at the center of the simulation domain with a “buffer” distance of 5 $μ$m (⁠$>10$ times gas mean-free path) from the simulation domain boundaries to avoid boundary condition effects as shown in Fig. 4. The ratio of DSMC gas particles to the real number of argon gas molecules in the simulation or a gas-FNUM of 10 000 was used for the simulations with a maximum octree refinement of five with a time step size of $5×10−9$⁠. The results were sampled for a total of 20 000 iterations after the flow reached steady state to reduce the stochastic noise associated with these simulations.

A monodispersive distribution of particles is used in this study. Spherical particles of uniform diameter are randomly placed inside the cubical area identified by PP$′$QQ$′$ as shown in Fig. 4 using the Monte Carlo algorithm. If particle overlap is detected, the Monte Carlo approach used to generate the particle location relocates the particle such that the minimum inter-particle surface distance is 0.1 $μ$m. Each of the 2 $μ$m diameter sized spherical particles used in the study is modeled using 192 triangular surface panels and initialized with a surface panel temperature of 300 K. Two particle volume fractions at 5% and 10% shown in Fig. 5 are considered with a corresponding number of particles of 96 and 191, respectively.

Starting with the gas flow features in the vicinity of the particle distribution, gas Z-velocity contours for the two particle distributions of particle volume fractions 5% and 10% are shown in Fig. 6 for a free-stream particle Reynolds number of 100. The isometric view of the simulation domain shows three planes of the fluid domain cutting through the particle distributions with interior planes (plane 2) for the shown separately in Fig. 7. We can see clearly that the higher particle volume fraction (10%) induces more blockage to the flow as compared to the lower volume fraction. The set of particles in the distribution circumscribed by sonic isoclines are identified as aerodynamically interacting particles (AIPs), whereas isolated particles are identified as independent particles (IPs). The collection or cluster of aerodynamically interacting particles are identified in these figures using the contour lines corresponding to gas Mach number $<1$⁠. The study reveals that more AIPs are observed for higher particle volume fraction as compared to the lower volume fraction case. From Fig. 7, we can see that at 10% particle volume fraction, more volume of gas is decelerated to subsonic regime due to the proximal presence of a large number of spherical particles. Furthermore, the higher blockage means a higher drag force on the aggregate due to higher momentum transfer between the gas and the particle distribution will be generated. It is also interesting to see that the detached bow shock (DBS), similar to that observed in earlier modeling,⁵⁵ becomes a “collective” shock in the case of the higher volume fraction as opposed to separate shocklets formed around particles in the low volume fraction case. Indeed, this collective shock formation plays a major role in blocking the flow through the particle distribution and increases the momentum transfer rate and hence the drag force.

With the high fidelity gas flow field, we can calculate an average drag force that will include contributions due to clusters of particles. The net force on the particle distributions shown in Fig. 8 is computed from the momentum transfer between the gas particles and particle surfaces. The force in the Z-direction, $FZ$⁠, identified as the drag direction, was found to be at least two orders of magnitude higher than the lateral directions, X and Y, hence the lateral components were neglected. The average modified drag force experienced by an individual particle, $Fmod$⁠, in the distribution due to aerodynamic interaction is deduced by dividing the net drag force in Z-direction with the total number of particles in the particular particle distribution, $Npart$⁠, as

Although this approach does not directly use the different morphologies identified in Figs. 6 and 7, we can assess the variation in $FZ$ due to variation in the local volume fraction shown in the green and magenta ellipses of Fig. 5. The  Appendix discusses further that analysis and the determination of that variability shown in the “error bars” of Fig. 8.

The particle cloud drag models of Di Felice et al.⁵⁶ and Tenneti et al.⁴¹ were considered in this study. In the former semi-empirical model, Di Felice et al.⁵⁶ developed a normalized drag model through comparisons with experiments of the form

The drag model for particle clouds was developed by Tenneti et al.⁴¹ through IBM gas-continuum simulations to give a normalized drag force of

where the remaining terms in the model are given as

and $Rep=U∞Dp/νg$ is the particle Reynolds number based on single particle diameter free-stream gas velocity $U∞$ and viscosity $νg$⁠. $ϕp$ is the particle volume fraction given by

where $N$ is the number of computational monodisperse spherical particles embedded in a computational cell and the ratio of number of true solid particulates representing a single DSMC computational particle is indicated by $pFNUM$⁠. Note that as the volume fraction approaches $→$ 0, the equations approaches $→Fisol(Rep)$(isolated sphere). At higher volume fractions, the equation predicts higher particle drag compared to single particle drag (⁠$Fisol(Rep)$⁠). Note that F$isol(Rep)$ is the drag on an isolated sphere given in Eq. (4) and recast here as

where $CD$ is given by Eq. (4) for the Di Felice model while $CD$ formulation of Naumann and Schiller⁵⁷ is used for Tenneti’s curve. We compare our new particle cloud drag results from the DSMC simulations with the drag models by Di Felice et al.⁵⁶ and Tenneti et al.⁴¹ for the two particle volume fractions of 5% and 10% in Fig. 8. Note that we normalize our particle drag force with that given by Stokes’s drag as

We see that for the two volume fractions considered, our simulation results match well with the drag model by Di Felice et al.⁵⁶ The particle cloud drag model by Tenneti et al.⁴¹ tends to over-predict the particle drag compared to our simulation result and that of Di Felici’s model. The root mean square error between the Di Felice and Tenneti models with our results for the 5% particle volume fractions are 1.6% and 2.7%, respectively, and 0.6% and 3.8% for Di Felici and Tenneti models, respectively, for the 10% volume fraction. According to Di Felice et al.,⁵⁶ the formulation is generally valid for particle volume fractions less than 50% and Reynolds number in the range $10−2$ to $104$⁠.

The error bars for the DSMC results show the maximum and minimum drag experienced by a particle within the distribution and the symbol represents the averaged value. It is to be noted that the upper bound of the error bar represents clustered particles or particles in the vicinity of a high local volume fraction. The lower bound generally represents an independent or a non-aerodynamically interacting particle and for all the cases considered here, the value is close to $f=1$⁠. For the 5% volume fraction case, in general, the drag force varies between $∼+77%and−73%$⁠. While for the 10% case, the variation is $∼140%and−80%$ about the average value.

With the confidence achieved from the gas kinetic simulation results and the good agreement with the Di Felice model, we use the Di Felice model to accommodate aerodynamic interaction effects in our PLBW simulations, discussed next.

A three-dimensional cubical simulation domain of edge length 1 m as shown in Fig. 1 is used in this study. The explosive charge is initialized as an isothermal sphere of gas of radius $Rc$⁠, centered at (⁠$a/2$⁠, $a/2$⁠, $a/2$⁠) in the simulation domain and “outflow” boundary condition is applied on all the boundaries of the domain. A diatomic model of air with a gas constant of 287.08 J/kg K corresponding to a molar mass of 28.96 g/mol is initialized in the domain with an ambient pressure, $P0$⁠, and temperature, $T0$⁠, of 1 atm and 300 K, respectively. The gas simulation parameters corresponding to original numerical studies by Brode et al.¹⁸ and later by Balachandar et al.³² are given in Table I, which gives the initialized temperature of the gas in the charge, $Tc$⁠, for the four cases as defined in Fig. 1 where the pressure of the charge gas, $Pc$⁠, is assumed to be 121 atm in all cases.

As discussed in Sec. II, the evolving flow field is obtained by solving the Navier–Stokes equations. The inbuilt Roe–Reimann solver⁵⁸ and the shock detection algorithm are used to capture and preserve the shock discontinuities in the gas flow field. The block-structured domain discretization provided by the PARAMESH⁵⁹ external library is used at an AMR (Adaptive Mesh Refinement) level of five, which was found to be sufficient in resolving the gas and particle phase features of the multiphase system under consideration. The computational blocks are further divided uniformly into 8000 (⁠$20×20×20$⁠) cells. Under this configuration, the minimum edge size of the block and the cells are 0.0625 and $3.125×10−3$ m, respectively. In order to enforce numerical stability in the simulations, a CFL criteria of 0.8 is used to determine the adaptive time step. For the conditions specified in this work, the maximum gas temperatures achieved are about 3000 K. Therefore, in contrast to real explosions that can experience temperatures as high as $∼10000$ K, high temperature effects such as gas radiation losses, dissociation, or ionization are not considered in this work.

Monodispersive distributions of spherical, irrotational, chemically inert point-force representations of particles of 100 $μ$m diameter are used in this study. The solid particles are distributed randomly inside the spherical region identifying the “explosion charge,” i.e., inside the spherical gaseous region of gas of radius $Rc$ in the Eulerian simulation domain as shown in Fig. 1. A constant material density of $1410kg/m3$ is used for all the particles in all the simulations, making the mass of each particle equal to $∼7.3×10−10$ kg. The temperature of the spherical particles is initialized to 300 K and since the maximum initial gas temperature considered in cases 2, 3, and 4 is 300 K, the particle to fluid temperature ratio $TR∼1$⁠. Therefore, the effect of particle heating or cooling on the dynamics of the particles can be neglected according to Nagata et al.⁶⁰ 

Table I summarizes the physical and numerical parameters for the four simulation cases considered in this work. Note that only cases 2, 3, and 4 have particles embedded in the flow field. For case 2, 500 representative spherical particles are used to emulate a particle volume fraction of 0.03%. Since this volume fraction corresponds to a dilute particle regime, dense particle effects such as inter-particle collisions and near-particle aerodynamic interactions will not be important, but the case provides a baseline to understand these dense particle effects. As shown in the table, the integer variable $pFNUM$ is used to represent the number of physical particles represented by a single computational particle. In case 3 where only inter-particle collisions are modeled, two volume fractions corresponding to dilute vs non-dilute are considered. Finally, in case 4, the same volume fraction as the non-dilute value for case 3 is used, but now both inter-particle collisions as well as near-particle aerodynamic interaction effects are modeled. For dense particle distributions in cases 3 and 4, the coefficient of restitution of the granular spherical particles is 0.6, a value chosen to approximately match the inelastic behavior of spherical glass beads of materials such as soda-lime or borosilicate.⁶¹ 

A clear understanding of the underlying gas field is important in predicting the evolution of the coupled particle phase. The unsteady, transient results obtained from case 1 simulation is compared to the scaling studies by Brode¹⁸ and Balachandar et al.³² for the same simulation parameters. Note that the latter mentioned studies were performed for a 1D spherically symmetric domain, whereas we performed a full three-dimensional simulation. To compare with their published $x$–$t$ diagrams, we follow their normalization approach. A non-dimensionalizing length scale (⁠$l∗$⁠) is obtained as a function of the initial total energy (⁠$Etotalg$⁠) and ambient pressure $P0$ as

then by assuming an ideal gas equation of state, the length scale can be further written as a function of the pressure ratio of the initialized charge pressure $Pc$ to the ambient gas pressure $P0$ and charge radius $Rc$ as

Non-dimensionalization of physical time is performed using the quantity $t∗$ which is given as

where $c0$ is the speed of sound in the ambient atmospheric conditions and the ratio of specific heats of the gas is represented by $γ$⁠.

The normalized $x$–$t$ diagram (phase diagram) for our case 1 simulation is in good agreement with the numerical results obtained by Balachandar et al.³² as shown in Fig. 9. We notice that the main shock (MS) moves radially outward and eventually reaches the acoustic limit defined by the local speed of sound. The second shock (SS) generated at the charge surface initially moves outward, but at later times move inward, back toward the charge center, and then reflects off the center as shown in the phase diagram (Fig. 9). Also marked in the figure are three specific points in time $A$⁠, $B$⁠, and $C$ corresponding to physical times 62.5, 187.5, and 250 $μ$s, respectively, where the additional macroparameter contours of the evolving flow will be presented.

For example, the gas pressure contours shown in Fig. 10 at these times identify the radially outward propagating Main Shock (MS) and the decreasing peak over pressure from $∼$8 atm at 62.5 $μ$s to $∼3$ atm at 250 $μ$s. The figures also show the SS moving outward at times $A$ and $B$⁠, reflecting at the explosion center at $∼250$ $μ$s $(C)$ and then moving outward at later times. We can clearly see that the SS has much lower peak pressure magnitudes compared to the MS. The temperature contours shown at different times after simulation initialization in Fig. 11 indicate that the high temperature gas is confined to the explosion center. The non-concentric gas features seen in these figures can be attributed to the instabilities developed in the vicinity of the explosion center due to gas density gradients. The gas temperature at the explosion core is increased by the focusing effect of the SS with time as seen in the sequence $B$ to $C$⁠. This effect was also observed and reported by Brode et al.¹⁸ using finite-difference numerical simulations. Since gas radiation losses are not included in this work, the main source of temperature reduction is due to gas conduction. This thermal energy is unavailable for increasing the shock strength and the high temperature core merely affects the speed with which the SS propagates. For the higher temperature case, the SS moves faster lowering the central intersection point represented by point $C$ in Fig. 9, causing the lower core temperature to decrease the shock speed and hence move the intersection point to a larger value on the time scale.

With the good agreement obtained in the gas phase modeling, we now introduce a dilute distribution of monodispersive spherical particles in the spherical gas charge as shown in Fig. 1. The gas phase and particle phase simulation parameters in the study corresponds to case 2 in Table I. Note that the initial gas temperature ratio $Tc/T0$ for this study is one as compared to ten for case 1. Figure 12 shows the $x$–$t$ phase diagram of various flow features developed in the multiphase system. The MS and SS in the gas phase follow the same temporal evolution characteristics as discussed for case 1. Since the initial temperature ratios, $Tc/T0$ is lower than case 1, the gas in the vicinity of the explosion center for case 2 is at a lower temperature than case 1 at all times, reducing the speed of SS. Therefore, the SS interacts with the center of explosion much later in case 2 at $τ∼1.5$ (labeled C in Fig. 9) as compared to $τ∼0.3$ in case 1. The particle front (PF), identified by the farthest particle from the core at every time instant, moves radially outward and can be observed to lag behind the MS due to the high inertia of the particles and is not significantly affected by the weak SS. Although the general characteristic of the PF line-plot is similar to the result by Balachandar et al.,³² a minor discrepancy can be observed between the two predicted PF’s. The discrepancy in the PF profiles can most likely be attributed to the difference in dimensionality of the two simulations. The original work by Balachandar et al.³² considered a one-dimensional simulation with spherical boundary conditions at the explosion center, while we have performed a complete three-dimensional computation. The three-dimensionality of the current work allows the particles to relax in all the three spatial directions, hence reducing the speed of the PF as compared to the previous one-dimensional study.

Figure 13 shows the instantaneous position of the dispersed particles on the quarter domain superimposed over gas pressure contour at different (normalized) times after explosion. Consistent with the phase diagram shown in Fig. 12, the PF shown in Fig. 13(a) lags behind the MS at all times. In Fig. 13(b), the SS can be seen to be ahead of the PF, while Fig. 13(c) indicates that the PF eventually overtakes the SS. Note that Fig. 13(c) corresponds to a time when the MS has already left the domain.

The distribution in terms of the fraction of particles at different radial locations, $r$⁠, with respect to the explosion center as the origin is shown in Fig. 14(a) for four normalized time instances $τ=0.12,0.88,1.65$⁠, and 1.91. The right end of each curve approximately captures the radial position of the PF with respect to the explosion center at that indicated time. We can clearly see that with progressing time, the particles are moving radially outward. The peak probability value of the distribution for all four time instances is relatively constant at $∼0.02$ and is located at the farthest radial location from the explosion center. This indicates that a major fraction of particles is located close to the PF at all times. A similar distribution of particle speeds is shown in Fig. 14(b) for the four time instances. The figure suggests that a majority of the particles accelerate between $τ=0.12$ and 0.88, reaching a maximum speed of $∼320$ m/s at $τ=0.88$⁠. At later times, $τ=1.65$ and 1.91, the drag force experienced by the particles decelerates the particle distribution to $∼240$ m/s. The peak of the velocity distributions for all four time instances shown in the figure indicates that a large portion of particles in the PLBW system have high speed magnitudes and hence high kinetic energies. In dilute particle distributions such as case 2, the fluid induced drag is the sole player responsible for the decrease in particle kinetic energy at later stages of system evolution.

Turning now to the dense particle cases, 3 and 4, the general characteristics of the evolving particle front is the same as discussed for cases 1 and 2. First, the effect of inter-particle collisions is studied in case 3 in the absence of near-particle aerodynamic effects. Figure 15 shows the evolution of the different kinetic energy components of the multiphase PLBW system for case 3. The particle kinetic energies are summed over the discrete particles in the system at any given time step. Only $∼0.5%$ of the deposited energy is converted to gas kinetic energy with a peak KE of $∼1,200$ J as the main shock is generated at $τ∼0.25$⁠. The gas KE continuously dissipates with decreasing gas velocity as the shock propagates and expands outward. A second smaller peak in gas kinetic energy is observed at $τ∼1.2$ corresponding to the generation of the outward expanding SS. The gas exchanges momentum with the particle phase to increase the kinetic energy in the particle phase as the embedded particles start moving. Since the particles considered in the study do not have any internal energy modes, such as rotation or vibration, the entire momentum transferred between the expanding gas and the particle phase is converted into kinetic energy in the latter. The temporal evolution of the total kinetic energy of the particle phase indicated in the legend as C for initial volume fractions of 7.5% and 22.5% is shown alongside simulations repeated without including particle inelastic collisions, denoted “collision-less” or CL. The peak of the particle phase kinetic energy lags that of the gas phase kinetic energy due to the high inertia of the particles compared to the gas phase. As we can clearly see, the difference in the particle phase kinetic energies for the 7.5% volume fraction case between the collisional (C) and collision-less (CL) simulations is negligible with only a peak difference of $∼1.8%$ between the curves. On the other hand, for the non-dilute case with a 22.5% initial volume fraction, the difference between the particle phase kinetic energy evolution curves is evident. The collision-less (CL) particle phase achieves a higher peak kinetic energy magnitude as compared to collisional (C) case.

The temporal evolution of the particle averaged kinetic energy dissipated at each time step, shown in Fig. 16, indicates that the inelastic particle collisions is an energy removal mechanisms only in the early stages of the system evolution (⁠$τ<.6$⁠). The fluctuating nature of the graph can be attributed to the randomness associated with the inelastic-DSMC collision mechanism. The rapid expansion of the particle phase along with the variation in the relative velocity between the particles generates spatially and temporally fluctuating collision frequencies according to Eq. (6), which, in turn, makes the dissipative particle kinetic energy fluctuate with time. During this transient period, the particle volume fraction, computed as the volume occupied by the particles in a spherical volume of gas circumscribing the distribution, also drops rapidly as indicated in Fig. 16. This early stage decrease in the particle kinetic energy reduces the average relative velocity between gas and particles in the later stages of the two phase blast wave flow compared to a collision-less system, as shown in Fig. 17. Therefore, at later times, the system KE remains higher for the collisional C simulation, compared to the collision-less CL case because the dissipation of system kinetic energy due to particle drag is lower in the collisional case as can be seen in Fig. 15.

Finally, case 4 considers the evolution of a dense particle phase of initial volume fraction 22.5% with both inter-particle collision and near-particle aerodynamic interaction effects modeled using the inelastic-DSMC approach and the drag correction proposed by Di Felice,⁵⁶ respectively. Figure 18(a) compares the distribution of the radial position of particles in the system at four different normalized time instances, $τ=0.12,0.88,1.65$⁠, and 1.91 with the case 2 results. We can clearly see that the PF for the dense particle system moves much faster than the dilute particle system, case 2, at all four time instances shown in the figure. Similarly, the distribution of the particle speed at three instances of normalized simulation time, $τ=0.12,0.71$⁠, and 1.91 is compared to its counterpart for the dilute particle system shown in Fig. 18(b). The dilute particle system results shown in the figure suggests that a major fraction of the particles in the distribution accelerates from $∼240$ m/s at $τ=0.12$ to $∼340$ m/s at $τ=0.71$⁠. At later time, the dilute particle system decelerates to $∼240$ m/s by $τ=1.91$⁠. In contrast, a major fraction of particles in the dense system (case 4) assumes higher velocities of $∼400$ m/s at early times of system evolution (⁠$τ=0.12$⁠). The particles in these systems continuously decelerate from $τ=0.12$ to 1.91 approaching peak particle speed magnitudes of $∼180$ m/s by $τ=1.91$⁠. The implication of near-particle aerodynamic interaction on the PLBW system is that the particles achieve high speeds at early stages of system evolution, but at later stages, the particle phase decelerates drastically to much lower speed magnitudes. This also implies that the total kinetic energy of the particle phase in dense systems is higher than its dilute counterpart at early stages of system evolution, but at later stages, it decreases to much lower values than the collisional dense case (case 3) or dilute particle phase (case 2) due to the higher drag resistance experienced by the particle phase in its early stage of evolution.

An Eulerian–Lagrangian solver was developed using the FLASH research code framework to study the contribution of non-dilute particle effects on particle-laden blast wave systems. The evolution of the gas-particle coupled system was analyzed for both dilute and dense particle volume fractions to reveal key aspects of these complex PLBW systems. Non-dilute particle effects such as inelastic inter-particle collisions and near-particle aerodynamic interactions were considered in this study by including inelastic collisions and a validated aerodynamic interaction model. The inter-particle collision was modeled using a stochastic inelastic-DSMC based approach, while the near particle aerodynamic interactions in these systems were simulated using the Di Felice drag correction approach with parameters obtained from separate DSMC simulations that consider non-local variations in particle volume fraction.

Four different cases of monodisperse solid particles embedded in a spherical gaseous explosive charge with varying particle volume fraction, temperature ratio of initialized gas in the charge region to ambient, and degree of inclusion of non-dilute particle interaction effects were studied. An analysis of the dilute particle system (case 2) revealed that the fluid induced drag force is solely responsible for the early stage acceleration as well as the later stage deceleration of the particle distribution. A case with inelastic collisions without near-particle aerodynamic interaction effects particle dense system (case 3) revealed that the system retains more kinetic energy as compared to the particle-dilute case. With the inclusion of aerodynamic interactions to the collisional system (case 4), we see two competing effects in the early stage of the PLBW system evolution when the local particle concentration is high. Inter-particle inelastic collisions tend to reduce the kinetic energy of the particle phase through irreversible inter-particle collisions during this stage, while the close proximity of particles in the dense particle regions enhances the drag force experienced by the particles which tends to increase the kinetic energy of the dispersed phase. Note that this effect is prevalent only in the early stages of system evolution $τ<0.56$ where the particle volume fraction is $>5%$⁠. In the later stages of system evolution, the system with aerodynamically interacting particles showed high levels of deceleration through drag resistance as compared to the cases with only inelastic collisions or the dilute flow, hence expending more particle kinetic energy. These conflicting trends can be attributed to the fact that the reduced particle velocity in the early stages of case 3 reduces the magnitude of particle drag deceleration at later stages. While in the case of active near-particle aerodynamic interaction systems, the early stage increase in particle velocity amplifies the drag resistance of the particle distribution in the later stages of particle evolution.

Our study shows that the inclusion of dense particle effects such as inelastic collisions and aerodynamic interactions are essential in correctly predicting the evolution of particle-laden blast wave systems. From a safety point of view, these results show that even dilute particle systems have the potential to cause higher damage at larger distances compared to dense particle systems. This is due to the fact that particle dense systems retain higher particle kinetic energy at smaller distances and readily dissipate kinetic energy at larger distances from the explosion center, reducing the possibility of damage at larger distances compared to dilute systems. Future work will target the evolution of polydisperse particle distributions that would describe a more realistic system. Since these systems will consider particles of diameter of $∼O(1)$ micrometers, the effect of unsteadiness and compressibility of the particle force models will also be taken into account.

All authors contributed equally to this work.

This work is sponsored by the Department of the Defense, Defense Threat Reduction Agency under Award No. HDTRA1-20-2-0001. The content of the information does not necessarily reflect the position or the policy of the federal government, and no official endorsement should be inferred.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

As mentioned in Sec. III, although we perform an averaging over the 5% and 10% volume fractions shown in Fig. 5, it can be seen that there are spatially localized clusters of particles where the local volume fraction is greater than the nominal (Marked: A) as well regions where the local volume fraction is lower than the nominal (Marked: B). Due to this internal variation in local volume fraction, the force experienced by individual particles within the distribution is also expected to vary. To understand the effect of that variation, we first determine the local particle volume fraction and its variation in the domain and then assess the statistical variation in drag force.

The local volume particle volume fraction in the vicinity of each particle is approximately estimated by sub-sampling the domain. The number of spherical particles contained in a spherical volume of radius $Rsam$ with a particle center of mass as the origin is determined and the local particle volume fraction (⁠$VolFraci$⁠) corresponding to that particle, $i$ is determined as follows:

where $Nsam$is the number of particles in the vicinity of particle $i$ determined during sub-sampling and $dp$ is the individual particle diameter. Note that this procedure yields a very approximate distribution of volume fraction in space. The distribution of local volume fraction of the 5% and 10% distributions shown in Fig. 19 indicate that the spread of local volume fraction in the 10% nominal distribution is greater than the 5% distribution.

Next, the distribution of fluid induced forces (drag) within the particle cloud or particle distribution is considered for incoming flows with Reynolds number 100. The flow is considered to be directed in the positive Z direction (see Fig. 5). The force experienced by each particle in the particle cloud is computed by numerically integrating the panel-wise forces corresponding to each particle. First, the forces in the lateral direction (X and Y) are considered in Fig. 20. For both the distributions, the lateral forces are equitably distributed about zero on the X-axis. Hence, as mentioned in Sec. III, the net effect of lateral forces on the distribution is negligible.

Figure 21 shows the distribution of drag force (Z-direction) for the two volume fractions. The force distributions are normalized with the maximum value of the distribution, i.e., 2.51 and 2.67 $μ$N for 5 and 10%, respectively. In the case of 5% volume fraction, the spread in local volume fraction is small with a maximum of 6.28% and a minimum of 1.94%. Due to this narrow spread in local volume fraction, the individual drag force distribution is also narrow and more uniform as compared to the 10% distribution. In the case of the 10% distribution, the local volume fractions vary from 21.2% to 1.3%, therefore, the corresponding drag force distribution is also spread across a wider range. We also notice that the particles in the high volume fraction zone experience higher drag as compared to the low volume fraction zones. For the 5% case, particles in the neighborhood of local volume fraction of $∼6.28%$ experience drag force of 2.51 $μ$⁠, while zones at $∼1.94%$ volume fraction experience only 1.445 $μ$N of force. Similar is the case with the 10% nominal volume fraction case, particles at $∼21.2%$ volume fraction experience $4.67μ$N force, while only $∼0.219μ$N force is experienced by particles in a local volume fraction of 1.3%.