A High Resolution Simulation of a Single Shock-Accelerated Particle

Particle drag is widespread in physical phenomena with parameter scales that span several orders of magnitude. For example, Crab Nebula particles can be seen in motion as electromagnetic radiation and gravity interact with them. Droplets in chemical weapons are deformed and broken up by a blast wave, such that it is difficult to determine how they are dispersed. Rotating detonation engines typically introduce droplets ahead of a detonation wave to break them up, evaporate them, and burn them to perpetuate their thermodynamic cycle. All situations described here benefit from a better understanding of the shock-driven multiphase instability (SDMI). The SDMI is a hydrodynamic instability that is studied by several groups [1,2], and exists due to the nonzero equilibration time of the solid or liquid particles. The equilibration time is a function of the particle and postshock fluid's material properties, and will give a distance the particle will lag by. Fully understanding the drag experienced by shock-accelerated particles is essential to understanding the fundamental physics and characteristic parameters of the SDMI, and in closing the larger problem of mixing in these situations.

The study of particle-laden flows is a rich area of study at the time of this paper. Many applications involve a cloud of particles, and recent papers focused on particle–gas interactions in a cloud of densely packed particles. Osnes et al. [3] has conducted a rigorous three-dimensional (3D) numerical study of a shock wave passing over a randomly distributed cloud of particles using large eddy simulation. Mehta et al. [4] has carried out a similar numerical study of a cloud of particles but without the consideration of viscosity. Others have studied vorticity production [5], turbulent stresses [6], and closure terms for the Navier–Stokes equations [7] in densely packed particle clouds. While these works have taken on the more challenging problem of resolving particle-to-particle interactions, there is still little highly resolved data capturing the particle material response in the simpler, single particle shock interaction.

For these many particle problems, it is still important to understand the interaction between a single particle and the external fluid. The governing equations for this interaction are a coupled set of nonlinear partial differential equations. In SDMI simulations, typically thousands or millions of particles are present [2], and these equations cannot be used to resolve a cloud of particles in flow. This is due to the computational cost associated with resolving the large range of relevant length scales. A simplification of these equations may be attributed to the scientists Maxey and Riley for their rederived form of the Basset-Boussinesq-Oseen equations [8]. This reduces the complexity of the equations to be solved to a single nonlinear integro-differential equation, by assuming the particle to be rigid and the fluid to be Newtonian.

This model is accurate for applications where compressibility is not a strong factor, and only depends on the Reynolds number (⁠$ReD=||vP−vPS||ρD/μ$⁠), where $vP, vPS$⁠, ρ, and μ are the particle velocity, the postshock velocity, the postshock density, and the postshock dynamic viscosity, respectively. However, for particles in high speed gas flow, the incompressible models will under-predict the drag as stagnation pressure increases with the Mach number, along with many other nonlinear effects. One model that incorporates the Mach number into its parameter space is that of Parmar et al. [9]. This model uses both the Reynolds and Mach numbers as inputs to yield the drag coefficient of the particle via a method of interpolating polynomials. Both of these models are considered in this work as a comparison to simulation results.

In Eq. (2), μ and μ_S are the air viscosities evaluated at the freestream and particle surface, respectively.

Experiments showed that, for a given set of flow and particle parameters, drag correlations of Parmar and Clift-Gauvin [10,17] severely under-predict the measured drag [18]. This was measured by tracking particle positions with a high-speed camera as the particle traversed down a shock tube. This motivated the application of a more complete set of physics to simulate this problem and determine if the results agreed with the models. This set of physics included the application of two-dimensional heat and momentum (through viscous effects) transfer and included the effect of particle motion (the particle was allowed to translate through the flow). It was found that the numerical solution agreed well with the models under the conditions of the experiments but failed to explain the disagreement between experiments and the models. The details of the numerical solver will be discussed in Sec. 2.

FLAG was the environment selected to simulate the shock-accelerated particle in this study, and is an arbitrary Lagrangian Eulerian (ALE) code. ALE codes provide the material-capturing abilities of Lagrangian codes, while also providing the some of the benefits associated with Eulerian codes (e.g., efficiency). They are executed by first solving the conservation laws while tracking the material, such that the mesh becomes deformed, then remapping the mesh using an optimization method. Remapping the mesh means to update the mesh but without many of the required steps to completely remesh the domain. This step is necessary to prevent the mesh from tangling in regions of high shear strain rates and to improve the quality of the solution. The mesh optimizer that was found to perform best, considering efficiency and effectiveness, was a condition number-based method. This method attempts to remap the mesh such that it returns cells that have been sheared back to rectangular. This is done by minimizing a global condition number objective function of the mesh, which is an average condition number of the Jacobi matrix of each element in the mesh. For more details to the implementation of this method into FLAG, refer to Dyadechko et al. [19]. It is one of the more robust mesh optimizers available for this application with minimal iterations necessary. Material is advected through the finite volume boundaries in this remapping step because the code is now leaving the Lagrangian step of its operation. FLAG uses the Flux-Corrected Transport method of Boris and Book [20] to yield second-order accuracy in space. This process is then repeated in a time-marching manner.

