High-frequency mechanical excitation has been shown to generate heat within composite energetic materials and even induce reactions in single energetic crystals embedded within an elastic binder. To further the understanding of how wave scattering effects attributable to the presence of an energetic crystal can result in concentrated heating near the inclusion, an analytical model is developed. The stress and displacement solutions associated with the scattering of compressional plane waves by a spherical obstacle (Pao and Mow, 1963, “Scattering of Plane Compressional Waves by a Spherical Obstacle,” J. Appl. Phys., 34(3), pp. 493–499) are modified to account for the viscoelastic effects of the lossy media surrounding the inclusion (Gaunaurd and Uberall, 1978, “Theory of Resonant Scattering From Spherical Cavities in Elastic and Viscoelastic Media,” J. Acoust. Soc. Am., 63(6), pp. 1699–1712). The results from this solution are then utilized to estimate the spatial heat generation due to the harmonic straining of the material, and the temperature field of the system is predicted for a given duration of time. It is shown that for certain excitation and sample configurations, the elicited thermal response near the inclusion may approach, or even exceed, the decomposition temperatures of various energetic materials. Although this prediction indicates that viscoelastic heating of the binder may initiate decomposition of the crystal even in the absence of defects such as initial voids or debonding between the crystal and binder, the thermal response resulting from this bulk heating phenomenon may be a precursor to dynamic events associated with such crystal-scale effects.

## Introduction

The study of the dynamic stress concentrations induced in an elastic medium by the scattering of waves by an embedded object has an extended history [110]. Ying and Truell's seminal work [1] investigated the specific case of incident compressional plane waves scattered by a spherical obstacle in an infinite elastic medium, with the obstacle taken as an elastic sphere, a rigid sphere, and a spherical cavity. The authors provided an analytical solution for the stress fields generated in the obstacle as well as in the surrounding medium. Pao and Mow [2] expanded and generalized this analytical solution of the stress fields to account for obstacles including a fluid-filled cavity and a rigid spherical obstacle undergoing translational motion. The analysis, specifically for a fluid-filled cavity, was extended to account for obstacles embedded in “lossy” viscoelastic media by Gaunaurd and Uberall [6]. This extended analytical solution can be used to predict the stress fields created in an infinite viscoelastic medium due to an incident compressional plane wave scattered by an embedded spherical obstacle. Additionally, the generation of heat within viscoelastic materials under cyclic loading, including at stress concentrations, has been widely investigated and the time-averaged volumetric heat generation due to the stress field produced by this mechanical loading can be predicted [1121]. Of especial note are the studies of Hinders et al. [21] and Chervinko et al. [17]. Hinders used an analytical stress analysis to develop modal responses for heat generation in a lossy inclusion. Chervinko used similar techniques to calculate temperature fields from two-dimensional finite element calculations of stresses in a viscoelastic medium around an inclusion. In this work, an analytical stress field solution for the scattering of compressional plane waves due to a spherical rigid obstacle within a viscoelastic material is extended to predict the resulting dynamic temperature field through the use of the Fourier law of conduction. Thus, the system thermal behavior is obtained from analytical stress expressions in three dimensions.

This treatment for predicting the heat generation within a particle–binder system under cyclic loading could prove particularly useful in the study of hot spot generation within composite explosives. It has long been suggested that the initiation of detonation of composite explosives is caused by the generation of numerous hot spots induced at locations of extreme stresses, most notably at defects or voids, during the mechanical deformation processes associated with initiation [2225]. The generation of such “hot spots” is difficult to observe experimentally, due to the extreme pressures, temperatures, and short time scales associated with these processes. Recently, ultrasonic excitation has been used to generate hot spots within samples consisting of an elastic binder material with embedded single energetic particles, or crystals [2629]. This experimental investigation of ultrasonically induced hot spots has been proposed as a valuable technique to investigate corresponding processes during detonation events.

In order to better understand the causation and locations of these hot spots occurring in single-crystal systems undergoing ultrasonic excitation and to address whether bonding defects or voids are required to yield sufficient heat generation to cause chemical decomposition, the stress and temperature fields of an idealized model are considered in this work. The stress fields created by the scattering of compressional plane waves at ultrasonic frequencies by a rigid spherical obstacle placed within a viscoelastic medium (where perfect bonding is assumed) are predicted and subsequently used to estimate the volumetric heat generation induced within the sample by the continuous harmonic excitation. This result is subsequently extended through the use of the Fourier law of conduction to predict the temperature evolution and distribution near the embedded crystal. The results indicate that significant heat generation occurs within the binder material near the crystal, which leads to temperature predictions that may reach or exceed the decomposition temperatures of many energetic materials over relatively short time scales. While the mechanism considered here may dominate the initial heating, this process may also be a precursor to other dynamic events occurring at the crystal–binder interface, such as debonding due to the effects of thermal expansion and phase changes.

## Stresses in a Lossy Medium With a Spherical Inclusion Under Compressional Plane Wave Excitation

The analytical solution for the stress fields in an infinite elastic medium with a spherical inclusion subjected to harmonic compressional plane wave excitation has been presented for numerous cases [110]. The case under consideration in this work is that of a movable rigid sphere in an infinite, homogeneous, isotropic, linear viscoelastic medium. For the sake of brevity, only the equations for the wave potentials and the solutions for the stress fields and particle displacement will be given here. Detailed derivations can be found in the works of Pao and Mow [2] and Sessarego et al. [9]. Since only the steady-state response is considered, those derivations for a lossless elastic medium can be easily extended to a linear viscoelastic medium by replacing the real elastic parameters with the corresponding complex material parameters [6,8], which accounts for the viscoelastic losses at the specified excitation frequency.

