A Finite-Deformation Constitutive Model of Particle-Binder Composites Incorporating Yield-Surface-Free Plasticity

Particle-binder composite materials consist of a large concentration of hard particles, called fillers, randomly dispersed in the matrix of a soft material. Generally, fillers are used to enhance the mechanical properties of the soft material. For example, filled elastomers such as carbon-black and silica filled rubbers [1,2] have been shown to possess superior stiffness, strength, and damping properties over natural rubber, making them suitable for application in automotive parts such as tires and bearing seals and in structures providing vibration and shock isolation to mechanical systems.

Another class of particle-binder composites, called energetic composite materials or plastic-bonded explosives (PBX), consists of explosive crystals and, in some formulations, metal fuel powder, embedded in a binder composed mainly of a soft polymer and a plasticizer. Since their initial development at Los Alamos National Laboratory in 1947, PBXs have been commonly used as ammunition in defense applications. Examples of explosives used include cyclotrimethylenetrinitramine (RDX) and cyclotetra-minetetranitramine (HMX), while the binder composition includes polymers like hydroxyl-terminated polybutadiene (HTPB) and Estane, and plasticizers like dioctyl adipate (DOA) and isodecyl pelargonate (IDP), along with smaller concentrations of various additives [3]. In a PBX binder, the polymer, apart from providing structural integrity to the explosive, reduces their impact and vibrational sensitivity [4–6], while the plasticizer improves their mechanical properties and processability [7]. The metal fuel, usually aluminum, is used to enhance the blast effects. Therefore, aluminized PBXs are typically used in naval weapons and missile warheads [8].

Energetic composites are commonly used in applications requiring high mechanical confinement or compression [9,10]. Additionally, these materials may be subjected to diverse loading conditions, ranging from low to high strain-rate cyclic, vibrational and impact loading during transport, storage and handling over their operational life. Since PBXs are designed to detonate under specific external energy stimulus, such loading conditions may alter their mechanical response, rendering them unpredictable and unsafe [11]. Therefore, understanding and predicting the mechanical response of PBX under different loading conditions is of particular interest to the defense and propulsion communities.

Earlier experimental studies on the mechanical behavior of PBX under uniaxial load have shown a dependence on strain rate and temperature [12–14]. Several constitutive models in the context of small and large deformation mechanics have been developed to model their rate and temperature-dependent behavior. For instance, Bardenhagen et al. [15] proposed a large deformation viscoelastic model for PBX binder materials using a Blatz-Ko hyperelasticity formulation and a Maxwell element. Clements and Mas [16] proposed a viscoelastic constitutive theory for filled polymers using the Boltzmann superposition principle, where the filler particles were modeled as randomly positioned elastic ellipsoidal particles. Composite stress relaxation functions were developed with their prony series coefficients dependent on filler concentration, ellipsoidal aspect ratio, and polymer moduli. The viscoelastic cracking constitutive model developed by Bennet et al. [17] combined a linear viscoelastic Maxwell’s model with the isotropic damage model by Addessio and Johnson [18]. It was based on the mechanism of microcracking and derived from the statistical crack mechanics approach proposed by Dienes [19] for brittle materials. The model was developed to predict the nonlinear stress–strain response of the material, as well as the strain-softening and nonlinearity observed due to extensive cracking at larger deformations. The model has since been extensively implemented in finite element analyses and used to predict the thermo-mechanical behavior of PBX and study hotspot generation [6,20]. In the context of PBX undergoing low-frequency vibrational loading, nonlinear viscoelastic models for a mass-material system undergoing base excitation were proposed by Paripovic and Davies [21,22].

