Localized Heating Near a Rigid Spherical Inclusion in a Viscoelastic Binder Material Under Compressional Plane Wave Excitation

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 [1–10]. 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 [11–21]. 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 [22–25]. 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 [26–29]. 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.

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 [1–10]. 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 ρ₁, 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 ρ₂ 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

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 $ũr$ (i.e., the radial component) and the displacement component in the in-plane rotational direction $ũθ$⁠. 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.

where the h_n are the spherical Hankel functions of the first kind, and the $Ã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 $Ã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.

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.

where the ε_** are the strain components and t₀ 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.

The symbol $∠$ denotes the phase of the argument. The complex Young's modulus is related to the Lamé parameters by $Ẽ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., q₂ = 0). Thus, the only bulk heating in the system due to material dissipation is the viscoelastic heating of the binder.

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 [31–34] 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.

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

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.

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 ρ₁ = 1030 kg/m³ [38], longitudinal wave speed v₁ = 1100 m/s [41], shear wave speed $vs1=570 m/s$ [41], longitudinal wave attenuation χ₁ = 2.4 dB/MHz/cm [42], thermal conductivity k₁ = 0.27 W/(m-K) [38], and thermal diffusivity γ₁ = 1.02 × 10^–7 m²/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 ρ₂ = 1910 kg/m³ [45]. The particle radius was set as a = 0.25 mm, and the thermal properties were specified as k₂ = 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 U₀ = 5 W/(m²-K), which is within the range for free convection of air [35], and the constant ambient (and initial) temperature T₀ 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$⁠).

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 rθ-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.

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.

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.

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.

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.

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.

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 χ₁ = 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].

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.

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.

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].