The utilized rectangular and spherical coordinate systems are shown in Fig. 1. The surrounding viscoelastic medium, or binder, is characterized by density ρ1, complex longitudinal wavenumber $α̃1$, and complex shear wavenumber $β̃1$. The first and second Lamé parameters, $λ̃1$ and $μ̃1$, as well as all other material moduli, can be subsequently calculated from the relations for linear viscoelastic media [30]: $λ̃1=ρ1ω2(1/α̃12−2/β̃12)$ and $μ̃1=ρ1ω2/β̃12$, where ω is the angular excitation frequency. The spherical particle, since it is modeled as rigid, is characterized mechanically only by the density ρ2 and particle radius a. The incident compressional plane wave, as shown in Fig. 1, travels in the positive z-direction in the surrounding medium and is characterized by its longitudinal wave potential

Fig. 1
Fig. 1
Close modal
$Φ̃inc=Φ̃0ei(α̃1z−ωt)$
(1)

where t is the time variable and $Φ̃0$ is the incident wave potential amplitude at z = 0. Since there is no transverse incident wave, the incident shear wave potential is zero.

The stresses in the viscoelastic binder medium are recovered through the application of the boundary conditions at the particle–binder interface, r = a, which require continuity of the displacement component normal to the interface $ũr$ (i.e., the radial component) and the displacement component in the in-plane rotational direction $ũθ$. The stress boundary condition is satisfied by Newton's second law: the force obtained by integrating the stresses over the surface of the particle equals the particle mass times its acceleration [2,9]. It should be noted that, due to the symmetry about the z-axis, there is no displacement in the $ϕ$-direction, and the shear stresses $σ̃rϕ$ and $σ̃θϕ$ are likewise zero.

In order to apply these conditions, the incident wave potential can be expressed in spherical coordinates as a summation [2]
$Φ̃inc=Φ̃0e−iωt∑n=0∞(2n+1)inPn(cos θ)jn(α̃1r)$
(2)
where the Pn are the Legendre polynomials and the jn are the spherical Bessel functions of the first kind. Similarly, the wave potentials for the scattered longitudinal and shear waves in the viscoelastic medium can be written, respectively, as
$Φ̃sctd=e−iωt∑n=0∞ÃnPn(cos θ)hn(α̃1r)Ψ̃sctd=e−iωt∑n=0∞B̃nPn(cos θ)hn(β̃1r)$
(3)

where the hn are the spherical Hankel functions of the first kind, and the $Ãn$ and $B̃n$ are the coefficients which are determined by the boundary conditions. The reader is referred to the work of Pao and Mow [2] for the expressions for the coefficients $Ãn$ and $B̃n$ (written in Ref. [2] without the tilde accents), as given for the case of a movable rigid inclusion. The refracted wave potentials do not need to be taken into account since the particle is modeled as rigid. Moreover, since no deformations are developed within a rigid particle, all of the strain components are identically zero in the region r < a, and the particle simply undergoes rigid body motion along the direction of the incident wave.

The solution for the stresses in the binder material (i.e., in the region r ≥ a), which allows for the translational motion of the spherical inclusion, follows from the wave potentials and can be explicitly expressed, in spherical coordinates, as [2,9]:
$σ̃rr,1=2μ̃1r2e−iωt∑n=0∞Pn(cos θ){Φ̃0in(2n+1)[(n2−n−β̃12r2/2)jn(α̃1r)+2α̃1rjn+1(α̃1r)]+Ãn[(n2−n−β̃12r2/2)hn(α̃1r)+2α̃1rhn+1(α̃1r)]−B̃n[n(n+1)(n−1)hn(β̃1r)−n(n+1)β̃1rhn+1(β̃1r)]}$
$σ̃rθ,1=2μ̃1r2e−iωt∑n=0∞dPn(cos θ)dθ{Φ̃0in(2n+1)[(n−1)jn(α̃1r)−α̃1rjn+1(α̃1r)]+Ãn[(n−1)hn(α̃1r)−α̃1rhn+1(α̃1r)]−B̃n[(n2−1−β̃12r2/2)hn(β̃1r)+β̃1rhn+1(β̃1r)]}$
$σ̃θθ,1=2μ̃1r2e−iωt∑n=0∞(Pn(cos θ){Φ̃0in(2n+1)[(n+α̃12r2−β̃12r2/2)jn(α̃1r)−α̃1rjn+1(α̃1r)]+Ãn[(n+α̃12r2−β̃12r2/2)hn(α̃1r)−α̃1rhn+1(α̃1r)]−B̃nn(n+1)hn(β̃1r)}+d2Pn(cos θ)dθ2{Φ̃0in(2n+1)jn(α̃1r)+Ãnhn(α̃1r)−B̃n[(n+1)hn(β̃1r)−β̃1rhn+1(β̃1r)]})$
$σ̃ϕϕ,1=2μ̃1r2e−iωt∑n=0∞(Pn(cos θ){Φ̃0in(2n+1)[(n+α̃12r2−β̃12r2/2)jn(α̃1r)−α̃1rjn+1(α̃1r)]+Ãn[(n+α̃12r2−β̃12r2/2)hn(α̃1r)−α̃1rhn+1(α̃1r)]−B̃nn(n+1)hn(β̃1r)}+cotθdPn(cos θ)dθ{Φ̃0in(2n+1)jn(α̃1r)+Ãnhn(α̃1r)−B̃n[(n+1)hn(β̃1r)−β̃1rhn+1(β̃1r)]})$
(4)
Finally, the displacement of the rigid particle in the z-direction can then be written as [2]
$ũz,2=ρ1ρ2a[3iΦ̃0j1(α̃1a)+Ã1h1(α̃1a)−2B̃1h1(β̃1a)]e−iωt$
(5)

It should be noted here that, though the binder medium was assumed to be infinite for the purpose of developing the analytical stress solutions, for the solution of the thermal response in Sec. 3.2, the binder will be treated as a large sphere, concentric with the particle, with a convective surface condition applied at the outer radius R. This outer radius at which the thermal boundary condition is applied must be large in comparison with the incident wavelength and particle radius, so that the computation of the stresses in an infinite medium remains a valid approximation.