While the majority of constitutive formulations proposed for the PBX model the material as viscoelastic due to the rubbery nature of the polymer binder, it has been shown that these materials may exhibit considerable plastic deformation with increasing confinement. Uniaxial confined compression tests carried out on inert mock sugar-based specimens for PBX 9501 [9,10] under different confining pressures revealed increasing plastic deformation and strain hardening with increasing pressure. Recently, Agarwal and Gonzalez [23] conducted unconfined compression tests at room temperature and low strain rate (10⁻³/s) on cylindrical specimens of three mock explosive formulations (Fig. 1) based on PBXN-109 [24], with sugar as a substitute for the explosive HMX crystals. The three formulations contained 85% w/w of solids (sugar and aluminum) but differed from each other by the amount of aluminum content in the solids composition, namely, 85–00 (0% Al in solids, 0% in total), 85–15 (15% Al in solids, 12.75% in total) and 85–30 (30% additive in solids, 25.5% in total). The monotonic nominal stress–strain response of the 85–00 specimen (Fig. 2(a)) exhibited a quasi-brittle behavior [25], with an initial nonlinear increase in stress until a peak stress level at around 10–11% strain, followed by strain-softening due to the evolution of microcracks into larger transgranular fracture leading to a loss of strength. Such behavior has been recorded and studied extensively for many non-aluminized PBX formulations [6,12,26]. However, the aluminized specimens, i.e., 85–15 (Fig. 2(b)) and 85–30 (Fig. 2(c)) exhibited a more ductile plastic flow behavior with strain hardening, indicating that the presence of aluminum could cause the material to deform plastically even in the absence of confinement.

Additionally, plastic deformation and damage in PBX have been shown to be associated with hotspot mechanisms [27,28], thus directly impacting their ignition characteristics. Therefore, the development of a constitutive model capable of characterizing the plastic behavior of explosives is extremely important for proper modeling and prediction of their mechanical behavior. Constitutive models for PBX incorporating plasticity have been proposed for cyclic deformation [29] and mild impact [30]. However, the models have been developed using small deformation theory and, therefore, are inapplicable to the large deformation mechanical behavior of explosives.

Lastly, Agarwal and Gonzalez [23] also conducted cyclic compression tests of the mock PBX specimens at room temperature and low strain rate (10⁻³/s). Figure 3 shows the nominal cyclic stress–strain response of the specimens. The strain-controlled tests were carried out by initially loading the specimen to a maximum strain level, and thereafter partially unloading–reloading between the maximum and a minimum strain. The optimal minimum-maximum strain levels were identified for each mock PBX composition through a systematic procedure detailed in Ref. [23] that was developed to ensure sufficiently large elastoplastic deformation while preventing the material to reach its ultimate strength as well as keep the specimen in compression throughout the loading process. Through this procedure, the cyclic strain levels were identified as 5–9% for 85–00, 8–16% for 85–15, and 6–9% for 85–30. Observable response attributes (indicated in the figures) include the following: (1) a nonlinear, continuous elastoplastic response without apparent yield, (2) hysteresis, and (3) cyclic stress-softening with eventual stabilization. It is also worth noting that the curvature of the initial loading path is different from those of subsequent reloading paths. The former is convex and the latter are concave, indicating irreversible changes in the material behavior during the initial loading itself. This behavior has also been previously observed for cyclic compression of aluminized RDX-based PBX in HTPB binder [31].

It is apparent that the observed mechanical response of mock PBX specimens exhibits similar characteristics to filled elastomers at large deformations. Constitutive models of filled elastomers have been proposed by adopting both phenomenological [32–37] and micromechanical [38–40] frameworks. For instance, Ayoub et al. [33] proposed a Zener-type visco-hyperelastic constitutive model of rubber-like materials that accounts for the Mullins effect [41,42], continuous stress-softening, and permanent residual strains by utilizing the network alteration theory [43,44]. Raghunath et al. [39] extended the micromechanical dynamic flocculation model developed by Klüppel [45] to include time-dependent effects typical of filled elastomers. It is worth mentioning that the majority of these models attribute the observed hysteresis and other inelastic phenomena to time-dependent viscous overstress in the rubber matrix. In contrast, only a limited number of models consider these phenomena as time-independent plastic deformation mechanisms. A case in point is the phenomenological elastoplastic constitutive model proposed by Kaliske and Rothert [46] that captures the rate-independent process of internal material friction due to irreversible polymer slippage on the filler surface and plastic deformation of the filler particles. This multi-yield-surface rheological model is comprised of a finite number of elastoplastic Prandtl elements arranged in parallel with an elastic member. More recently, yield-surface-free endochronic plasticity [47] was utilized by Netzker et al. [48] to model the smooth and hysteretic cyclic stress–strain response of carbon-black filled rubbers, with higher computational efficiency and fewer parameters.

