DPDE-based mesoscale simulations of shock response of HE composites

The dissipative particle dynamics with energy conservation (DPDE) method is extended to simulate the shock response of high energetic (HE) materials at micron length scales. The symmetrical impact of an RDX impactor and target plates with 1μm diameter spheres is simulated at planar impact velocities of 208 m/s and 876 m/s with a Lennard-Jones-like potential, dissipative, and random forces, and artificial viscosity force between particles. The in situ shock quantities were obtained using Hardy's averaging. In situ longitudinal stresses from simulations were 0.84 and 3.82 GPa. Values from the literature are 0.81 and 2.92 GPa at the two impact velocities, respectively. The uniaxial strain condition was predicted with equal lateral stresses and negligible shear stresses. The higher stress value at 876 m/s may be due to lack of inelasticity in the interparticle potential. The temperature increases of 5.5 0K and 93.7 0K, respectively, were predicted assuming dissipation of a fraction of the potential energy. It is concluded that the DPDE method holds promise for a unified computational framework for multi-scale simulations of HE.


Introduction
The dissipative particle dynamics with energy (DPDE) conservation method [1,2] may successfully simulate the shock response of condensed phase high energy (HE) materials due to the built-in conservation of mass, momentum, and energy. Developed as the isothermal dissipative particle dynamics (DPD) method and used at molecular scales for gases and liquids [3,4], addition of particle internal energy (DPDE) as a variable allows the temperature variations needed for shock propagation simulation without the need to solve heat transport equations. This method has been extended to coarse grained molecular scales for HE [5] and to nanoscales for polymers [6]. It has yet to be applied at the micron and sub-micron length scales for solids. The shock response of HE, leading to detonation, has a transitory state where (a) the solid phase undergoes deformation, failure, and heating at heterogeneous mesoscale [7]; (b) a decomposed state with energy release and heating; and (c) a gas phase with energy release and heating. Each of these phases coexists with different length scales. Finite element methods have been used for inhomogeneous mesoscale simulation of solids [8], and fluid dynamics methods are used to study detonation under the detonation shock dynamics (DSD) field. No computational tools exist, to the best of our knowledge, that can model and simulate the coupled response of the heterogeneous transitory state with a mixture of phases and length scales.
With the recent focus on particle methods for condensed phase HE [9], the DPDE method may be better suited to this task, given that, on the one hand, it can be extended to the molecular length scale from which it originated, and that evolving failures and phases, on the other hand, can be modeled naturally on-the-fly. For example, using DPDE, solid failures (inter-and intra-grain fracture, voids, etc.) can be modeled by dissociating a particle from its neighbors; different phases can be modeled by transitioning to a different set of interparticle potentials; and permeation can be modeled by decomposing a particle into sub-particles. The present work explores whether the DPDE method is applicable at the micron and sub-micron mesoscale lengths of solids while still preserving its core features.

Objectives
The ultimate motivation is to arrive at a unified multiscale, multiphase, simulation framework for HE, while the present objectives are to extend the DPDE method to simulate the shock response of materials at a micron length scale, determine simulation parameters, retrieve in situ averaged shock response parameters for comparison with experimental data, and quantify the in situ shock response for a given impact loading.