The viscosity of the air, μ is calculated as a function of the temperature, T given a reference viscosity μ₀ measured at a given temperature T₀. All capabilities described here make FLAG ideal for tracking the motion of a shock-accelerated particle while resolving the heat transfer and mechanical response of the particle due to the shock.

In addition to the viscous diffusion capabilities added to FLAG, a thermal diffusion solver has been added as well. This was done so that heat would transfer into the particle from the flow, which results in a thermal boundary layer. The resultant decrease in temperature in the gas immediately surrounding the particle drops the viscosity as well, in accordance with Sutherland's law (Eq. (3)). In the case where heat transfer is not considered, the no-slip condition for the flow at the surface of the particle forces there to be a region of large shearing with no mass advection. The viscous work done causes an artificial increase in temperature, and therefore viscosity, at the surface which drastically increases drag. The heat diffusion solver is the only one in FLAG that is implicit, which implies that stability will not be affected by imposing the heat transfer solver to the list of physics being considered.

While FLAG has a built in mesh generator, Ingen is capable of generating meshes that are more intricate and often better suited for our simulations. Ingen is a Python library that has the ability to construct rectangular-structured, radial-structured, fully unstructured, and combination meshes from a command line interface. Ingen, like FLAG, is also a code that is under development at LANL. Since its start, its operation has been overseen by the LANL Setup team. The domain can be optionally subdivided into regions where different meshing rules may be imposed. Voronoi tessellations and Delaunay triangulation are supported as unstructured meshing rules. The Delaunay triangulation works on a set of discrete points when no point lies inside the circumscribed circle of any triangle. The Voronoi tessellation is the mathematical dual of the Delaunay triangulation of a point set. Furthermore, functions for dendritic derefinement are supported within Ingen. This work makes use of rectangular-structured, radial-structured, and Voronoi tessellation meshes to resolve the domain. Dendritic derefinement is implemented in both the rectangular and radial structured meshes. The finalized mesh is exported as an ×3D file so that it can be imported into FLAG.

Before simulating the shock-accelerated particle, the computational domain had to be built. It was decided that the three dimensional problem was to be reduced to a two-dimensional axisymmetric one to reduce the computational burden associated with solving the equations in three dimensions. In this decision, the assumption of perfect symmetry was imposed, however the Reynolds number for the problem was low enough (⁠$ReD$ = 40.8) that turbulence is not expected to onset, and therefore three-dimensional effects are not expected to be significant. This assumption may be checked in future work.

Kinematic measurements of the 4 μm particles were taken by the Extreme Fluids team [23] with given ambient conditions of the lab, the shock strength, and the particle properties also provided. Postshock properties are computed with the shock tube relations, given by Anderson [24], which derive from conservation of mass, momentum, and energy for an ideal, calorically perfect gas. The preshocked and postshocked gas properties are listed in Table 1. The particle is modeled as an isotropic solid, with the yield strength given as an input to observe the deformation, and possibly failure, of the particle.

The domain was built as rectangular with three wall boundary conditions and one inflow boundary condition. Because it is difficult to implement an outflow condition for ALE codes that do not reflect shocks, the domain had to be extended in the streamwise direction such that the reflected shock would not interfere with the accelerating particle. The distance was calculated using the initial and reflected wave speeds so that the particle could translate at least a full diameter into the postshock flow before reshock occurred. The downstream domain was derefined by a factor of sixteen to reduce the computational cost without losing important physics of the problem. The width of the domain was constrained to be ten times the radius of the particle in an attempt to avoid reflecting boundary effects. These effects are in the form of an artificial shock, which are large enough to affect the simulation result. The fluid domain was meshed to be mostly square-structured with a transition region from the particle that is resolved with a Voronoi tessellation and a radial-structured mesh. These descriptions are illustrated in Fig. 1.