In this paper, we present an elastoplastic constitutive formulation with isotropic damage capable of modeling the cyclic response of mock energetic composite materials. Specifically, a Lagrangian finite-deformation formulation based on the additive decomposition of strain energy [49–53] into elastic and plastic parts is adopted. The formulation uses Ogden’s hyperelastic model [54,55] to predict the rubber-like nonlinear elastic response of the polymeric binder and a hereditary (path-dependent) yield-surface-free endochronic plasticity theory [47], based on the concept of internal state variables [56], to account for irreversible deformations. A discontinuous isotropic damage model [57,58] is utilized to model the stress-softening that occurs during unloading, and an endochronic material scale function is utilized to model progressive cyclic stress-softening and attainment of stable cyclic response. It is worth mentioning that while the material scale function has been utilized in previous works to describe the cyclic hardening or softening behavior of metals and alloys [47,59–61], its application to filled polymeric materials with complex microstructure like particle-binder composites has not been demonstrated in the open literature, and thus, it is one of the primary contributions of this paper. The number of model parameters is a function of the number of active Ogden terms and endochronic branches; therefore, a significant number of parameters may need to be identified. Accordingly, we develop a parameter identification method based on a nonlinear multivariate least-squares problem, which is expected to be non-convex and affected by multiple local optima. The range of behavior predicted by the proposed model is demonstrated by calibrating the model parameters for the three mock PBX material compositions [23] under both monotonic and cyclic compressive loading. Finally, the versatility and the capability of the model to represent the material, specifically, the material composition, is demonstrated by identifying a correlation between the calibrated model parameters and the unique stress–strain response of each mock PBX specimen.

The paper is organized as follows. The constitutive model is presented in Sec. 2, followed by Sec. 3 which provides a discrete numerical procedure to solve for stresses along a loading path. Section 4 presents the model parameter identification procedure and the calibrated material properties of the three mock PBX formulations. Section 5 shows a detailed comparative analysis of the three specimen formulations and presents the correlation identified between the estimated mechanical properties and the mechanical behavior, as well as the composition of the specimens. Finally, a summary and concluding remarks are presented in Sec. 6.

Additionally, filled rubber or soft polymeric materials in general exhibit a stress-softening behavior when unloaded from a maximum previous deformation. This phenomenon, widely known as the Mullins effect [41,42], has been discussed extensively in the literature. The micromechanics behind the occurrence of this phenomenon is highly debated, and numerous theories have been suggested, such as rupture of filler-polymer bonds [69,70], molecular slipping on the filler surface [71,72], and breakage of filler clusters [73,74]. However, most phenomenological descriptions of this effect employ the thermodynamic framework of continuum damage mechanics (CDM) [57,58], where the loss of the material’s ability to sustain load due to damage is described by applying a scalar reduction factor, also known as a damage variable, to the stress in a hypothetical, undamaged material. Various models based on this principle have been proposed, generally classified according to the evolution of damage being discontinuous [49,75–77] or continuous [78–81] (for more information on the physical interpretation and classification of phenomenological models for Mullins effect [82]). While discontinuous damage models predict the same reloading response as the unloading response during cyclic loading (referred to as idealized Mullins effect), continuous models have the capability to model diverging unloading–reloading responses. Since the different unloading–reloading responses as seen in energetic composites can be effectively modeled by the proposed path-dependent constitutive formulation, the damage evolution can be defined as discontinuous and occurring only during unloading for simplicity and computational efficiency.

