Numerical investigation of laser-driven shock interaction with a deformable particle

Shock interactions with non-uniform media are encountered in a variety of physical situations including astrophysical flows,^1–3 multiphase explosives,^4,5 shock propagation through bubbles,^6,7 and inertial confinement fusion.⁸ Shock–particle interactions are a general class of problems that involve shock propagation through media characterized by density or temperature inhomogeneities. The interaction is highly transient and nonlinear. A complex system of regular and irregular shock-wave reflection, diffraction, and shock focusing may exist, even in an idealized interaction between a shock and an isolated spherical non-deforming particle.⁹ 

The dynamics of a particle accelerating behind a shock wave is commonly described with its drag history. Several attempts have been made to develop analytical force models that predict the drag on a particle in an incompressible flow, some of which are reviewed by Michaelides et al.¹⁰ Parmar et al.^11,12 presented a force model for unsteady compressible flows and applied it to a shock traversing a spherical particle, demonstrating the importance of the inviscid unsteady contribution in capturing the peak force on particles behind the shock wave. Later, Ling et al.¹³ used the force model proposed by Parmar et al. and evaluated the unsteady forces over a range of shock-wave Mach numbers (M), particle Reynolds numbers (Re), and particle-to-host density ratios. The unsteady contribution had a short-time (i.e., order of acoustic timescale) influence on the particle evolution but captured the peak force during shock passage. Furthermore, its effect on the particle motion was inversely proportional to the particle-to-host density ratio.

However, these analytical models were formulated under the simplified assumptions of vanishing Re and M on a non-deformable spherical particle. Hence, the interactions of a deformable particle at finite Re and M must be studied experimentally or computationally due to the unavailability of theoretical solutions in the non-linear regime. In the weak shock regime, Haas and Sturtevant¹⁴ performed shock tube experiments to study the interaction of a shock wave (M < 1.3) with a single gaseous inhomogeneity. Multiple wave interactions including transmitted, reflected, and refracted waves were observed using optical diagnostics. Later, Igra and Takayama,¹⁵ Sun et al.,⁹ Martinez et al.,¹⁶ and Bordoloi et al.¹⁷ calculated time-dependent drag on a single particle under shock loading using experiments. Sun et al.⁹ presented time-resolved drag force measurements of a particle made of aluminum alloy subjected to a M = 1.22 shock. The peak unsteady force was unaffected by the viscosity of the surrounding flow and was an order of magnitude larger than the peak steady force. This was an important experimental observation which later reinforced Parmar et al.'s force model.¹¹ High-resolution numerical simulations^18–21 have explored pressure evolution around the non-deformable particles and computed the time-dependent drag coefficient during the passage of a shock over a spherical particle. Sridharan et al.¹⁹ conducted numerical simulations of shock propagation in air initially over a single aluminum particle and demonstrated that the computed drag coefficient decreases with increasing M.

Particle deformation is important in the context of multiphase explosives and high energy density (HED) systems where shock pressures typically vary from O(10) to O(100) GPa. Experimental data at these conditions are limited due to difficulty in conducting experiments. Hence, direct numerical simulations were utilized to observe the deformation of metal particles and measure the particle acceleration imparted by strong shocks.^4,5,22 Zhang et al.⁴ showed that metal particles, such as aluminum, beryllium, and magnesium, achieved 60% to 100% of the shocked explosive's velocity and were severely deformed. It was suggested that the particle drag model should account for the history-dependent particle shape as such deformation could modify the particle acceleration. More recently, numerical simulations of shock interaction with a deformable aluminum particle in nitromethane were performed for post-shock pressures up to 10 GPa.^23,24 Although particle deformation was modest during shock passage, it had a major influence on the particle drag behavior at later times, highlighting the need to include particle deformation in the force models.²⁴ At shock pressures higher than 500 GPa, Klein et al.¹ performed experiments on the Nova laser to understand the evolution of a high-density copper sphere embedded in a low-density plastic medium after the passage of a M = 10 shock as a model for shock–cloud interactions. 2D hydrodynamic simulations reproduced the experimental images. The sphere underwent considerable deformation at the initial stage and broke up at later times. Klein et al. revealed that at very high laser drives, the solid copper sphere was significantly preheated before the arrival of the shock, resulting in the interaction of a strong shock with a gaseous body rather than a cold solid. They also indirectly observed 3D vortex ring instabilities and demonstrated their role in late-stage cloud destruction, which was later confirmed by experiments on the OMEGA.^2,3