A Delaunay triangulation pattern was attempted to resolve the transitional region, but when the triangular cells deform as the simulation runs, the volume of each cell cannot be held constant by altering another degree-of-freedom, and the pressure in the cell artificially changes as a result. The mesh constructed with Ingen used both structured and unstructured components. This was done with the intent of building a mesh that would be robust as resolution increased while the meshing rules stayed the same. The sound speed in the solid region is much larger relative to the fluid region, so the CFL limitation typically occurs within the particle. Small cells in the interior of the particle will drive the time-step to be extremely small values unless intentional action is taken to prevent this. The mesh in the particle is dendritically derefined by a factor of four to prevent cells from shrinking and resulting in an un-necessarily small time-step. A small region of high resolution is present along the rim of the particle to have mesh continuity at the boundary for calculation of the fluid stresses at the boundary. The particle is resolved with a rectangular-structured mesh so mesh stiffeners may be used in conjunction with the hydrodynamics to treat stress deviators.

Here, M is the postshock Mach number, $M=vp/camb$⁠. A general rule is that when the Knudsen number is much less than 1, typically less than 0.01 [25], the continuum assumption is a good approximation. Conversely, when the Knudsen number becomes large, free-molecular flow is a good approximation. When the Knudsen number approaches 1 from either side, the flow will behave in a less predictable manner, and adjustment models are oftentimes imposed to reduce the complexity to free molecular or continuum flow. The Knudsen number for this case is at 0.014 and is at the bounds of the continuum flow regime. It is still a strongly debated topic on which value the continuum assumption fails at; however, it is observed that the continuum assumption for particles is still appropriate for larger Knudsen numbers [26] and is argued here to be an effective model.

To guarantee a mesh-independent solution, a resolution study was carried out. Convergence of the solution was verified for both kinematic and thermal measurements. Ingen was utilized to develop four different mesh resolutions for the simulation to run on. The refinement of the mesh was measured by the number of faces present on the simulated particle's half-circumference. All meshing rules were kept constant otherwise to keep a fair comparison between cases. The computational cost of each of the four cases are given in Table 2, which shows that for accuracy and associated computational cost, the 645 case is the strongest. The resolution was scaled up with the length scale of the elements to make comparisons about the timing statistics and accuracy of solution. The simulations were run out to a stopping time of five hundred nanoseconds and the results were compared, as transient effects have become negligible at this simulation time. For the sake of comparison to other studies, the time has been nondimensionalized as $τP=Dt/(Mcamb)$⁠, where c_amb is the sound speed at the ambient conditions, $camb=γRTamb$⁠, γ is the ratio of specific heats, R is the gas constant for the ambient species, and T_amb the ambient temperature.

Here, A is the area of the computational cell, and r is the variable ranging over the radial direction of the particle. This is repeated for every time-step in each of the four resolution study simulations, and the average particle temperature is plotted against simulation time. The results are shown in Fig. 2. Interestingly, from the coarsest mesh resolution, the average particle temperature is relatively converged to the same answer.

Clear convergence can be seen in the position results as mesh resolution moves from the 101 faces per half-diameter case to the 645 one, as illustrated in Fig. 3. The relative error in the particle position between the coarsest two cases is over 50%; however, the relative error between the most refined cases is less than 3%. This demonstrates the law of diminishing returns, and is consistent with the second-order accuracy in time and space. It is assumed here that the solution is reasonably converged with the case of highest resolution. The limiting factor clearly is set by the particle kinematics agreement and not by the average particle temperature. The 645 faces case was chosen to be the one that was run out in time, due to the minimal relative error between it and the 432 faces case, and its reasonably affordable computational cost.

Here, $α=κ/(ρpc)$ is the thermal diffusivity of the particle, c is the specific heat of the particle, $h¯$ is the surface-average heat transfer coefficient, and $T∞$ is the temperature of the freestream. N is the number of subdividing shells the particle is divided into in the one-dimensional model. Convergence to the final answer occurs at 200 subdividing shells. As parameters are matched, and convergence of the lumped capacitance method of shells is realized, the convection correlation is isolated as the model to test FLAG against. This correlation is compared to the values of the FLAG particle in simulation time. Temperature gradients are found to appear in the particle over time. This is expected as the computed Biot number of the particle is 0.6. Thus, the average temperature of the particle is computed over time from the simulation and compared to that of the one-dimensional model. It is seen that the model is in good agreement with the data provided in FLAG (Fig. 4). Given that compressibility is not taken into consideration with the model, and that heat transfer correlation of Whitaker has experimental scatter of 30% [12], the conduction heat transfer physics in FLAG is within the range of validity for the model.