Therefore, during the initial loading of the virgin material when the current value of the deformation remains equal to its maximum value throughout the material’s history, the damage function does not evolve, and the material’s strain energy remains unchanged. However, during unloading–reloading from a previous maximum deformation, the current deformation measure is lower than the maximum value, and the damage function evolves monotonically between 0 and 1. This in turn reduces the effective strain energy in the material and produces the softening effect.

The damage function is then a monotonic function of the deformation measure in the interval $ϕ∈[3,ϕmax]$⁠, with positive damage parameters m and p.

In the above improved scale function, the additional parameter z_ref, defined as the value of z at which the last reversal of the strain path occurred, adds the capability to model “fading memory” effects exhibited by materials with memory [99]. In turn, the saturation value at steady-state depends on the reference time scale z_ref. It is equal to c for z_ref = 0, i.e., before the first reversal, and approaches c/s for z_ref → ∞, i.e., at cyclic stabilization. It is worth noting that for c > s > 1, the initial saturation value c is greater than the saturation value c/s at cyclic stabilization, i.e., the scale function describes the progressive cyclic softening typically observed in particle-binder composites. Figure 4 shows a schematic representation of the scale function proposed by Lin et al. [61]. It is evident from the figure that the rate of evolution of the material scale function is controlled by β_c, whereas the rate of evolution of the saturation value is controlled by β_s.

The value of intrinsic time scale increment, Δz, is obtained by using a midpoint rule to discretize Eq. (29) together with Eq. (31), i.e., from

The previous nonlinear equation is solved numerically using a Newton-Raphson method.

The range of behavior predicted by the proposed constitutive model is demonstrated by calibrating the model parameters for the three mock PBX formulations (85–00, 85–15, and 85–30) under monotonic and cyclic compression [23]. Therefore, Sec. 4.1 presents the incremental procedure for uniaxial loading under unconfined lateral conditions and incompressible material assumptions. Section 4.2 presents a parameter identification procedure based on a nonlinear multivariate minimization problem and the least-squares principle. Finally, Sec. 4.3 shows the results of the parameter identification.

Following the incremental procedure presented in Sec. 3, we assume that the state of the material, $Cn$⁠, $S~j,ni$⁠, z_n, z_ref,n, and $ϕnmax$⁠, is known at loading step n, and that $Cn+1$ at loading step n + 1 is given. The maximum deformation level $ϕn+1max$ is determined from Eqs. (18) and (19) using

We have developed a parameter identification method based on a nonlinear multivariate minimization problem. The proposed constitutive model has 2M + 2P × N + N + 6 material parameters, namely

1. 2M parameters in the elastic branch: μ_k and α_k (k = 1, …, M).

2. 2P × N elastic parameters in the yield-surface-free endochronic branches: μ_ij and α_ij (i = 1, …, P; j = 1, …, N).

3. N kernel parameters in the yield-surface-free endochronic branches: γ_j(j = 1, …, N).

4. Four material scale function parameters: c, s, β_c, and β_s.

5. Two isotropic damage parameters: m and p.

An iterative procedure is followed to determine the number of terms M and P in the Ogden strain energy series functions for the elastic branch and the endochronic branches, respectively, as well as to determine the number of endochronic branches N. The procedure starts from an initial value of M = P = N = 1 and calculates the least-squares error given by Eq. (45). During each iteration, three error values are computed by independently incrementing the values of M, P, and N by 1. The increment corresponding to the minimum of the three errors is adopted, only if the minimum error is reduced by a specified percentage tolerance (taken as 5% in this study) as compared to the minimum error computed in the previous iteration. Otherwise, the values of M, P, and N obtained during the previous iteration are adopted and the procedure is concluded. This procedure aims to improve the estimation of the model without overly increasing its complexity.