To accurately simulate material behavior at relatively low sub-eV temperatures, we require models with strength properties in the solid/liquid regime. Such models are not often included in most radiation-hydrodynamics simulation codes including FLASH,^25,26 which are used as tools to design HED experiments. In addition, high-temperature equation of state (EOS) models used for simulations become less predictive of the thermodynamic material properties necessary for describing the hydrodynamic processes taking place at low temperatures. For instance, FLASH only uses thermal pressure to compute the local sound speed of material. This neglects the existence of nonthermal pressure, which determines the behavior of shock-compressed solids. Theoretically, this leads to higher material compressibility even at low shock pressures.

In this paper, we develop a technique to implement a simple modification of ideal gas EOS for modeling solid materials in hydrocodes that lack material strength. Using this model, we study the problem of shock–particle interaction in solid targets relevant to laser-driven shock experiments under relatively less extreme conditions (⁠$≤100$ GPa). We aim to study the interaction of a laser-driven shock with an isolated particle embedded inside a low-density plastic target with an incident shock (IS) pressure high enough to exceed the yield strength but not necessarily melt or break apart the embedded particle during or after the shock passage, unlike shock–cloud experiments.^1,2

The paper is organized as follows. The target description, numerical methods, and EOS models used in this study are described in Sec. II. The results for a laser-driven shock interacting with an isolated particle are presented in Sec. III. We examine the pressure field, particle drag coefficient, and particle kinematics. We also discuss the deformation of the particle over time and quantitatively compare the evolution of particle diameter with that using the IONMIX²⁷ EOS tables. Finally, we study the dependence of particle size, particle density, and particle acoustic impedance on the velocity transmission factor. The concluding remarks are given in Sec. IV.

We consider a laser-driven shock propagation through a solid target consisting of a solid metal particle embedded inside a plastic host medium. The laser ablation-driven shock first traverses the host medium as shown in Fig. 1(a). The nature of shock refraction and shape of the transmitted shock (TS), as seen in Fig. 1(b), depends on parameters, such as the shock-impedance ratio and shock speed-ratio between materials across the interface.^28,29

The target configuration and temporal laser profile are shown in Figs. 2(a) and 2(b), respectively. The pre-shock material properties and shock Hugoniot parameters used for our simulations are listed in Table I. The target was driven by an incident laser beam of spot size of 800 μm diameter with spatio-temporally uniform intensity distribution. We employed IONMIX EOS tables to model the plastic ablator and the aluminum heatshield. We used the modified ideal gas EOS to describe plastic (CH), aluminum (Al), titanium (Ti), and tungsten (W), as will be described below. We comment here that widely used SESAME³⁰ EOS tables (e.g., SESAME 3720 for Al) were not compatible with the hydrodynamic code at pressure/temperature conditions considered in our work. The simulated SESAME data resided into the regimes of negative pressures, for which the code predicted nonphysical sound speeds as addressed by Farmakis et al.³¹ 

Plasmas contain electrons, ions, and thermal radiation due to the high-temperature field involved. The electron temperature is not necessarily equal to the ion temperature. Thus, the plasma is described by the “three temperature” (3T) approximation. The governing equations of the evolution of unmagnetized multi-temperature plasma are discussed in Tzeferacos et al.²⁶ and also described in the FLASH code user's guide.³² 

The multiphysics radiation-hydrodynamics FLASH code was used to carry out 2D Cartesian simulations to study the laser-driven shock interaction with a deformable particle. The governing equations are solved on an adaptive mesh refinement (AMR) grid using FLASH's unsplit scheme,³⁵ a finite-volume Godunov method consisting of a single-step, second order in time, directionally unsplit multidimensional data reconstruction-evolution algorithm, based on the corner transport upwind (CTU) method.³⁶ A third order in space reconstruction (piecewise parabolic method³⁷) is carried out using a minmod slope limiter along with a flattening technique to treat shocks.²⁵ The time advanced fluxes are computed using a HLLC Riemann solver³⁸ with second order accuracy. A Courant–Friedrichs–Lewy (CFL) number of 0.4 is used for numerical stability. The schematic of the computational domain is shown in Fig. 2(a). Outflow (zero-gradient) boundary conditions were imposed on all of the domain boundaries. A special treatment to the boundary condition on the particle interface is discussed in Sec. II C. We have neglected radiation in our study since it had an insignificant effect on the particle's overall hydrodynamic response.