## Viscoelastic Heating of the Lossy Medium

### Volumetric Heat Generation.

The mechanical energy dissipated in a viscoelastic medium per unit volume can be approximated by the losses in the strain energy density over each harmonic loading cycle [1113]. Under the assumption that the temperature variation is on a much slower time scale than the mechanical loading, this energy dissipated per cycle can be time-averaged to give the rate of energy dissipation. If it is further assumed that the mechanical energy is dissipated entirely as heat, then this rate is equal to the volumetric heat generation
$q=ω2π∫t0t0+2π/ω(σrr∂ϵrr∂t+σθθ∂ϵθθ∂t+σϕϕ∂ϵϕϕ∂t+2σrθ∂ϵrθ∂t+2σrϕ∂ϵrϕ∂t+2σθϕ∂ϵθϕ∂t) dt$
(6)

where the ε** are the strain components and t0 is the initial time for the specified cycle. Note that the real parts of the stress and strain components are used in the integral computation, as the phase lags introduced by the material losses yield the net mechanical energy losses.

By using the linear viscoelastic model employed in Sec. 2, the heat generation in the binder material can be written explicitly in terms of the complex material moduli [30]. The generalized Hooke's law for isotropic media can be applied to relate the strain components to the nonzero stress components, which yields the strain magnitudes and relative phase lags. The volumetric heat generation in the binder medium for this particular loading case is then the sum of four nonzero contributions to the net change in strain energy density
$q1=ω2|Ẽ1|(ξrr,1+ξθθ,1+ξϕϕ,1+2ξrθ,1)$
(7)
where $Ẽ1$ is the complex Young's modulus and the scaled contributions ξ**,1 are given by
$ξrr,1=sin(δrr,1)|σ̃rr,1||σ̃rr,1−ν̃1(σ̃θθ,1+σ̃ϕϕ,1)|ξθθ,1=sin(δθθ,1)|σ̃θθ,1||σ̃θθ,1−ν̃1(σ̃rr,1+σ̃ϕϕ,1)|ξϕϕ,1=sin(δϕϕ,1)|σ̃ϕϕ,1||σ̃ϕϕ,1−ν̃1(σ̃rr,1+σ̃θθ,1)|ξrθ,1=sin(δrθ,1)|1+ν̃1||σ̃rθ,1|2$
(8)
where $ν̃1$ is the corresponding complex Poisson's ratio and the δ**,1 are the phase differences between the respective stress and strain components
$δrr,1=∠Ẽ1+∠σ̃rr,1−∠[σ̃rr,1−ν̃1(σ̃θθ,1+σ̃ϕϕ,1)]δθθ,1=∠Ẽ1+∠σ̃θθ,1−∠[σ̃θθ,1−ν̃1(σ̃rr,1+σ̃ϕϕ,1)]δϕϕ,1=∠Ẽ1+∠σ̃ϕϕ,1−∠[σ̃ϕϕ,1−ν̃1(σ̃rr,1+σ̃θθ,1)]δrθ,1=∠Ẽ1−∠[1+ν̃1]$
(9)

The symbol $∠$ denotes the phase of the argument. The complex Young's modulus is related to the Lamé parameters by $Ẽ1=μ̃1(3λ̃1+2μ̃1)/(λ̃1+μ̃1)$, and the complex Poisson's ratio by $ν̃1=λ̃1/[2(λ̃1+μ̃1)]$ [30]. It should be noted that the stress components at any given point are, in general, not in phase with one another, which is taken into account in Eqs. (8) and (9).

The heat generation term specified by the relations in Eqs. (7), (8), and (9) represents the extension of loss predictions based on the hysteresis loop for one-dimensional stress-strain states [12] to the three-dimensional state under consideration in this work. Direct application of the relations, along with the nonzero stress components, yields the volumetric heat generation induced in the binder material. Since the particle is treated as rigid, there is no viscoelastic heat generation in the spherical inclusion (i.e., q2 = 0). Thus, the only bulk heating in the system due to material dissipation is the viscoelastic heating of the binder.

### Heat Transfer Equation.

The temperature distribution and evolution in the system are modeled in this work by application of the Fourier law of conduction and are solved by using a finite difference approach. Thermal isotropy is assumed in both the binder and the particle. The strains due to thermal expansion and the corresponding losses due to thermoelastic damping [3134] are neglected in the temperature predictions. This approximation is valid only in the limit of small temperature increases, or, equivalently, over a time interval sufficiently small such that only small temperature rises occur, with the result that the thermal strains are negligible in comparison to those from the applied mechanical excitation. The residual thermal stresses at the particle–binder interface, which develop due to the mismatch of the thermal expansion coefficients and would introduce a nonzero offset into the harmonic stresses, are therefore also not computed here. By the same considerations, the temperature-dependence of the material properties is also neglected. The magnitudes of the stresses in Eq. (4) and the time-averaged heat generation term in Eq. (7) are thus approximated as stationary.

Since there is no variation in the stress distribution in the $ϕ$-direction, and hence, no variation in the heat generation in that direction, heat diffusion takes place only in the -plane. The Fourier law of conduction in spherical coordinates [35] thus simplifies, in each medium m, to
$ρmcpm∂T∂t=km[1r2∂∂r(r2∂T∂r)+1r2 sin θ∂∂θ(sin θ∂T∂θ)]+qm$
(10)

where k and cp denote, respectively, the thermal conductivity and specific heat capacity at constant pressure.