The constitutive model and parameter identification method presented earlier were used to calibrate the experimental monotonic (Fig. 2) and cyclic (Fig. 3) compressive response of 85–00, 85–15, and 85–30 mock PBX cylindrical specimens. The nonlinear multivariate minimization problem was solved in matlab^® [100] using the constrained optimization function fmincon with the default interior-point algorithm. The model-calibrated nominal stress–strain curves are presented in Figs. 5 (monotonic) and 6 (cyclic), and the estimated material parameters for each specimen formulation are listed in Table 1. Comparison of experimental and model curves suggests an excellent agreement between the two for all three formulations. The strain offset was calculated as 2.25% for 85–00, 1.03% for 85–15, and 2.51% for 85–30, and it is denoted in the figures by a dashed line.

For detailed comparative analysis, experimental and model-estimated values of the cumulative energy dissipation (Fig. 7), peak stress (Fig. 8), and apparent modulus (Fig. 9) are compared for cyclic compression of the specimens. The apparent modulus [23,37] is calculated as the slope of the line connecting peak and valley stresses in each cycle, while the cyclic energy dissipation [23] is calculated as the difference between the energy supplied to the material (i.e., the area under the loading path) and the energy recovered after a cycle (i.e., the area under the unloading path). Again, an excellent agreement between the experimental and model values is obtained for the three specimens, demonstrating the capability of the proposed model to describe the strength, stiffness, and softening characteristics of particle-binder composites.

From the observation of calibrated properties of mock PBX specimens in Table 1, it is possible to identify the correlation between the property values and the unique mechanical behavior of each specimen and, therefore, establish the relationship between material properties and material composition. For instance, the 85–00 specimen is accurately modeled by two endochronic branches, while 85–15 and 85–30 specimens require three endochronic branches to fully represent their mechanical response, which signifies the higher material nonlinearity and more extensive inelastic behavior of the aluminized specimens. Specific trends in individual property values are also observed, which are discussed in detail below.

In the context of representing the deformation mechanics of particle-binder composites by a linear rheological model (Fig. 10) with respect to the constitutive formulation presented in Sec. 2, the elastic branch can be seen as a spring of stiffness 2μ^e acting in parallel with N endochronic branches, each of which consists of a spring of stiffness $2μji$ and a friction element with threshold strain γ_j, where j = 1, …, N. The total stiffness μ of the material is represented by a sum of the elastic and endochronic moduli, i.e., $μ=μe+∑j=1Nμji$⁠, while the threshold strains represent the macroscopic strains at the activation of inelastic processes, such as irreversible slipping of binder molecular chains on the filler (solid) surface and plastic deformation of the solid particles [46,48].

From the calibrated parameter values listed in Table 1, the total stiffness of 85–00, 85–15, and 85–30 specimens is found to be 5.1941 MPa, 12.3132 MPa, and 8.2868 MPa, respectively. This suggests that the addition of aluminum in the PBX composition increases material stiffness; however, the stiffness starts to reduce with an increase in the amount of aluminum in the composition. This postulation is also supported by the evolution of the apparent modulus during the cyclic loading of the three compositions (Fig. 9). The figures show that the apparent modulus of the cyclic loops ranges between 5 MPa and 7 MPa for 85–00, 6 MPa and 11 MPa for 85–15, and 4 MPa and 11 MPa for 85–30. Therefore, the experimentally determined average modulus follows the same order with respect to the material composition as the model-estimated modulus. Another interesting observation lies in the contribution to the total modulus by the endochronic branches, which from Table 1 is 0.0013 MPa (0.025%) for 85–00, 8.2364 MPa (66.89%) for 85–15 and 7.5831 MPa (91.51%) for 85–30. Therefore, it is evident that with the activation of inelastic deformation mechanisms, the elastic material stiffness decreases with increasing aluminum content, which is also consistent with the observation of predominantly nonlinear mechanical response and ductile plastic flow to larger strain levels during monotonic compression in 85–15 and 85–30 as compared to 85–00 (Fig. 2).