As discussed in Sec. I, the absence of strength models along with the use of high-temperature EOS models in simulations overestimate the material compression and deformation at low temperatures and pressures. Hence, we sought to mitigate compression and deformation by employing a modified form of ideal gas EOS to model both the host and the particle,

where P, γ, ρ, and ε are total pressure, adiabatic index, mass density, and total internal energy, respectively. We have defined a constant average ionization inside the materials (e.g., a value of 1 is used for the Al particle). We believe that this is a good approximation for a particle that is heated by the TS. The temperature of the particle compressed by the TS is less than 1 eV in our study. The Saha ionization model shows that for Al at solid-like densities, average ionization is close to 1 at sub-eV temperatures.³⁹ Therefore, our choice of average ionization is justified by the theoretical model.

For a limiting case of strong shocks in an ideal gas, Rankine–Hugoniot relations imply that compression η and $us−ups$ ratio on the Hugoniot depend only on γ, where u_s and u_ps denote the shock velocity and post-shock particle velocity, respectively,

Therefore, we adjusted γ to match the experimental Hugoniot data.³⁴ In Figs. 3(a) and 3(b), we compare the Hugoniots generated with our modified EOS against those generated by IONMIX EOS. Our simulated Hugoniots agree well with the experiment, and in Fig. 3(b), the compression in Al shifted from the IONMIX-predicted value of ∼3.8 to ∼1.6 at 85 GPa. The stiffened gas EOS is another widely used model to describe liquids and solids under high pressures.^33,40,41 However, Fig. 3(a) shows that this model deviates from the experimental curve as the shock strength in the material increases. Conversely, the modified ideal gas EOS model is seen to perform better for a broader range of shock strength. In addition, Fig. 3(c) shows that the simulated $us−ups$ Hugoniots of CH, Ti, and W modeled using this technique agree well with their respective experimental Hugoniot curve.

Once γ is determined, the initial pressure is chosen using Eq. (4) to match the material sound speed, given as follows:

where P₀ and ρ₀ are initial the pressure and density of the material, respectively. We should note that, to appropriately model c₀, both the particle and the host could not be kept at equilibrium pressure at the initial state. However, to delay particle expansion until the arrival of shock, we numerically “freeze” (i.e., apply reflecting boundary condition at the particle interface) the grid cells that lie inside the particle. Once the shock meets the leading edge of the particle, we “unfreeze” those grid cells. This generates a pressure gradient across the interface and results in an increase in the particle width as discussed in Sec. III C.

The modified ideal gas EOS model is fully based on adjusting γ and c₀ to match the Hugoniot data from the LASL database.³⁴ However, due to the unavailability of experimental data in the database at much higher pressures than that presented in this work, we were not able to model our materials—Al, CH, Ti, and W—at such conditions. For example, the single-shock Hugoniot data for Al in the database spans up to ∼120 GPa. Figure 3(b) shows that our model provides excellent agreement with the experimental data up to ∼120 GPa. Ju et al.⁴² showed that the shock melting of Al occurs between 93 and 140 GPa. As such, we expect our model to be suitable for modeling materials at solid/liquid regimes (i.e., temperature and pressure of sub-eVs and few hundreds of GPa, respectively). Beyond these conditions, we should also account for other important physics in the EOS model such as temperature-dependent average ionization³⁹ and radiation effects, which are not accounted in the current model. We also expect the preheat effects to be negligible in our study. As shown in Nilsen et al.,⁴³ preheating is not an issue in low-drive experiments at pressures less than 130 Mbar. The shock pressure in our CH is ∼0.6 Mbar; hence, we are certainly in a low-drive regime where preheat should be negligible.

Time scales relevant to a shock–particle interaction problem are discussed by Mehta et al.⁴⁴ The time scale for a laser-driven planar shock passing through an isolated particle of initial diameter d_p is defined as follows:

The non-dimensional time t′ is defined as $t′=(t−ta)/τs$⁠, where t_a is the time the incident shock arrives at the leading edge of the particle.

Particle position (⁠$xp¯$⁠), particle velocity (⁠$up¯$⁠), and particle pressure (⁠$Pp¯$⁠) were numerically computed as mass-averaged quantities defined as follows:

where $ϕp¯$ refers to any field variable, such as particle pressure, and V_p is the volume of the particle.

Similarly, the inviscid force on a particle is defined as follows:

where P is the pressure acting on a surface of unit normal $n→$⁠, and A is the cross-sectional area of the particle.

The total inviscid force is computed in the streamwise direction as follows:

where k is the index of each cell-face that makes up the surface of the particle and $j→$ is the unit normal in the streamwise direction.

Finally, the unsteady drag coefficient for the particle is computed as follows:

where ρ_ps is the post-shock host density measured by probing a location that is $∼3dp$ distance away to the left of the particle.

A grid convergence test was performed for the case of an Al particle embedded in CH at a shock pressure of 55 GPa by increasing the refinement levels on the AMR grid. The number of grid points across d_p is denoted by N. The grid size was chosen by varying $dp/N$ and C_D was computed using Eqs. (7)–(9) at various grid sizes as shown in Fig. 4(a). The case with 70 points across the particle diameter captured the peak drag to less than 2% of the finest resolution tested. Hence, we performed our simulations at this grid refinement.

To verify our numerical schemes, we simulated a 6 GPa shock propagating in nitromethane over an embedded Al particle, as studied by Sridharan et al.²⁴ Figure 4(b) compares the time histories of C_D obtained from Sridharan et al.'s simulations and theoretical model with our numerical results. Although the theoretical model predicted a higher peak drag value than in simulations, both of the simulations showed good agreement in capturing the drag minimum due to shock focusing.⁹ The differences between the theoretical model and the simulations could be attributed to the zero deformation assumption made by the model.²⁴ 

We show pressure contours in Fig. 5 to highlight the flow dynamics as the shock propagates through the particle. The incident shock (IS) is planar before it reaches the particle, as seen in Fig. 5(a). As the shock meets the front of the particle, it experiences an impedance mismatch at the interface. Hence, a reflected shock (RS) travels back into the compressed CH and a transmitted shock (TS) travels downstream. In Fig. 5(b), the TS is planar due to the TS traveling at equal speed to the unrefracted shock outside. As the TS reaches the downstream end of the particle, it experiences an impedance mismatch at the interface. This generates a TS into the downstream CH and a reflected expansion wave back into the particle. In the CH, the diffracted shocks around the particle meet toward the center of the downstream end of the particle and further strengthen the shock. This phenomenon of shock focusing has been discussed in past works.^9,45 As this strong shock reflects back into the particle, it competes with the reflected expansion to produce temporal variations in pressure inside the particle [Figs. 5(c) and 5(d)]. Such pressure variations are plotted as mass-averaged particle pressure in Fig. 6(a) from t′ = 1 ns to t′ = 2.8 ns. When the laser is turned off at $t′=1.9$⁠, rarefaction waves start to catch up to the flow along the downstream direction. Once these waves meet the particle, the particle pressure monotonically decreases after t′ = 2.8 due to the influence of rarefaction stretching⁴⁶ on the compressed particle. Hence, the wave interactions are highly unsteady, thereby not allowing the particle pressure to equilibrate within a few τ_s.

The time history of the streamwise force on the Al particle is presented in Fig. 6(b). C_D quickly rises to its maximum value of 4.1 at $t′=0.55$ as the shock passes over the particle. As time progresses, C_D decreases and becomes negative due to shock focusing, indicating that the pressure downstream of the particle is larger than the pressure in front of it. The drag reaches a minimum value of $CD=−0.9$ at $t′=1.8$⁠. After that, the pressure inside and around the particle tends to equilibrate making the drag coefficient positive again. Once the laser is turned off at $t′=1.9$⁠, the compressed CH upstream of the particle starts to decompress due to rarefaction stretching. This decelerates the particle and results in another drag minimum at $t′=3.2$⁠.