From Fig. 5, the thermal and momentum boundary layers may be compared qualitatively. Note that the boundary layers develop as expected, with the ratio of the thickness of the momentum to thermal boundary layers roughly matching the value of the Prandtl number (0.71). This provides some confidence that, with the thermal model agreement with simulation results and qualitative agreement seen between thermal and momentum boundary layers, the models and simulation results for the particle kinematics should be in agreement. This will be covered in Sec. 3.3.2.

Numerical integration is carried out by decomposing the equation into two first order ordinary differential equations (ODEs) and integrating with the Runge–Kutta based ODE solver built into Scipy [27]. The comparison of the model to the simulated particle is shown in Figs. 6 and 7. The model of Clift-Gauvin slightly under-predicts the simulation, while the simulation is seen to agree well with the Parmar model. The positions yielded by the models are seen to agree with one another within 5%. The agreement of the models suggests that compressibility for the given Mach number of roughly 0.4 is not a strong factor, which is important in the validation of the heat transfer model, as compressibility is not considered in its parameter space. However, compressibility may be the difference between a fairly close answer and the correct answer. The previous models were constructed using data taken from stationary particles in both simulations and experiments. Here, our simulations allowed the particle to move, yet the results were still in agreement with these models. The source of the deviation in particle position and acceleration originates during the time the shock wave is passing over the particle, before significant motion of the particle occurs. Thus, the motion of the particle at early times does not seem to have a strong effect on the particle drag. While we neglected the viscous unsteady component of Parmar's model (for simplicity), this force may be able to account for the remaining deviation between our simulation results and the model.

Given that the models agreed well with the simulation results, it was desirable to determine what might have caused the experimental results to disagree. One possibility is that the particle size reported in the experiment may have been incorrect. To determine what particle size would have brought the experiments and models into agreement, the particle diameter was treated as a variable and numerically solved for using the drag coefficient correlation of Parmar et al. [17] and the position data from the Extreme Fluids team [23]. Using the Maxey-Riley equation, and assuming only the forces of quasi-steady and inviscid-unsteady drags acting on the particle, a nonlinear first order ordinary differential equation appears, and is given by Eq. 12.

Rather than solving for the velocity of the particle as a function of time, the particle diameter may be treated as the variable if given a position at a postshock time. In this work, the numerical time integration is done with a variable time-step Runge–Kutta method, and the numerical solution for particle diameter is done with a Newton–Raphson root-finding method [27]. The particles were reported to be nylon [18], fixing their density at 1140 kg per cubic meter. This information, coupled with the postshock properties of the flow given in Table 1, lead to the conclusion that the particles must have a diameter of 0.53 μm to agree with the model proposed by Parmar et al. and Clift and Gauvin. Unfortunately, this value is well outside the regime of validity for continuum flow, and FLAG is exclusively a continuum hydrocode.

The viscous and heat transfer physics solvers of the hydrocode FLAG have been shown to agree with validated models in their respective fields. This has been done by executing a direct numerical simulation of a single shock-accelerated particle and comparing kinematic and thermal results to experimentally validated one-dimensional models. Because more relevant physics have been considered in this work, such as the material response and dynamics of the particle from the shock, it may be concluded with better certainty than before that the drag models perform well within their parameters' regions of validity. The agreement may additionally suggest that continuum applies well for Knudsen numbers slightly larger than previously suggested, for both models and simulations in the field of particle-laden flows. Drag models considered here have been developed without acknowledging the acceleration of the particle into the flow, but results obtained via the simulations suggest that the motion of the particle into the flow does not alter the acceleration. The solution of the kinematics would suggest that, for time-scales of the order of microseconds for micrometer-sized particles, viscous unsteady effects of the shock need to be accounted for in models. Furthermore, the convection correlation model of Whitaker has been shown to agree reasonably well with the heat transfer simulated with FLAG, considering the experimental scatter of the data the model was built with. These results have been found in the pursuit of the initial goal of this work: to replicate the experiments from LANL. With the agreement seen between the one-dimensional models, a possible explanation is that the particles that were previously used in the experiments were smaller than reported. We look forward to future experiments from the Extreme Fluids team that may shed light on this question.

This work was funded by contract number 89233218CNA000001 and subcontract number 527153 of Triad National Security. The authors would like to thank the Extreme Fluids Team at LANL. In particular, Kyle Hughes and John Charonko for their ongoing guidance and collaboration along the way. Approved for unlimited release: LA-UR-20-26838.

• Los Alamos National Laboratory (Grant No. 89233218CNA000001; Funder ID: 10.13039/100008902).