With regard to the threshold strains γ_j, the value of strain γ₁ in the first endochronic branch is very small and comparable in magnitude for all PBX compositions (0.0017 for 85–00, 0.001 for 85–15, and 0.0018 for 85–30). This suggests that mock PBX material starts exhibiting inelastic behavior from almost the beginning of the loading process regardless of aluminum content, which rationalizes the selection of a yield surface-free plasticity theory to model these materials. However, the value of threshold strain in the second endochronic branch (γ₂) is much larger for 85–00 (0.1723) as compared to those of 85–15 (0.0358) and 85–30 (0.0622), indicating that more extensive inelastic processes start occurring at smaller macroscopic strains in the aluminized specimens. When comparing 85–15 and 85–30 specimens, we observe that although the second threshold strain is smaller for 85–15, the threshold strain in the third branch (γ₃) is smaller for 85–30 (0.2388) as compared to 85–15 (0.3413). Therefore, it can be concluded that higher aluminum content leads to the material exhibiting predominant inelastic deformations at a smaller strain. This conclusion is also supported by observations of higher residual strain in 85–30 (∼2%) as compared to 85–15 (∼1.57%) after the recovery period following the virgin cyclic loading tests as observed by Agarwal and Gonzalez [23].

As explained in Sec. 2.5, the material scale function f(z), which is a function of the intrinsic time scale z, is used in the constitutive model to represent the progressive cyclic stress-softening behavior observed during cyclic compression of particle-binder composites. The function is given by Eq. (31) and consists of four material parameters, namely, c, s, β_c, and β_s. The parameter c is the saturation value of the function during initial loading from zero deformation when the reference intrinsic time scale, z_ref, is equal to zero. If the loading remains monotonic, the simpler scale function proposed by Wu and Yip [59] and given by Eq. (30) is recovered. According to this simpler function, the value of f(z) evolves from 1 as z increases and approaches c at z → ∞, with the rate of evolution governed by β_c. If c > 1, then f(z) ≥ 1. Therefore, during a loading step, the increment in the intrinsic time scale z is greater than the increment in deformation measure $|C|$ according to Eq. (29), which produces a strain-hardening behavior. Conversely, if c < 1, a strain-softening response is obtained. This is also true during cyclic loading–unloading; however, the saturation value changes as z_ref increases, and approaches c/s at z_ref → ∞ which signifies cyclic stabilization. If c/s < c, then a progressive stress-softening response is obtained, which, as already stated, is one of the primary deformation mechanisms observed in mock PBX.

The value of parameter c for the studied mock PBX specimens is seen to increase with increasing aluminum content in the composition. Specifically, it is much higher for 85–30 as compared to 85–00 and 85–15 specimens. The reason for this difference lies in the monotonic compression response of the specimens. Figure 11 shows the influence of variation in c from its estimated value on the monotonic compression response of 85–00 specimen. It is evident that c controls the occurrence of maximum or ultimate stress in the response, as well as the magnitude of the ultimate stress, which stems from its representative role as the magnitude of strain hardening in the material. The ultimate stress is observed in the monotonic response of 85–00 and 85–15 specimens, with the value of the stress and the strain at which it occurs is higher for 85–15 as compared to that for 85–00. Therefore, the estimated value of parameter c for 85–15 is comparatively higher than for 85–00. However, the 85–30 specimen does not show an ultimate stress and continues to harden until the maximum applied deformation of $30%$⁠; hence, in accordance with Fig. 11, the estimated value of c for 85–30 is much higher than both 85–00 and 85–15. With regard to parameter β_c, a correlation between the estimated values and PBX aluminum content is not observed; rather, the values are closely related to the evolution of the scale function during monotonic loading and the estimated values of c. Figure 12 shows the influence of variation in β_c on the monotonic compression response as well as the evolution of scale function f(z) with intrinsic time z for the three specimens. It is evident that for 85–00 and 85–15 specimens, the scale function f(z) reaches its saturation value c during the monotonic loading process, while the same does not occur for the 85–30 specimen as the material continues to harden beyond the maximum applied strain. For 85–15 (Figs. 12(c) and 12(d)), the saturation value is reached much before the attainment of ultimate stress, and therefore, any variation in β_c has only a slight influence on the loading response at small strains (∼1%–5%). For 85–00 (Figs. 12(a) and 12(b)), the saturation value is reached near the end of the loading process, and therefore, the influence of variation in β_c on the stress–strain response is comparatively larger than 85–15. For 85–30 (Figs. 12(e) and 12(f)), however, any variation in β_c has a significant impact on the stress–strain response. It is specifically observed that β_c controls only the evolution of stress, with the stress at any strain increasing with increasing β_c and vice versa. With regard to the magnitude of β_c, a higher value of β_c increases the rate of evolution of f(z) and therefore the rate of attainment of the saturation value. Hence, 85–15 has the highest value of β_c (944.2924), followed by 85–00 (125.3399), while 85–30 assumes a much lower value (2.0906).