These results confirm the importance of unsteady forces to the bulk motion of the particle. If we compare C_D from our simulations to that from previous works,^19,22,24 we find qualitative agreement of the drag history which features a peak drag coefficient and the first drag minimum due to shock focusing. However, due to the propagation of a rarefaction from the ablation front in our laser-driven system, the drag coefficient features a second minimum associated with laser shut-off. This indicates that rarefaction stretching will contribute to the unsteady drag coefficient as time progresses.

The mass-averaged particle velocity is shown in Fig. 6(c). The velocity history displays three distinct phases: an acceleration phase that occurs over a time scale of τ_s, followed by a phase where velocities tend to level off, and finally a deceleration phase, where the particle velocities decrease significantly due to rarefaction stretching.

We observed roughly fourfold streamwise compression of the particle along with significant deformation, as shown in Fig. 7(a). By applying the modified ideal gas EOS as described in Sec. II, we were able to mitigate compression and deformation of the particle compared to that resulting from using IONMIX EOS.

As the IS propagates through the particle, the particle is compressed in the streamwise direction. The temporal variation of the length (i.e., streamwise) and width (i.e., transverse) of the evolving particle interface is plotted in Fig. 7(b) to quantitatively characterize the particle evolution. The length of the evolving particle decreases quickly at the early stages due to shock compression and reaches a minimum value. At the intermediate stages from t′ = 1 to 3, the length gradually decreases. In tandem, particle width increases at the early stages. This is caused by limitations in keeping both the particle and the host at equilibrium pressure, as discussed in Sec. II. During the intermediate stages, we observe modest growth in particle width. At the later stages, rarefaction catches up to the flow ahead, causing an increase in both the length and the width of the particle. We should remark here that although the particle length reduces by ∼25% and width increases by ∼30% by $t′=3.5$⁠, the particle remains intact without any rollups or interface distortion.

If we observe the interface evolution of the particle modeled with IONMIX EOS, we see a greater reduction in particle length due to shock compression which is inconsistent with the shock Hugoniot data shown in Fig. 3(b). This is caused due to material being modeled with lower γ (i.e., lower sound speed) that resulted in ∼4× streamwise compression of the particle as compared to ∼1.5× with the modified model. We should also note that the width of the particle modeled with the IONMIX EOS slightly reduces during shock propagation. This is due to the unrefracted shock in CH traveling faster than the TS inside the particle. This creates a pressure difference across the particle interface from the shocked CH into the inside of the unshocked particle in the lateral direction as the TS travels through the particle. Such effects were mitigated using the modified EOS model for materials. Figure 3(b) shows how the simulated Hugoniot compression of Al particle using the modified model matches the experimental data.

Figure 8(a) provides the time histories of mass-averaged particle velocity for three metal particles: Al, Ti, and W. The lighter Al particle accelerates more rapidly and to a higher maximum velocity than Ti and W particles during the shock–particle interaction. For instance, the Al particle accelerates to a velocity of ∼4 μm/ns that corresponds to ∼91% of the surrounding shocked fluid velocity. Ling et al.,¹³ from their analytical study, showed that the velocity gained by the particle from only the pressure gradient force scales inversely with particle-to-host density ratio. In the previous works, for the case of Al particles in air, the particle-to-host density ratio is $O(103)$⁠. This results in lower velocity gain in the particle from the shocked air. In our study, the particle-to-host density ratio is O(1). Hence, the gain in particle velocity from the medium during and after the shock interaction was much larger than that in the case of Al in air. To quantify the velocity gain for different cases studied here, we calculate the velocity transmission factor α as the ratio of the peak mass-averaged particle velocity after the initial shock interaction to the shocked surrounding velocity. In Table II, α ranges from 0.29 to 0.96 for particles of different density.

For the same initial particle density, we observe in Fig. 8(b) that the smaller particle accelerates quickly to reach the peak velocity.