The boundary conditions in the θ-coordinate, in both the particle and the binder, are specified by the requirement of zero heat flux in the θ-direction at θ = 0 and θ = π, due to the symmetry of the stresses and heat generation about the θ = 0 axis. At the particle–binder boundary, r = a, the conductive interface requires continuity of the temperature and heat flux. At the outer boundary of the binder, r = R, the convective surface condition is applied, with convection coefficient U0. Finally, the remaining boundary condition in the radial coordinate in the particle is given by continuity of the heat flux at r = 0. The boundary conditions in the r- and θ-coordinates can thus be summarized, respectively, as
$k2∂T∂r(0,θ,t)=−k2∂T∂r(0,θ+π,t)T(a−,θ,t)=T(a+,θ,t)k2∂T∂r(a−,θ,t)=k1∂T∂r(a+,θ,t)k1∂T∂r(R,θ,t)=U0[T0−T(R,θ,t)]km∂T∂θ(r,0,t)=0km∂T∂θ(r,π,t)=0$
(11)
where the superscripts − and + denote the limits from the negative and positive radial directions, respectively, and T0 is the constant ambient temperature. The initial condition is specified as an isothermal state at the ambient temperature
$T(r,θ,0)=T0$
(12)

The heat transfer equation, Eq. (10), was solved in this work by a finite difference method, implemented in MATLAB®. The Crank-Nicholson method [36], an implicit, second-order scheme in time, was used to solve for the evolution of the temperature. Similarly, for the spatial derivatives, a second-order central difference approximation was used for all interior mesh points, and also for the boundary condition at r = R. For the boundary conditions at the particle–binder interface, r = a, second-order backward and forward difference approximations were used for the negative and positive limits, respectively. Due to the singularities in the analytical heat equation at r = 0, θ = 0, and θ = π, first-order finite difference approximations were employed for the application of the boundary conditions at those points.

Additional effects which are not taken into account in the proposed model, including interface effects (such as debonding and intermittent contact), chemical decomposition or phase changes of the energetic crystal, and, as previously mentioned, thermal strains and the temperature-dependence of the material properties, are expected to play significant roles in the thermomechanics of energetic composites subjected to high-frequency excitation. However, since the thermal model presented here is valid in the limit of small temperature rises and perfect bonding of the crystal and binder, the analytical equations and temperature predictions may provide insight into the initial heating of these composite systems. Moreover, the relations provide a basis for understanding the quantitative effect of the various material and excitation parameters on the initial thermomechanics of such systems.

## Numerical Results and Discussion

In order to illustrate the model predictions for a typical energetic crystal embedded in a viscoelastic medium, the case of a spherical cyclotetramethylene-tetranitramine (HMX) crystal situated in a Sylgard® 184 binder is considered here. Sylgard® is a silicone-based elastomer, developed by the Dow Corning Corporation, Midland, MI and is typically employed as an encapsulant in electrical applications [37,38]. It has, however, also been used as the binder material in polymer-based energetic composites [39,40]. Moreover, for the purpose of the future experimental validation of model predictions, the optical transparency of Sylgard® may prove useful in providing a pathway for measurements on embedded crystals in configurations similar to that considered in this work.

The material properties of the Sylgard® 184 binder were taken as: density ρ1 = 1030 kg/m3 [38], longitudinal wave speed v1 = 1100 m/s [41], shear wave speed $vs1=570 m/s$ [41], longitudinal wave attenuation χ1 = 2.4 dB/MHz/cm [42], thermal conductivity k1 = 0.27 W/(m-K) [38], and thermal diffusivity γ1 = 1.02 × 10–7 m2/s [37]. The specific heat capacity is then specified by the relation $cp1=k1/(ρ1γ1)$. No suitable value could be found in the literature for the shear wave attenuation in Sylgard®. In order to specify the missing material parameter, the assumption was made that the imaginary part of the Poisson's ratio is negligible (i.e., $Im[ν̃1]≈0$), an approximation that is supported by measurements on other polymers [43,44]. All other material moduli can then be computed through the relations for linear viscoelastic media [30].

Since the crystal is modeled as a rigid particle, only the density of HMX is required to characterize its mechanical properties, which was taken as ρ2 = 1910 kg/m3 [45]. The particle radius was set as a = 0.25 mm, and the thermal properties were specified as k2 = 0.4184 W/(m-K) and $cp2=1015 J/(kg‐K)$, taken at 21 °C [39]. For the heat transfer equation, the convection coefficient of the surrounding fluid was set as U0 = 5 W/(m2-K), which is within the range for free convection of air [35], and the constant ambient (and initial) temperature T0 was assumed to be 21 °C for consistency with the property specifications. The outer radius of the binder was specified as 20 times the particle radius: R = 5 mm. Finally, except where the excitation frequency and amplitude are varied in Sec. 4.3, the incident wave frequency is given as f = 500 kHz (f = ω/[2π]) and the amplitude is specified as $Φ′=1 μm$ at the Sylgard® outer boundary z = −R, x = 0 (thus, the wave potential amplitude at z = 0 is $Φ̃0=Φ′eiα̃1R/α̃1$).

### Induced Stresses.

The analytical model for the stresses induced by the harmonic ultrasonic excitation was applied to the HMX–Sylgard® system. The magnitudes of the nonzero stress components in the -plane are presented in Fig. 2 for the radial and in-plane shear stresses, and in Fig. 3 for the polar and azimuthal stresses. Note that only the stresses induced in the surrounding binder medium are computed, and the rigid crystal is shown as a null region in the figures. It should also be noted that an arbitrarily fine grid was used for the generation of these distributions, since the relations are analytical, but this grid differs from the mesh used for the finite difference method to solve the heat transfer equation, which is detailed in Sec. 4.2. As is evident in Fig. 2(a) the radial stress shows a maximum near the front edge of the crystal, at which point the constructive interference due to scattering is maximal, and also a local maximum in stress just behind the crystal, attributable to the harmonic motion of the crystal and binder in the z-direction. Additional stress concentrations of lower magnitudes are observed at greater distances in front of the crystal, which are caused by the constructive interference of the incident wave and the waves scattered from the crystal–binder interface.