DPDE Methodology and In Situ Averaged Shock Quantities
For a comprehensive description of the DPDE method, the reader is referred to previous work [2]. Key features of extending this method to the micron length scale and averaging in situ quantities, however, are presented below. It is assumed that the stochastic fluctuation-dissipation theorem is applicable at the micron length scale. Analogous to the molecular dynamics (MD) method, the DPDE method simulates the shock response of materials by integrating a set of equations-of-motion, where the forces on particle i from neighbouring particles, j, are given as, The force components, shown in brackets above from left to right, are the interparticle potential, dissipative, random, and external forces. The left superscript, t, denotes time. The last term represents the artificial viscosity force component used in the DPDE method for the first time. A Lennard-Jones like potential is adopted in this work from [10] and is given as, 1 1 where A is the cross-sectional area of the particle; α, m, n are the material parameters; and r 0 is the equilibrium separation between particles. Denoting the position vector of particle i by r i and the position vector from particle j to particle i by r ij = r i -r j , the interparticle distance t r ij represents the particle separation. The unit vector at time, t, is t n ij = t r ij / t r ij . Dissipative and random forces are, where γ ij ,ω D and a ij , ω R are the amplitudes and weight functions of the dissipative and random forces respectively, p i is the momentum vector of particle i with p ij = p i -p j , m ij (= m i = m j ) is the particle mass, ς ij is a random number from the standard normal distribution with zero mean and unit variance, and t ∆ is the time step. The final part of equation (3) represents the relationship between ω D and ω R from [2] in which r c is the interaction range for the dissipative and random interactions. The artificial viscosity force, as a combination of linear and quadratic terms, is given as, where 1 C and 2 C are constants. The stochastic and deterministic integration of the equation-of-motion (1) is carried out using a modification to scheme [2], given in condensed form as, 2 2 0.5 0.5 where C i F is the sum of the interparticle and artificial viscosity forces. Integration is carried out for j > i and the negative of the stochastic incremental velocities are updated for particle j. Finally, the particle position is updated as, The particle internal energy (u) calculation is performed (as in equation (21) of [2]) for the conductive and mechanical contributions, while the interparticle potential, equation (2), is integrated to determine potential energy. The internal temperature i θ of particle i is calculated from the The in situ averaged shock quantities are calculated following the Hardy method [11] in current coordinates using a 3D cosine function as the localization function Δ.

Simulation Procedure
The shock response of RDX was simulated for an assumed planar symmetric plate impact experiment. Figure 1(a) shows the simulation representative volume element (RVE) of a 20 μm x 20 μm crosssection with a 50 μm thick RDX plate impacting a 100 μm thick RDX sample. Both the impactor and sample were modeled as regularly stacked 1μm diameter spheres with a center-to-center equilibrium distance of 1 μm. The total mass of the impactor and sample was lumped equally among their particles. The longitudinal (X-) impact velocity to the impactor particles and 300 0 K temperature to all particles were assigned as initial conditions, while the lateral velocities were constrained for particles at the RVE boundaries as boundary conditions. Instead of using a contact algorithm, the interaction between the impactor and the sample particles was analyzed using the interparticle potential-treating the two as one body. For this reason, results are presented for the time when the release wave from the sample's free surface reaches the impactor-sample interface. Simulations were carried out at impact velocities of 208 m/s and 876 m/s to explore the DPDE method for an elastic particle velocity of 104 m/s and an inelastic particle velocity of 438 m/s, respectively, in the (100) direction of Hooks et al. [12] data. The free surface velocity was calculated by averaging the longitudinal velocity of the particles over a 4 μm x 4 μm surface area, shown in red; the in situ shock parameters were obtained by averaging in a 4 μm x 4 μm x 2 μm volume, shown in green in the zoomed figure 1(b). The averaging volume was centered with the cross-section and longitudinally followed a specified particle center.

Materials Properties
The properties of materials used in the simulations are summarized in Table 1. The mass density and specific heat values are in the range of experimental values [13,14]. The parameters of the dissipative and random forces were adapted from Lisal et al. [2]. The interparticle potential parameter α [10] was calculated for RDX by adopting an average Young's modulus of 19.29 GPa [15]. Simulations considered the particles' first nearest-neighbor at 1 μm, the second nearest neighbor (face diagonal) at 2 1/2 μm, and the third nearest neighbor (body diagonal) at 3 1/2 μm. A constant cutoff distance of 1.5 μm from the equilibrium position in the attraction (or tensile) phase was used for all three neighbor types. The (m,n) parameters of the interparticle potential given in equation (2) were varied as (2,4), (2,3), and (1,2) reducing the slope of the curve (or, a reduction in the materials stiffness) in the repulsion (or compressive) phase. The value of (2,4) is the same as that used by Tang et al. [10] for HMX. The use of the same α and (m, n) parameters resulted in different potentials for the second and third nearest neighbor types.