We discuss the dependence of the velocity transmission on the particle acoustic impedance $γρP$⁠. The calculations were performed with particles of same initial density and pressure but different acoustic impedances by modifying γ. Figure 9(a) compares the peak velocity attained by particles as a function of γ immediately after the emergence of the TS at the trailing edge of the particle. Note that the maximum particle velocity (black cross) decreases with increasing acoustic impedance. In particular, a fourfold increase in γ (i.e., γ = 20) results in ∼40% reduction in the peak particle velocity. Despite the differences during early times (⁠$t′≤1$⁠), however, Fig. 9(b) shows that the mass-averaged particle velocities attained long after the shock interaction look identical, resulting in a similar α as shown in Table II. To study an effect of γ on the shock velocity, Fig. 10 plots density field showing the shape of the TS inside the particles. For all the cases studied, the TS is convex in shape and runs ahead of the unrefracted IS. Furthermore, the TS propagates faster but imparts lower material velocity in particles with increasing γ [see also Fig. 9(a)]. These observations, more importantly, provide evidence of how modification in γ controlled the stiffness of the material in our simulations. However, its effect to the bulk particle motion is seen for a very short time [i.e., O(τ_s)]. Eventually, particles of the same density but varying acoustic impedance attain the same peak velocity. Therefore, once the shock completely traverses the particle, the particle-to-host density ratio determines the inviscid peak terminal velocity of the particle.

A numerical investigation of laser-driven shock propagation through an isolated particle embedded in a plastic target is presented using the radiation-hydrodynamics code FLASH. The predicted evolution of the particle modeled with IONMIX EOS did not reproduce the experimental shock Hugoniot. Hence, we developed a technique to implement a modified form of ideal gas EOS to model the materials and study the dynamics of the embedded particle. The simulated shock Hugoniots of multiple materials, modeled using this technique, compared well with experimental data. We then examined the flowfield and observed that the wave interactions were highly unsteady to allow the particle pressure to equilibrate within a few τ_s. We also demonstrated that the unsteady drag coefficient for the particle featured a peak drag due to an unsteady interaction with the transmitted shock and a drag minimum due to shock focusing at the rear end of the particle. However, unlike previous studies performed without laser drives, the particle drag coefficient featured a second minimum due to rarefaction stretching associated with laser shut-off. Furthermore, to quantitatively characterize the particle deformation, we plotted temporal variation of length and width of the deforming particle. Although a ∼30% lateral expansion and ∼25% streamwise compression is observed, the particle maintained integrity without any rollups and significant interface distortion. We then conducted numerous simulations and investigated the particle response for a range of particle densities, sizes, and acoustic impedances. The results revealed that lighter particles, such as Al, gained significant momentum up to 96% from the shocked CH, compared to 29% in the case of heavier W. Finally, we studied the effect of particle acoustic impedance on the bulk particle response. Despite differences observed in the early stage of shock interaction, the acoustic impedance did not have an effect on the peak particle velocity. This also identified particle-to-host density ratio as a dominant factor in determining the inviscid terminal velocity of the particle.

Time scale analysis in previous works have pointed out that the shock–particle interaction time scale could be of the same order as the viscous time scale, particularly for condensed-matter systems.⁴ Hence, viscous effects coupled with rarefaction stretching effect could be important for particle drag calculation in the intermediate to later stages of shock interaction. To this end, future work should include viscous models in the simulations to accurately calculate the particle response in such systems. Finally, preheat effects should be negligible due to relatively low-drive conditions studied in this work. Nevertheless, we hope to extend our modified EOS model in the future toward a much higher laser drive and provide a temperature-dependent γ to the model.

We thank Bertrand Rollin for valuable advice during the early stages of this work. We also thank three anonymous reviewers for their thorough reading and helpful suggestions. This work was performed under the auspices of the U.S. Department of Energy under Grant No. DE-SC0019329 within the joint HEDLP program and was supported by the Laboratory Basic Sciences program administered by UR/LLE for DOE/NNSA. H.A. was also supported by U.S. DOE Grant Nos. DE-SC0014318 and DE-SC0020229, NSF Grant Nos. PHY-2020249 and OCE-2123496, U.S. NASA Grant No. 80NSSC18K0772, and U.S. NNSA Grant Nos. DE-NA0003856 and DE-NA0003914. J.K.S. was also supported by NSF Grant No. PHY-2020249 and NNSA Grant No. DE-NA0003914. The software used in this work was developed in part by the DOE NNSA- and DOE Office of Science-supported Flash Center for Computational Science at the University of Chicago and the University of Rochester.

The authors have no conflicts to disclose.

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