Fig. 2
Fig. 2
Close modal
Fig. 3
Fig. 3
Close modal

The maxima of the magnitudes of the other stress components are lower than that of the radial stress. The shear stress in Fig. 2(b) also shows stress concentrations near the crystal interface, but at locations offset from the front and rear of the crystal, at which points the shear wave scattering is strongest. The shear waves then propagate at the scattered angles through the lossy medium. The polar stress in Fig. 3(a) exhibits maximum values at larger angles with respect to the incident wave's propagation direction, since the induced normal stresses in the polar direction are strongest at these angles. Finally, in Fig. 3(b), the variation in the azimuthal stress magnitude, as for the radial stress, shows stress concentrations at several locations of constructive interference in front of the crystal, as well as over a small region immediately behind the crystal.

### Heat Generation and Thermal Response.

The time-averaged volumetric heat generation was computed by application of Eqs. (7), (8), and (9) to the stresses computed in Sec. 4.1. The result is presented as Fig. 4, where the same arbitrarily fine grid as was used for the stresses has been employed to show the analytical solution. The maximum in the heat generation is observed directly in front of the crystal, and the elevated heating induced at several additional points of constructive interference in the z-direction in front of the crystal is also evident. Though the heat generation distribution is largely dominated by the radial stress component, the effects of the in-plane shear stress and polar stress are readily apparent as well, with appreciable heating observed along the shear wave scattering angles and along the vertical direction near the crystal. The effects of the azimuthal stress are less prominent, but also influence the topology of the heating distribution, particularly near the crystal–binder interface.

Fig. 4
Fig. 4
Close modal

The heat transfer equation was solved by the finite difference scheme described in Sec. 3.2. The spatial mesh was specified to consist of 241-points in the r-coordinate from r = 0 to R, and 121-points in the θ-coordinate from θ = 0 to π. The solution was then mirrored over the line of symmetry θ = 0 for the purpose of plotting the full distributions. The condition of temperature continuity at the particle–binder interface was enforced by specifying r = a as one of the mesh points. For the implicit Crank-Nicholson method [36], 100 time steps were used to advance the temperature solution from t = 0 to 0.5 s. The temperature results with respect to both the spatial and temporal mesh sizes were observed to converge to within a tolerance of 1%.

The predictions for the maximum transient temperatures of the crystal and binder are presented in Fig. 5. Though the maximum heat generation is induced near the crystal–binder interface, noticeably larger temperature rises are observed in the binder due to its lower thermal conductivity as compared to the crystal. The location of the maximum crystal temperature is, as expected, at the nearest point on the incident side (x = 0, z = −0.25 mm), and the location of the maximum binder temperature is a small distance along the z-direction in front of the crystal, at x = 0, z ≈ −0.38 mm for this particular excitation frequency. The maximum crystal temperature rises by about 32.6 °C over the 0.5 s of excitation, and the rate of temperature increase for the crystal approaches approximately 55 °C/s.

Fig. 5
Fig. 5
Close modal

The computed temperature distribution at t = 0.5 s is shown as Fig. 6. The point of the maximum thermal response in the binder, located close to the crystal interface, is clearly evident. The local maxima in temperature that appear at the stress concentrations along the z-direction in front of the crystal are also apparent, as are the lesser temperature rises at angles offset from the incident wave propagation direction, where the shear stress and polar stress contribute significantly to the heat generation. The minimum temperature in the system is located at a point within the back half of the crystal, where the thermal diffusion induced by the greater heat generation near the front surface of the crystal meets the diffusion induced by the smaller heat generation near the back surface. The results presented here indicate that viscoelastic heating of the binder material induced by the applied excitation is significant, particularly near the crystal, and that this heating mechanism is likely to play an important role in the formation of hot spots at the crystal–binder interface. It should be further noted that this heating is predicted in the absence of voids or delamination between the crystal and binder, but that these interface effects may also contribute substantially to the response through frictional heating.

Fig. 6
Fig. 6
Close modal

### Effects of Excitation Amplitude and Frequency.

In the context of the cyclic loading of energetic composites, the excitation parameters, in contrast with the material properties in the system, are considered tunable to a degree, and so the quantitative effect of these parameters on the thermal response is of interest here. Specifically, the effects of varying the incident wave amplitude and frequency on the maximum crystal temperature are investigated.

Figure 7(a) shows the effect of varying the excitation amplitude on the maximum crystal temperature at t = 0.5 s and the corresponding rate of temperature increase. The amplitude $Φ′$ is, as before, specified at the Sylgard® outer boundary z = −5 mm, x = 0. Since the dynamic model for the stresses is linear, each of the stress components scales directly with the amplitude and, therefore, the volumetric heat generation scales with the square of the amplitude. As such, the temperature increases and the corresponding heating rates exhibit a simple quadratic variation with the incident wave amplitude, which is evident in Fig. 7(a). Though the predicted thermal response shows large temperature increases for higher excitation amplitudes, the resulting thermal strains, changes in material properties, interface effects, and physical and chemical changes would have to be taken into account to accurately assess the temperature evolution and distribution.

Fig. 7
Fig. 7
Close modal

The effect of varying the excitation frequency is presented in Fig. 7(b), again considering both the maximum crystal temperature at t = 0.5 s and the corresponding rate of temperature increase. The variation with frequency also resembles a quadratic dependence, but the relation is more complicated due to the frequency-dependence of the Bessel and Hankel function terms in the analytical solution for the stresses. As a result, the phase differences given in Eq. (9) for the volumetric heat generation term also vary with frequency, in addition to the linear dependence of the heating shown in Eq. (7) and the variation of the stress magnitudes. It should be noted, too, that for the case of linear variation of the wave attenuation coefficients with frequency, specified in this work by the longitudinal coefficient χ1 = 2.4 dB/MHz/cm and the condition $Im[ν̃1]=0$, the phases of the material moduli are independent of frequency. But if the frequency-dependence of the attenuation coefficients is not linear, then the changes in the phases of the material moduli would further contribute to the variation of the phase differences in Eq. (9). As for the amplitude variation, additional effects related to the temperature increases, as well as to the frequency levels, should be considered to accurately assess the thermal response at the higher excitation frequencies.