As discussed before, the role of parameters s and β_s in the context of PBX material is to reduce the saturation value of the scale function f(z) from c to c/s during the cyclic loading process in order to achieve progressive stress-softening until cyclic stabilization. While s primarily controls the magnitude of the softening, β_s controls the softening rate or the rate of attainment of cyclic stabilization. From the estimated values of s for the three specimens, the value of c/s is obtained as 11.9027, 14.4811, and 99.0748 for 85–00, 85–15, and 85–30 respectively, resulting in a reduction of 20.58%, 50.94%, and 83.24% in their respective saturation values at cyclic stabilization as compared to c. This implies that the magnitude of softening is the highest for 85–30, followed by 85–15, and the least for 85–00. Concurrently, the estimated values of β_s for the three specimens indicate the highest rate of stabilization for 85–00 (212.6798), followed by 85–15 (16.7109), and the lowest for 85–30 (5.7283). These findings with reference to the model parameters are in agreement with the experimental data, as shown in Fig. 13, where the ratio of cyclic peak stresses with respect to the first peak stress (Fig. 13(a)) drops the highest for 85–30, followed by 85–15 and then 85–00, while the change in peak stress with each cycle (Fig. 13(b)) approaches zero (i.e., cyclic stabilization) fastest for 85–00, followed by 85–15 and then 85–30.

As shown in Secs. 2.1 and 2.3, the primary function of the isotropic damage model is to represent the overall softening between loading and unloading responses during cyclic loading of particle-binder composites. The model consists of two positive parameters m and p, each of which, according to Eq. (17), increase (decrease) the softening effect as their values decrease (increase). Figure 14 shows the influence of variation in m and p from their estimated values on a simulated load–unload response from 16% strain for the 85–15 formulation. It is observed that p predominantly affects the softening rate at higher strains, more specifically, near the beginning of unloading. On the contrary, m exerts a predominant control at lower strain levels. A closer observation of Fig. 14 indicates that the effects of p are dominant during the initial rapid softening (∼15–16% strain in the figure), while the influence of m predominates thereafter when the softening rate decreases. Figure 15 shows the evolution of the slope of the first cyclic unloading curve for the experimental cyclic compressive response of the three mock PBX specimens. The slope is plotted for the first 2% strain after the start of unloading since it is seen that within this range the unloading curves of all three specimens begin to curve and the softening rate declines. Therefore, the effects of both m and p can be studied within this range for the three specimen compositions. From the figure, it is evident that the 85–30 specimen begins unloading with the highest slope and, therefore, exhibits the highest rate of initial softening, followed by 85–15 and then by 85–00. This is in agreement with the decreasing trend of the estimated p values in Table 1 with increasing aluminum content. With further unloading, the slope of the stress–strain curve drops rapidly for all specimens and then approaches a stable value, indicating more linear stress–strain response beyond ∼0.5% unloaded strain. In this region, the influence of parameter m is predominant and a lower slope indicates a higher softening effect (Figure 14(a)). Therefore, 85–30 again experiences the highest softening in this strain range since it has the lowest slope. However, the slope of 85–00 is lower than 85–15 in this range and, therefore, exhibits a higher softening. This is again in agreement with the estimated values of m in Table 1, which is the highest for 85–15, followed by 85–00 and lowest for 85–30.