Results and Discussion
In the initial simulations, tolerance issues and optimal integration time step (Δt) determination were targeted. Simulations were stable and yielded identical results for Δt ≤ 0.14x10 -9 seconds. All results presented below are from simulations with a constant Δt of 0.14x10 -10 seconds. Figure 2 compares the free surface velocity profiles at 208 m/s impact velocity simulated with and without artificial viscosity. Similar results were obtained at an impact velocity of 876 m/s. As found in shock propagation simulations using finite difference/element   methods, the linear and quadratic artificial viscosities effectively damp oscillations in the free surface velocity and other in situ shock quantities predicted by the DPDE method. Figure 3 shows the average free surface velocities at the two impact velocities for the three sets of (m, n) parameters. The simulations predicted a steady peak shock state, independent of the potential parameters, in qualitative agreement with experimental data [12]. The velocity profiles agree with the shock propagation theory. The free surface velocity increases on arrival of the shock wave from the impact surface, attains peak, and remains constant representing the shock state (plateau) under the combined influence of the shock and release waves reflected from the free surface. The velocity begins decreasing to zero after the rarefaction wave from the impactor free surface arrives at the target's free surface. Shock velocity was calculated from the time of the half-rise of the free surface velocity. Shock velocity diminishes and shock rise time increases as the interparticle stiffness drops from (2,4) to (1,2) parameter values. Figure 4 shows the in situ average mass density, longitudinal stress, and lateral stresses at 876 m/s impact velocity with positive values representing compression. The stresses begin to increase on arrival of the shock wave; they peak, and then remain constant (plateau) until the arrival of the release wave. Due to the proximity of the averaging domain to the impact surface, the rarefaction wave from the impactor's free surface reduces the stress to zero. The stresses become tensile (negative) on arrival of the release wave from the target's free surface. The simulations predicted nearly constant peak values of the shocked mass density increasing with the reduction in the interparticle stiffness. The reduced shock speed and increased compression are associated with lower in situ longitudinal stress.
The method did not predict the measured [12] two wave profiles at 876 m/s attributed to inelastic (plastic) deformation, shown in Figures 3 and 4. The onset of inelastic deformation reduces the slope of the longitudinal stress vs. strain curve and, in turn, reduces the wave speed (termed inelastic wave speed). The data shows two wave profiles due to a slower inelastic wave following the elastic wave. The single wave profile represents either the elastic wave at lower impact stresses (as in the case of 208 m/s) or the elastic wave overtaken by the inelastic wave at higher impact stresses. The Lennard-Jones like potential used in the simulations has a continuously increasing force-displacement slope without any reduction expected at the elastic-inelastic transition. This may be why the measured two wave profile was not predicted at 876 m/s. The peak longitudinal stress values at the two velocities for the (2,3) parameters were 0.84 and 3.82 GPa compared to 0.81 and 2.92 GPa, respectively [12]. The over-prediction may further indicate that the elastic response at 876 m/s was predicted by the simulation because inelasticity limits deviatoric stress, which, in turn, lowers longitudinal stress below the elastic value.
As an approximation, 90% of the potential energy was dissipated as heat to explore heat dissipation and conduction in the DPDE method. Table 2 lists calculated temperature and other in situ quantities in the peak shock state. Moreover, the reduction of in situ longitudinal and lateral stress difference as the (m,n) parameters varied from (2,4) to (1,2), shows the strength of the prediction capability of this method. The stress difference is the material's strength under the uniaxial strain condition as per the von-Mises plasticity theory.

Conclusions
The present work successfully extends the DPDE method to simulate the shock response of materials to the micron scale. This method predicted the free surface velocity and in situ average shock response in qualitative agreement with the uniaxial shock wave propagation theory. The longitudinal stress value at the two impact velocities was close to the literature value, however quantitative agreement cannot be established due to the lack of inelasticity in the interparticle potential. The heat dissipation and conduction built into DPDE was verified through an approximation of dissipating a fraction of the potential energy. The calculated temperature rise, therefore, should be treated as a model calculation. Nevertheless, this first use of artificial viscosity led to prediction of a smooth shock wave profile. The present work is being extended to incorporate a contact algorithm to model the impactor-sample interface and interparticle potential with inelasticity. It can be concluded that the present DPDE implementation has potential to link the continuum scale to the micron and sub-micron length scales, preserving the emphasis on interparticle interactions based on potential rather than continuum phenomenological models.