In addition to the effect of the excitation frequency on the thermal response of the crystal, the effect on the translational motion of the rigid inclusion was also investigated. Figure 8 shows the displacement amplitude of the crystal, which is undergoing harmonic motion in the z-direction, as a function of the excitation frequency, with the incident wave amplitude again set as $Φ′=1 μm$ at the Sylgard® outer boundary z = −5 mm, x = 0. At low excitation frequencies, the particle is observed to translate with an amplitude that is slightly less than the prescribed displacement of the incident compressional wave due to the decrease in amplitude associated with the viscoelastic losses of the binder material, and this difference approaches zero in the zero frequency limit (since the attenuation in the binder is assumed to scale with frequency). Above this region, the displacement amplitude of the crystal increases with frequency to a maximum value near 276 kHz, and then decreases rapidly with further increasing frequencies. The characteristics of this displacement curve are similar to those of the curves presented by Oien and Pao [5].

Fig. 8
Fig. 8
Close modal

## Conclusions

A thermomechanical model of a rigid moveable sphere embedded within an infinite viscoelastic material undergoing mechanical excitation via a harmonic compressional plane wave has been developed to predict the stress, heat generation, and temperature distribution surrounding the inclusion. This model was then examined for the case of a spherical energetic particle of HMX embedded within a binder material of Sylgard® 184 at an excitation frequency of 500 kHz and an excitation amplitude of 1 μm. The results for this case yield predictions of significant heat generation near the incident-side of the particle corresponding to the locations of major stress concentrations. The viscoelastic heating effect results in a rate of temperature increase of approximately 55 °C/s for the energetic crystal. These rates can be compared to those estimated for the particle surface by Mares et al. [26]. For a 750–950 μm HMX crystal excited at 215 kHz, the temperature increase was estimated to be around 37 °C/s, and for a 400–500 μm ammonium perchlorate crystal, it was estimated as high as 125 °C/s. The analytical results agree well with these in order of magnitude although differences in frequency, particle size (and morphology), as well as amplitudes of excitation make direct quantitative comparison difficult.

This model was presented in order to explore the nature of the heat generation due to stresses created by the interference of the scattered and incident waves near a rigid inclusion and, specifically, to investigate how this mechanism of heat generation can be utilized within the context of the ultrasonic excitation of energetic systems. This phenomenon may help to explain how localized hot spots initially develop within crystal–binder systems undergoing ultrasonic excitation, even in the absence of initial flaws or voids. It should, however, be further emphasized here that the scope of the physics addressed by this model is expected to be adequate in characterizing the heating in the system only for short excitation times. Once certain temperature thresholds are reached, additional mechanisms associated with heat generation such as chemical decomposition or phase changes of the energetic crystal, significant thermal stresses at the crystal–binder interface, changes in material properties, and debonding between the crystal and binder may become substantial. As such, the phenomena associated with these additional heating mechanisms lie beyond the scope and validity of the proposed model. A qualitative diagram representing the temperature and temporal regions in which the proposed model is expected to remain valid is presented in Fig. 9. The temperatures and times beyond which those additional mechanisms are expected to become significant are shown conceptually in the diagram as well.

Fig. 9
Fig. 9
Close modal

It is also important to note that the particle morphology of the embedded inclusion has been experimentally shown to greatly affect the stress distributions as well as the heating rates associated with the ultrasonic excitation of similar energetic systems [29]. Since the model with a spherical inclusion does not account for those stress concentrations associated with irregular particle morphologies, the predictions shown in the present work may represent a conservative estimate of the heating magnitudes in these systems. By adapting this model to account for such particle-specific effects, a more accurate representation of the underlying physics may be captured.

Future efforts should attempt to address the time scales at which thermal stresses become significant, and at which debonding occurs, as well as the effects of nonlinear binder behavior and temperature dependent properties. Moreover, the heating mechanisms associated with frictional effects at the crystal–binder interface (intermittent contact, rubbing, etc.), which may dominate the thermal behavior of the system after debonding, should be investigated. Additionally, experimental validation of the predicted harmonic particle displacement and temperature distribution in an idealized particle–binder system should be demonstrated through a wide range of applicable excitation parameters utilizing various materials. As the stress solution is analytical, it affords itself well to the alteration of parameters including excitation frequency without incurring lengthy simulation times due to the small time steps required to resolve periodic stresses at high frequency by the finite element method. This characteristic may also prove useful in determining viscoelastic material properties of binder materials under various excitation frequencies and amplitudes.

## Acknowledgment

This research is supported by the U.S. Office of Naval Research through ONR Grant Nos. N00014-10-1-0958 and N00014-11-1-0466, by the U.S. Air Force Office of Scientific Research through Award No. FA9550-15-1-0102, and by the National Science Foundation Graduate Research Fellowship Program through Grant No. DGE-1333468. A preliminary version of this work appeared in the proceedings of the ASME 2016 International Mechanical Engineering Congress & Exposition [46].

## References