We have presented a constitutive model for particle-binder composites that accounts for finite-deformation kinematics, nonlinear elastoplasticity without apparent yield, cyclic hysteresis, and progressive stress-softening before the attainment of stable cyclic response. The model is based on an additive decomposition of strain energy into elastic and inelastic parts, where the elastic response is modeled using Ogden hyperelasticity while the inelastic response is described using yield-surface-free endochronic plasticity based on the concepts of internal variables and of evolution or rate equations. Stress-softening is modeled using two approaches; a discontinuous isotropic damage model to appropriately describe the overall softening between loading and unloading responses, and a material scale function to describe the progressive cyclic softening until the attainment of stable response. The constitutive model is then based on the deformation mechanisms experimentally observed during cyclic compression of PBX at large strain.

Furthermore, we have proposed a discrete numerical procedure to solve for stresses along a loading path, and we have developed a parameter identification method, based on a nonlinear multivariate minimization problem, to determine material properties from experimental nominal stress–strain data. Specifically, we have used monotonic and cyclic compression data for three mock PBX formulations based on PBXN-109 with varying aluminum content [23] to demonstrate the range of behavior predicted by the proposed constitutive model and the effectiveness of the parameter identification method. This is evident from the remarkable model calibration results that yielded an excellent agreement between experimental data and model response. Finally, we have identified the correlation between the mechanical properties and the unique stress–strain response of each specimen, which contributes toward establishing and validating specific trends observed in the property values with respect to the specimen composition, specifically, the amount of aluminum in the mock PBX. We, therefore, conclude that each of the estimated material parameters is more than a mere number and has a significant contribution toward the establishment of the highly complex and nonlinear mechanical behavior of particle-binder composites under finite strain. Based on these conclusions, the constitutive model presented in this paper serves as a reliable tool to characterize different compositions of particle-binder composites, and it especially serves the defense and propulsion communities in their efforts to characterize the mechanical properties of both existing and new energetic material formulations.

We close by pointing out some limitations of our approach and possible directions for the future extension of our analysis. First, we have restricted our attention to quasi-static, i.e., low strain-rate behavior of particle-binder composites, and thereby neglected any strain-rate-dependent effects in the constitutive model. However, it is well known that viscous effects become dominant in these materials at moderate to high strain rates. Second, the proposed constitutive model carries an assumption of elastoplastic incompressibility, which limits the analysis to particle-binder composites with moderate levels of compressibility. Therefore, the extension of the model to rate-dependent viscous effects and compressible behavior is a worthwhile direction for future research.

The authors wish to acknowledge Professor Patricia Davies, Professor Jeff Rhoads, Professor Steve Son, Professor Stuart Bolton, Jelena Paripovic, Allison Range, Nick Cummock, and Tim Manship for their feedback and interesting discussions.

There are no conflicts of interest.

The data sets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request. Data provided by a third party are listed in Acknowledgment.

This research is supported by the Air Force Research Laboratory through Grant No. FA8651-16-0287 entitled “Near-Resonant Thermomechanics of Energetic and Mock Energetic Composite Materials,” Grant No. FA8651-16-D-0287 entitled “Near-Resonant Thermomechanics of Energetic and Mock Energetic Composite Materials, Part II,” and Grant No. FA8651-17-S-0003 entitled “Exploring the Thermomechanics of Energetic and Mock Energetic Composite Materials Under Quasi-Static and Near-Resonant Excitations.”

Approved for public release. Distribution unlimited (AFRL-2021-2732).