1.
Ying
,
C. F.
, and
Truell
,
R.
,
1956
, “
Scattering of a Plane Longitudinal Wave by a Spherical Obstacle in an Isotropically Elastic Solid
,”
J. Appl. Phys.
,
27
(
9
), pp.
1086
1097
.
2.
Pao
,
Y. H.
, and
Mow
,
C. C.
,
1963
, “
Scattering of Plane Compressional Waves by a Spherical Obstacle
,”
J. Appl. Phys.
,
34
(
3
), pp.
493
499
.
3.
Thau
,
S. A.
,
1967
, “
Radiation and Scattering From a Rigid Inclusion in an Elastic Medium
,”
ASME J. Appl. Mech.
,
34
(
2
), pp.
509
511
.
4.
Mow
,
C. C.
, and
Pao
,
Y. H.
,
1971
, “
The Diffraction of Elastic Waves and Dynamic Stress Concentrations
,” United States Air Force Project, RAND Report No. R-482-PR.
5.
Oien
,
M. A.
, and
Pao
,
Y. H.
,
1973
, “
Scattering of Compressional Waves by a Rigid Spheroidal Inclusion
,”
ASME J. Appl. Mech.
,
40
(
4
), pp.
1073
1077
.
6.
Gaunaurd
,
G. C.
, and
Uberall
,
H.
,
1978
, “
Theory of Resonant Scattering From Spherical Cavities in Elastic and Viscoelastic Media
,”
J. Acoust. Soc. Am.
,
63
(
6
), pp.
1699
1712
.
7.
Beltzer
,
A.
,
Robinson
,
B.
, and
Rudy
,
N.
,
1979
, “
The Effect of Random Compressional Waves on a Rigid Sphere Embedded in an Elastic Medium
,”
J. Sound Vib.
,
66
(
4
), pp.
513
519
.
8.
Beltzer
,
A. I.
,
1980
, “
Random Response of a Rigid Sphere Embedded in a Viscoelastic Medium and Related Problems
,”
ASME J. Appl. Mech.
,
47
(
3
), pp.
499
503
.
9.
Sessarego
,
J. P.
,
Sageloli
,
J.
,
Guillermin
,
R.
, and
Uberall
,
H.
,
1998
, “
Scattering by an Elastic Sphere Embedded in an Elastic Isotropic Medium
,”
J. Acoust. Soc. Am.
,
104
(
5
), pp.
2836
2844
.
10.
Avila-Carrera
,
R.
, and
Sanchez-Sesma
,
F. J.
,
2006
, “
Scattering and Diffraction of Elastic P-and S-Waves by a Spherical Obstacle: A Review of the Classical Solution
,”
Geofís. Int.
,
45
(
1
), pp.
3
21
.
11.
Tauchert
,
T. R.
,
1967
, “
Heat Generation in a Viscoelastic Solid
,”
Acta Mech.
,
3
(
4
), pp.
385
396
.
12.
Brinson
,
H. F.
, and
Brinson
,
L. C.
,
2008
,
Polymer Engineering Science and Viscoelasticity: An Introduction
,
Springer
, New York.
13.
Dinzart
,
F.
,
Molinari
,
A.
, and
Herbach
,
R.
,
2008
, “
Thermomechanical Response of a Viscoelastic Beam Under Cyclic Bending: Self-Heating and Thermal Failure
,”
Int. Appl. Mech.
,
60
(
1
), pp.
59
85
.
14.
Dimarogonas
,
A. D.
, and
Syrimbeis
,
N. B.
,
1992
, “
Thermal Signatures of Vibrating Rectangular Plates
,”
J. Sound Vib.
,
157
(
3
), pp.
467
476
.
15.
Chervinko
,
O. P.
, and
Senchenkov
,
I. K.
,
2002
, “
The Coupled Thermomechanical State of a Notched Viscoelastic Rectangular Plate Under Harmonic Loading
,”
Int. Appl. Mech.
,
38
(
2
), pp.
209
216
.
16.
Chervinko
,
O. P.
,
2004
, “
Calculating the Critical Parameters Characterizing the Thermal Instability of a Viscoelastic Prism With a Stress Concentrator Under Harmonic Compression
,”
Int. Appl. Mech.
,
40
(
8
), pp.
916
922
.
17.
Chervinko
,
O. P.
,
Senchenkov
,
I. K.
, and
Yakimenko
,
N. N.
,
2007
, “
Vibrations and Self-Heating of a Viscoelastic Prism With a Cylindrical Inclusion
,”
Int. Appl. Mech.
,
43
(
6
), pp.
647
653
.
18.
Katunin
,
A.
, and
Fidali
,
M.
,
2012
, “
Self-Heating of Polymeric Laminated Composite Plates Under the Resonant Vibrations: Theoretical and Experimental Study
,”
Polym. Compos.
,
33
(
1
), pp.
138
146
.
19.
Woods
,
D. C.
,
Miller
,
J. K.
, and
,
J. F.
,
2015
, “
On the Thermomechanical Response of HTPB-Based Composite Beams Under Near-Resonant Excitation
,”
ASME J. Vib. Acoust.
,
137
(
5
), p.
054502
.
20.
Miller
,
J. K.
,
Woods
,
D. C.
, and
,
J. F.
,
2014
, “
Thermal and Mechanical Response of Particulate Composite Plates Under Inertial Excitation
,”
J. Appl. Phys.
,
116
(
24
), p.
244904
.
21.
Hinders
,
M. K.
,
Fang
,
T. M.
,
McNaughton Collins
,
M. F.
, and
Collins
,
J. J.
,
1994
, “
An Analytic Solution for Energy Deposition in Model Spherical Tumors Undergoing Ultrasound-Hyperthermia Treatments
,”
Phys. Med. Biol.
,
39
(
1
), pp.
107
132
.
22.
Bowden
,
F. P.
, and
Yoffe
,
A. D.
,
1952
,
Initiation and Growth of Explosion in Liquids and Solids
,
Cambridge University Press
, New York.
23.
,
C. L.
,
1963
, “
Shock and Hot Spot Initiation of Homogeneous Explosives
,”
Phys. Fluids
,
6
(
3
), pp.
375
381
.
24.
Chaudhri
,
M. M.
, and
Field
,
J. E.
,
1974
, “
The Role of Rapidly Compressed Gas Pockets in the Initiation of Condensed Explosives
,”
Proc. R. Soc London A
,
340
(
1620
), pp.
113
128
.
25.
Swallowe
,
G. M.
, and
Field
,
J. E.
,
1982
, “
The Ignition of a Thin Layer of Explosive by Impact; the Effect of Polymer Particles
,”
Proc. R. Soc. London A
,
379
(
1777
), pp.
389
408
.
26.
Mares
,
J. O.
,
Miller
,
J. K.
,
Gunduz
,
I. E.
,
,
J. F.
, and
Son
,
S. F.
,
2014
, “
Heat Generation in an Elastic Binder System With Embedded Discrete Energetic Particles Due to High-Frequency, Periodic Mechanical Excitation
,”
J. Appl. Phys.
,
116
(
20
), p.
204902
.
27.
Chen
,
M. W.
,
You
,
S.
,
Suslick
,
K. S.
, and
Dlott
,
D. D.
,
2014
, “
Hot Spots in Energetic Materials Generated by Infrared and Ultrasound, Detected by Thermal Imaging Microscopy
,”
Rev. Sci. Instrum.
,
85
(
2
), p.
023705
.
28.
You
,
S.
,
Chen
,
M. W.
,
Dlott
,
D. D.
, and
Suslick
,
K. S.
,
2015
, “
Ultrasonic Hammer Produces Hot Spots in Solids
,”
Nat. Commun.
,
6
(
4
), p.
6581
.
29.
Miller
,
J. K.
,
Mares
,
J. O.
,
Gunduz
,
I. E.
,
Son
,
S. F.
, and
,
J. F.
,
2016
, “
The Impact of Crystal Morphology on the Thermal Responses of Ultrasonically-Excited Energetic Materials
,”
J. Appl. Phys.
,
119
(
2
), p.
024903
.
30.
Borcherdt
,
R. D.
,
2009
,
Viscoelastic Waves in Layered Media
,
Cambridge University Press
, New York.
31.
Zener
,
C.
,
1938
, “
Internal Friction in Solids II: General Theory of Thermoelastic Internal Friction
,”
Phys. Rev.
,
53
(
1
), pp.
90
99
.
32.
Ting
,
E. C.
,
1973
, “
Dissipation Function of a Viscoelastic Material With Temperature-Dependent Properties
,”
J. Appl. Phys.
,
44
(
11
), pp.
4956
4960
.
33.
Bishop
,
J. E.
, and
Kinra
,
V. K.
,
1996
, “
Equivalence of the Mechanical and Entropic Descriptions of Elastothermodynamic Damping in Composite Materials
,”
Mech. Compos. Mater. Struct.
,
3
(
2
), pp.
83
95
.
34.
Lakes
,
R.
,
1997
, “
Thermoelastic Damping in Materials With a Complex Coefficient of Thermal Expansion
,”
J. Mech. Behav. Mater.
,
8
(
3
), pp.
201
216
.
35.
Incropera
,
F. P.
,
DeWitt
,
D. P.
,
Bergman
,
T. L.
, and
Lavine
,
A. S.
,
2007
,
Introduction to Heat Transfer
,
Wiley
, Hoboken, NJ.
36.
,
B.
,
2006
,
A Friendly Introduction to Numerical Analysis
,
Pearson Prentice Hall
37.
Dow Corning Corporation
,
1986
, “
Information About High Technology Materials Sylgard® 182 & 184 Silicone Elastomers
,” Data Sheet No. 61-113C-01.
38.
Dow Corning Corporation
,
2014
, “
Product Information, Sylgard® 184 Silicone Elastomer
,” Data Sheet No. 11-3184B-01.
39.
Gibbs
,
T. R.
, and
Popolato
,
A.
, eds.,
1980
,
LASL Explosive Property Data
,
University of California Press
, Berkeley, CA.
40.
Shepodd
,
T.
,
Behrens
,
R.
,
Anex
,
D.
,
Miller
,
D.
, and
Anderson
,
K.
,
1997
, “
Degradation Chemistry of PETN and Its Homologues
,” Sandia National Laboratories, Report No. SAND–97-8684C.
41.
Millett
,
J. C. F.
,
Whiteman
,
G.
,
Stirk
,
S. M.
, and
Bourne
,
N. K.
,
2011
, “
Shear Strength Measurements in a Shock Loaded Commercial Silastomer
,”
J. Appl. Phys.
,
44
(
18
), p.
185403
.
42.
Garvin
,
K. A.
,
Hocking
,
D. C.
, and
Dalecki
,
D.
,
2010
, “
Controlling the Spatial Organization of Cells and Extracellular Matrix Proteins in Engineered Tissues Using Ultrasound Standing Wave Fields
,”
Ultrasound Med. Biol.
,
36
(
11
), pp.
1919
1932
.
43.
Pritz
,
T.
,
1998
, “
Frequency Dependences of Complex Moduli and Complex Poisson's Ratio of Real Solid Materials
,”
J. Sound Vib.
,
214
(
1
), pp.
83
104
.
44.
Pritz
,
T.
,
2000
, “
Measurement Methods of Complex Poisson's Ratio of Viscoelastic Materials
,”
Appl. Acoust.
,
60
(
3
), pp.
279
292
.
45.
Agrawal
,
J. P.
,
2010
,
High Energy Materials: Propellants, Explosives, and Pyrotechnics
,
Wiley
, Weinheim, Germany.
46.
Mares
,
J. O.
,
Woods
,
D. C.
,
Baker
,
C. E.
,
Son
,
S. F.
,
,
J. F.
,
Bolton
,
J. S.
, and
Gonzalez
,
M.
,
2016
, “
Localized Heating Due to Stress Concentrations Induced in a Lossy Elastic Medium via the Scattering of Compressional Waves by a Rigid Spherical Inclusion
,”
ASME
Paper No. IMECE2016-68219.