Peridynamic modeling of dynamic damage of polymer-bonded explosive (PBX) with sensitivity analysis

Abstract

Polymer-bonded explosive (PBX) is one of the typical heterogeneous composite materials and dynamic damage processes are crucial for understanding its behavior. The present study establishes a comprehensive peridynamic (PD) framework to model dynamic damage of PBX with sensitivity analysis of model parameters. The PD model successfully captured the various characteristics of PBX damage evolution: energetic crystal damage dominated under low-velocity impacts (20 m/s), while interfacial damage contributions increased substantially at higher velocities (40–60 m/s). Moreover, the importance of the total nine PD model parameters belonging to three different materials regarding total bond damage and damage modes (trans-granular/inter-granular) under 20–60 m/s impact loading is quantified based on the established three-dimensional PD model of steel-encased PBX. The Kriging surrogate model from PD model parameters to PBX damage is also constructed to decrease the computational cost. The sensitivity analysis revealed that binder failure strain exerted significant influence on total mechanical damage at different velocities. Notably, Young's modulus of energetic crystal demonstrated strong cross-mode effects at elevated velocities, governing both trans-granular crack initiation and inter-granular debonding, thereby emerging as a critical constitutive parameter for relatively high impact velocity. The current study offers a useful computational framework to evaluate how the uncertainty of PBX properties impacts its dynamic damage response.

Keywords

Peridynamics dynamic damage polymer-bonded explosives sensitivity analysis

Introduction

Polymer-bonded explosive (PBX) is typically a heterogeneous composite material composed of energetic crystals such as HMX (octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine) and RDX (1,3,5-trinitroperhydro-1,3,5-triazine), dispersed within polymer binders like hydroxyl-terminated polybutadiene and Estane. They are widely used in both military and industrial fields because of their combination of good detonation performance and mechanical properties (Kim et al., 2018). PBX may be exposed to various external mechanical stimuli during transportation, storage, or handling (Gruau et al., 2009), and low-velocity impact represents a common yet potentially hazardous scenario. Although such impacts usually involve relatively low energy levels, they may still induce the unintended explosion of PBX due to the formation of hotspots caused by the localized energy deposition (Field, 1992), thereby posing significant safety risks due to accidental initiation.

Experimental observations also show that explosives are first squeezed, and the occurrence of chemical reaction and generation of hotspots only after being squeezed to a certain extent under impact loads (Napadensky et al., 1960). Understanding the mechanical damage process and its mechanisms is important for revealing the failure mechanism of PBX and improving its safety. Numerical simulations can reveal the dynamic processes of crack propagation and hotspot formation. Finite-element methods are previously often used to simulate the dynamic response of PBX, and the simulation results can help analyze the internal evolution processes of PBX micro-structures that are difficult to observe experimentally, as well as the influence of different material parameters on the ignition mechanism (Barua et al., 2012a, 2013). Bennett et al. (1998) developed and demonstrated the use of the viscoelastic response and statistical fracture mechanics in the numerical constitutive model of PBX-9501 to predict its mechanical response. Since the extension of the crack surface plays a key role in ignition, Dienes et al. (2006) further developed statistical crack models to explain the characteristics in crack growth and aggregation. Xiao et al. (2024) also developed a mechanical–thermal constitutive relation including viscoelastic responses to predict the mechanical response and thermal response of PBX. Because the singularity at the crack tip is difficult to capture precisely, a finite-element method has certain limitations in crack simulation (Madenci and Oterkus, 2013). Improvements such as extended finite-element method (Wyart et al., 2009) and cohesive finite-element method (Barua et al., 2012a, 2012b, 2013) provide a better description of crack-related problems, but these methods have significantly increased complexity and computational efficiency.

Peridynamics (PD), proposed by Silling (2000) and Silling et al. (2007), is a continuum mechanics model based on non-local theory that overcomes the assumption in traditional continuum mechanics that the displacement field needs to be continuously differentiable by introducing an integral form of the motion equation instead of a differential equation. The theory can naturally describe the spontaneous initiation and dynamic propagation of discontinuous problems such as cracks in materials. Compared with traditional numerical methods, PD can naturally describe crack initiation and propagation without special treatment of crack-tip singularities, pre-defined crack paths, or additional crack-tracking criteria. It is therefore suitable for simulating complex fracture processes involving crack nucleation, branching, coalescence, and interaction. This feature is particularly useful for PBX, whose dynamic damage behavior involves the coupled evolution of crystal fracture, binder deformation, and interfacial debonding. PD formulations can generally be classified into bond-based and state-based models (Ahmadi et al., 2022). A classical bond-based PD model is subject to a fixed Poisson's ratio limitation, which is 1/4 for three-dimensional isotropic materials. A state-based PD model provides a more general constitutive framework and can overcome the fixed Poisson's ratio limitation of the classical bond-based formulation. However, a bond-based PD model is simpler and more efficient in numerical implementation, which is advantageous for repeated simulations in parametric studies. In the present work, a bond-based PD model is adopted.

Currently, PD has been widely applied to various areas such as glass (Ha and Bobaru, 2011) and composites (Diyaroglu et al., 2016). Therefore, PD is well suited for investigating the dynamic damage behavior of PBX under impact loading, where multiple damage modes may coexist and evolve simultaneously. Moreover, Deng et al. (2024) constructed a mechanical–thermal–chemical coupled multi-physics PD framework to simulate the damage and ignition response of PBX under low-velocity impact loads. The study analyzed the ignition mechanisms related to trans-granular and inter-granular damage modes, providing key information for establishing an ignition model. Huang et al. (2022) investigated the dynamic damage behavior of PBX confined in spherical steel shells using bond-based PD, highlighting the influence of impact velocity and steel-shell confinement on PBX damage. Prakash and Seidel (2018) proposed a numerical simulation of coupled electromechanical formulas based on PD non-local continuum theory to explore structural health monitoring of PBX based on nano-composite piezo resistance rates. Li and Hao (2023) investigated the elastic deformation and brittle fracture behavior of PBX simulators using a modified three-body potential PD model.

The dynamic damage of PBX involves complicated processes such as particle–binder interface debonding and micro-crack nucleation and propagation, thus the simulation results are highly dependent on material and structural parameters (e.g. elastic modulus, adjustable interface parameters, loading rate, etc.) (Deng and Wang, 2020; Huang et al., 2022). As a heterogeneous brittle material, the particle distribution, porosity, and explosive grain content of PBX may lead to differences in mechanical properties. Similar micro-structure effect analyses have also been reported for polymer composites. For example, Huang et al. (2024) investigated carbon black-filled rubber composites by combining mechanical experiments with finite-element representative volume element simulations, and showed that filler volume fraction, particle orientation, and agglomeration can significantly affect the stress–strain response and local stress concentration of polymer composites. In addition, the mechanical response of PBX is also related to other factors such as strain rate, temperature, confinement pressure, and so on. Overall, it is difficult to accurately determine the specific values of model parameters. Therefore, even for the same PBX, different values of simulation parameters are used in the previously published literature.

Therefore, it is first necessary to consider the influence of uncertainty in model parameters on the simulation results regarding the simulations of dynamic mechanical response of PBX. Moreover, it is also necessary to conduct sensitivity analysis on parameters to systematically quantify the extent to which different parameters affect the numerical simulation results (Jiang et al., 2016). By doing this, the coupled relationship between the material model parameters and the mechanical response of PBX can be understood more deeply. Furthermore, the simulation results can also be better used to design the experiments or explain the experimental results. Uncertainty quantification has already been introduced into the PD framework in previous studies. For example, Franzelin et al. (2014) applied a non-intrusive uncertainty quantification method based on spatially adaptive sparse-grid stochastic collocation to multi-variate PD simulations, showing that uncertainty propagation in PD models can be efficiently analyzed without directly relying on large-scale Monte Carlo simulations. In addition, Kriging-based surrogate modeling combined with Sobol sensitivity analysis has been successfully used in other damage and mechanical systems. Denimal and Sinou (2021) developed advanced Kriging-based surrogate modeling and sensitivity analysis for uncertain rotor dynamics. And the similar surrogate-based sensitivity strategies (Li et al., 2025) have also been applied in coupled damage modeling contexts. These studies provide important methodological support for using surrogate models and global sensitivity analysis in complex damage simulations.

However, existing uncertain quantification and sensitivity studies are still limited for the specific problem of PBX dynamic damage under impact loading. In particular, few studies have systematically combined PD simulation, Kriging surrogate modeling, damage-mode quantification, and Sobol sensitivity analysis to evaluate how uncertain PD model parameters affect total damage and the transition between trans-granular and inter-granular damage modes. Unfortunately, as far as the uncertainty and sensitivity analysis of PBX dynamic damage is concerned, the previous research is very limited and has mainly focused on the effects of randomness of micro-structure on the simulation results (Kim et al., 2014). The lack of systematic analysis on the uncertainty and sensitivity of PBX dynamic damage inhibits the further progress of numerical simulations.

In this study, the influence of uncertainty in model parameters involved in PD models for PBX impact damage and damage modes was studied in a systematic manner. By constructing a Kriging surrogate model, the parameters that have a significant impact on damage modes and total bond damage were identified based on the Sobol global sensitivity analysis method. The current results are of great significance for an in-depth understanding of the PBX dynamic damage behavior.

PD theory and sensitivity analysis of model parameters

PD theory

PD is a novel non-local continuum mechanics theory based on integral equations, suitable for modeling discontinuous problems such as material damage and crack propagation. PD theory discretizes a continuous medium into a series of material points, each with mass, position, and volume properties. The material points interact with each other through a response function that contains all constitutive information related to the material. Figure 1 shows a schematic diagram of the PD theory, where the material points that can interact are represented by $H_{x}$ , with a radius of $δ$ . The mechanical behavior of each material point depends not only on directly adjacent points but also on the collective action of all points within the horizon. In PD, the interactions between material points are referred to as “bonds,” and the spatial non-locality of these interactions is governed by the horizon parameter, which defines the maximum interaction range of a material point within its neighborhood. As the horizon shrinks, the interaction becomes more localized. Therefore, the classical elasticity theory can be regarded as a limit case of the PD theory when the horizon approaches zero.

Figure 1.

Schematic diagram of PD theory. PD: peridynamics.

The basic equation of motion of PD is:

ρ \ddot{u} (q, t) = \int_{H_{x}} f (u (q, t), u^{'} (q^{'}, t), q, q^{'}, t) d V_{q^{'}} + b (q, t)

(1)

where

ρ

is the density,

H_{x}

is the horizon with radius

δ

as the center q,

f

is the force density function of the bond,

u

is the displacement vector (i.e.

\ddot{u}

is the corresponding acceleration vector), and

b

is the volume force. Because the right-hand side of the equation is in the integral form, the PD theory avoids the singularity that discontinuity may cause in the numerical simulation process. The horizon is usually a sphere centered on the particle with a radius of three times the particle length (

δ = 3 Δ q

). PD describes the interaction between material points through the force density function of the bond, so the construction of the force density function

f

is crucial. The commonly used force density function is the prototype micro-elastic brittleness model, and its constitutive model is shown as follows:

f = c s μ (t, ξ) \frac{ξ + η}{‖ ξ + η ‖} V_{c o e}

(2)

s = \frac{‖ ξ + η ‖ - ‖ ξ ‖}{‖ ξ ‖}

(3)

where c is the micro-modulus, s is the relative elongation, and

μ

is the damage function.

V_{c o e}

is a volume correction factor, and the

V_{c o e}

can be expressed as:

V_{c o e} = {\begin{matrix} (\frac{δ - ‖ ξ ‖}{Δ x} + \frac{1}{2}) \begin{matrix} if (δ - 0.5 Δ x) \leq ‖ ξ ‖ \leq δ \end{matrix} \\ \begin{matrix} 1 \\ 0 \end{matrix} \begin{matrix}  \end{matrix} \begin{matrix} if ‖ ξ ‖ \leq (δ - 0.5 Δ x) \\ otherwise \end{matrix} \end{matrix}

(4)

where

ξ

represents the size of the relative position of the undeformed material point and the adjacent material point (

ξ = q^{'} - q

). The relative displacement of the two material points after the deformation of the material point and can be used to represent

η = u^{'} - u

. So

ξ + η

can represent the size of the relative position of the current two material points after the deformation. The micro-modulus c can be determined based on the equivalence of strain energy density with that of the classical continuum. In three-dimensional simulation, c is a function of Young's modulus E and horizon radius

δ

given by the following formula:

c = \frac{12 E}{π δ^{4}}

(5)

In PD theory, damage to a material point is represented as bond breakage. When the elongation s is greater than the critical relative elongation $S_{c}$ of the bond, it indicates a bond breakage between two material points, and the bond force between the two points will not be calculated during the simulation. The damage function is given by (Silling and Askari, 2005):

μ (t, ξ) = {\begin{matrix} 1, & if s < S_{c}, 0 < t^{'} < t \\ 0, & otherwise \end{matrix}

(6)

where the critical elongation

S_{c}

is related to the fracture energy

G_{0}

, which represents the energy per unit fracture area required for all the bonds to break on the fracture surface, usually

s_{c} = \sqrt{5 G_{0} / 6 E δ}

To measure the damage state of the material point, the local damage function $φ$ is defined as:

φ (q, t) = 1 - \frac{\int_{H} μ (t, ξ) d V^{'}}{\int_{H} d V^{'}}

(7)

According to the above definition, when the damage function value is 0, it indicates that the material point has bond forces with other material points in the horizon, that is, the material point is in an undamaged state. In contrast, when the damage function value is 1, it indicates that all bonds associated with the material point have been broken, that is, the material point is in a severely damaged state. If cracks appear inside the material, the object will split into two parts along the crack surface, meaning that approximately half of the bonds within the neighborhood of material points near the crack boundary will fail. Accordingly, a local damage value of about 0.5 is often used in the PD literature as an indicator of possible crack evolution (Zhang and Bobaru, 2016).

Global Kriging surrogate model

When conducting sensitivity analysis on the constitutive model parameters regarding the mechanical damage of PBX, it is necessary to quantify the uncertainty of damage. This can be achieved by propagating the uncertainty of constitutive model parameters to the mechanical damage through uncertainty propagation methods such as Monte Carlo simulation. However, PD needs to calculate all material points and the interaction forces within the horizon of that material point at each time step. The process is usually time-consuming, and therefore it is not realistic to conduct PD simulations based on a large number of parameter samples across the entire domain. As a statistical model, the Kriging surrogate model technique can fit mathematical models to a limited number of observation points as approximations. The low-cost Kriging surrogate model can significantly reduce the computational time for uncertainty propagation analyses, such as Monte Carlo sampling (Liu et al., 2025).

The core assumption in the theoretical framework of the Kriging surrogate model is to reinterpret the originally deterministic output function ( $g (X)$ ) as a concrete implementation of a random process ( $g_{k} (X)$ ). Figure 2 shows how the method models the output response by combining the global trend model and the local deviation model. The global trend model, which is constructed using polynomials or other basis functions such as linear or quadratic polynomials, is used to capture the overall variation pattern of the output response at the macro level. The local deviation model models the residual of the trend model through the spatial correlation function, characterizing the local deviation between the trend prediction and the actual observation point.

Figure 2.

Schematic diagram of the Kriging surrogate model.

For a single-output mathematical model $Y = g (X)$ , where $X$ is an n-dimensional input variable and Y is a one-dimensional output variable, the Kriging surrogate model used to approximate $Y = g (X)$ can be represented as:

g_{k} (X) = τ (X) + ε (X)

(8)

where

τ (X)

represents the global trend component of the Kriging surrogate model, also known as the regression predictor. The regression model

τ (X) = h^{T} (X) β = β_{1} h_{1} (X) + \dots + β_{l} h_{l} (X)

is a linear combination of l selected basis functions

h (X)

and coefficients

β

. For the local spatial deviation model,

ε (X)

is generally a Gaussian distribution with 0 mean and standard deviation

σ

. Its covariance is given by:

C o v [ε (x_{i}), ε (x_{j})] = σ^{2} R_{X} (x_{i}, x_{j})

(9)

where

ε (x_{i})

and

ε (x_{j})

denote the residuals at different input samples. And

R_{X}

is the correlation function between input samples, typically in Gaussian form:

R_{X} (x_{i}, x_{j}) = \exp (- \sum_{p = 1}^{n} θ_{p} {| x_{i}^{p} - x_{j}^{p} |}^{2})

(10)

where

θ_{p} (p = 1, \dots, n)

are the hyper-parameters of the Kriging surrogate model, and

x_{i}^{p}

denotes the

p

-th dimension of the sample

x_{i}

Given $N_{k}$ training samples $X_{t r a i n} = {x_{1}, \dots, x_{N_{k}}}$ and the corresponding outputs $y = {y_{1}, \dots, y_{N_{k}}}$ , the following likelihood function can be constructed to determine the Kriging surrogate model parameters to minimize the mean squared prediction error:

L (β, σ^{2}, θ) = - \frac{N_{k}}{2} \ln (2 π σ^{2}) - \frac{1}{2} \ln (| R_{X} |) - \frac{1}{2 σ^{2}} (y - H β)^{T} R_{X}^{- 1} (y - H β)

(11)

where

R_{X}

is the correlation matrix constructed from the input training samples, and

H = [h (x_{1}), \dots h (x_{N_{k}})]^{T}

is the regression function matrix based on the input training samples.

By performing maximum likelihood estimation on the above likelihood function, the estimates $\hat{β}$ , ${\hat{σ}}^{2}$ , and $\hat{θ}$ of the Kriging surrogate model parameters can be obtained. The explicit form of the Kriging surrogate model can then be written as:

g_{k} (X) = h^{T} (X) \hat{β} + r^{T} (X) R_{X}^{- 1} (y - H \hat{β})

(12)

where the vector

r (X) = [R_{X} (X, x_{1}), \dots, R_{X} (X, x_{N_{k}})]^{T}

denotes the correlation terms between an arbitrary input

X

and each training sample.

When building a Kriging surrogate model, a large number of samples need to be drawn as a candidate sample pool $X_{t o t a l}$ based on the distribution and range of the uncertain parameters, but the actual training sample set $T$ for building a Kriging surrogate model is a portion of the samples drawn from the candidate sample pool. As can be seen from Figure 2, the Kriging surrogate model deviates from the predictions at different points, so the key to building the mathematical surrogate model of the original model is whether the points with the greatest deviation in the candidate sample pool can be added to the training samples; ultimately, the accuracy of the surrogate model at each input variable sample point can meet the requirements.

Recently, the adaptive sampling strategy has received increasing attention. This sampling strategy selects new points sequentially based on an approximate model and information from the data itself of previous iterations. The overall procedure of the adaptive sampling strategy used in this study is illustrated in Figure 3. The autonomous learning function based on the leave-one-out cross-validation (LOOCV) can effectively select the training samples that contribute the most to fitting the original model (Liu et al., 2017). Extracting $N_{k}$ training samples from $X_{t o t a l}$ as the training set $T - = [(x_{1}, y_{1}), \dots, (x_{N_{k}}, y_{N_{k}})]$ , the basic principle of LOOCV is to set the samples $(x_{i}, y_{i})$ ( $i$ is a positive integer on the interval $[1, N_{k}]$ ) as the validation set, and then take $N_{k} - 1$ samples $T -^{\sim i} = (x_{\sim i}, y_{\sim i})$ ( $\sim i$ denotes all indices except $i$ ) as the training set to build the Kriging surrogate model ${\hat{g}}_{k}^{\sim i} (\cdot)$ . The constructed Kriging surrogate model calculates the deviation between the predicted value ${\hat{g}}_{k}^{\sim i} (x_{i})$ and the true value $y_{i}$ at $x_{i}$ , which serves as an evaluation metric for the surrogate model. This deviation is termed the leave-one-out error $E^{L O O}$ :

E^{L O O} (x_{i}) = (y_{i} - {\hat{g}}_{k}^{\sim i} (x_{i}))^{2} (i = 1, \dots, m)

(13)

Figure 3.

Flowchart of the adaptive sampling strategy for Kriging surrogate model construction.

Combining the $E^{L O O}$ for each sample can determine whether the accuracy of the current Kriging surrogate model still needs to be fitted by expanding the training sample size. However, this process needs to be repeated m times to obtain the $E^{L O O}$ for each sample. The calculation also becomes difficult when the sample base increases. Sundararajan and Keerthi (2001) avoided the computational cost of re-estimating the correlation parameters in the Kriging surrogate model under the same sample pool by using the random correlation matrix output from the Kriging surrogate model and the regression model parameters. A $E^{L O O}$ training set is formed $T_{c} = [(x_{1}, E^{L O O} (x_{1})), \dots, (x_{N_{k}}, E^{L O O} (x_{N_{k}}))]^{T}$ after obtaining the $E^{L O O}$ for each training sample. By treating the training sample errors as outputs, the Kriging surrogate model $E_{k}^{L O O} (x)$ constructed from the training set $T_{c}$ can be used to compute the $E^{L O O}$ for the remaining samples $X_{r e m a i n} = X_{t o t a l} - T$ in the candidate pool. The point with the maximum predicted value $max (E_{k}^{L O O} (X_{r e m a i n}))$ , which indicates thatthe current Kriging surrogate model may have a greater prediction deviation at that point, is then selected as the next candidate sample. Update the current Kriging surrogate model by adding the sample with the largest residual error value in the candidate sample pool as the best point to the training sample set until the root mean square error of the Kriging surrogate model at the training sample is less than the set value.

Sobol global sensitivity analysis

Sobol global sensitivity analysis method (Sobol, 2001) has remained widely used since its introduction because it applies to any input and output. The Sobol indices it produces offer robust and reliable assessments of parameter sensitivities over a range of variations in input parameters without being affected by minor changes. In this study, Sobol indices are employed to represent the importance of each parameter. The steps for calculating Sobol indices are described below.

The evaluation function can be described as a numerical simulation process with n-dimensional inputs and one-dimensional output:

Y = g (X)

(14)

where

X = [X_{1}, X_{2}, \dots, X_{n}]^{T}

and Y are the input and output vectors, respectively. Based on the variance decomposition formula, the first-order Sobol index

S_{i}

and the total index

S_{i}^{T}

X_{i}

can be represented as:

S_{i} = \frac{V_{i}}{V} = \frac{V a r [E (Y | X_{i})]}{\sum_{i = 1}^{n} V_{i} + \sum_{i_{1} = 1}^{n} \sum_{i_{2} = i_{1} + 1}^{n} V_{i_{1} i_{2}} + \dots + V_{1, 2, \dots, n}}

(15)

S_{i}^{T} = 1 - \frac{V_{\sim i}}{V} = 1 - \frac{V a r [E (Y | X_{\sim i})]}{\sum_{i = 1}^{n} V_{i} + \sum_{i_{1} = 1}^{n} \sum_{i_{2} = i_{1} + 1}^{n} V_{i_{1} i_{2}} + \dots + V_{1, 2, \dots, n}}

(16)

In the above equations, V denotes the total (unconditional) variance, and $V_{i}$ , $V_{1, 2, \dots, n}$ , and $V_{\sim i}$ are the conditional variance associated with variables, where $i \in [1, n]$ and $\sim i$ refers to all variables excluding the $i$ -th variable. Moreover, $V a r (\cdot)$ represents the variance operator and $E (\cdot)$ is the expectation operator. The first-order Sobol index $S_{i}$ represents the contribution of a single input variable alone to the output variance, and the total index $S_{i}^{T}$ represents the sum of the interaction indices of all variables associated with that input variable, which is an important reference metric. When the variables are correlated, such as HMX–HMX, HMX–binder, and binder–binder, the number of broken bonds is a function of the PD model parameters which have a certain connection. At this point, the correlation between the variables should be considered. The first-order Sobol index $S_{i}$ should decompose the total contribution rate of the input variable to the variance of the output response quantity into independent contribution $S_{i}^{U}$ and correlation contribution $S_{i}^{C}$ according to the theory proposed by Xu and Gertner (2008):

S_{i} = S_{i}^{U} + S_{i}^{C}

(17)

where the independent contribution rate and the related contribution are calculated by the following formulae:

S_{i}^{U} = 1 - \frac{V a r (E (Y | X_{\sim i}))}{V}

(18)

S_{i}^{C} = S_{i} - S_{i}^{U}

(19)

In the present study, the three correlated variables are taken as the broken numbers of HMX–HMX, HMX–binder, and binder–binder bonds, denoted by $Y_{1}$ , $Y_{2}$ , and $Y_{3}$ , respectively, and the total damage is defined as:

Y_{t} = Y_{1} + Y_{2} + Y_{3}

(20)

To implement the above decomposition numerically, 100,000 samples of the uncertain PD model parameters were generated using a Sobol low-discrepancy sequence, and the trained Kriging surrogate model was employed to predict the corresponding values of $Y_{1}$ , $Y_{2}$ , and $Y_{3}$ . For each bond type $Y_{i}$ , its total contribution was estimated from the fitted response of $Y_{t}$ regressed on $Y_{i}$ , and its independent contribution was estimated by first removing the part of $Y_{i}$ explained by the other two bond types, that is:

Z_{i} = Y_{i} - {\hat{Y}}_{i} (Y_{\sim i})

(21)

where

Z_{i}

denotes the residual component of

Y_{i}

that cannot be explained by the other damage variables. The total damage

Y_{t}

was then regressed on

Z_{i}

, and the variance of the corresponding fitted response was taken as the independent contribution. Accordingly, the total contribution, independent contribution, and correlation contribution were calculated as:

S_{i} = \frac{V_{i}}{V}, S_{i}^{U} = \frac{V_{i}^{U}}{V}, S_{i}^{C} = S_{i} - S_{i}^{U}

(22)

where

V = V a r (Y_{t})

V_{i}

is the variance of the fitted response associated with the total contribution of

Y_{i}

, and

V_{i}^{U}

is the variance of the fitted response associated with the independent contribution of

Y_{i}

. It should be pointed out that calculating the independent contribution rate and the related contribution rate of each type of broken bond separately can help determine whether each type of damage affects the total damage more individually or collectively, and also helps analyze the coupling relationship among different damage types.

Simulation model and parameters

Dynamic damage model of PBX

The PBX dynamic damage calculation model is shown in Figure 4. The inside of the steel shell is wrapped with PBX, and both are cubic structures; the impact loading surface is not covered by the steel shell. The PBX is subjected to a constant-speed impact along the x direction, and the base of the steel shell in the x direction is a fixed structure. The internal PBX measures 32 mm × 22 mm × 22 mm in the x, y, and z directions. The grid size of the three-dimensional PBX computing domain is set at 1 mm and contains a total of 15,488 PBX material points. Since PBX is a poly-crystalline structure, the Voronoi method is used to randomly generate seed points within the material and assign which grain a material point belongs to the distance between the material point and the seed point. As can be seen from Figure 5, PBX is composed of HMX crystals and binder under discretization modeling, so there are three different bonds inside the PBX: HMX–HMX, HMX–binder, binder–binder. A total of 30 grains are set up inside the PBX, so the average size of the grains is 516 μm. The outer layer of the PBX is wrapped in a 3 mm steel layer, leaving the impact load loading surface uncovered. The grid size of the steel layer is also 1 mm, and there are approximately 12,000 steel material points in total. During the simulation, the three layers of PD particles along the x direction of the PBX on the impact load loading surface moved at a constant speed for 50 $μ s$ , and the simulation continued until the damage evolution of the PBX stabilized. In the simulation, except for the steel layer particles at the bottom of the x direction being fixed, the rest of the steel layer and PBX particles are free particles. The time step is set to $Δ t = 2 \times 10^{- 7}$ s, which is smaller than the critical time step estimated as $Δ t_{c r} = Δ x \sqrt{ρ / E}$ using the HMX properties, and therefore satisfies the stability requirement of the explicit PD calculation. To further assess the influence of spatial discretization, a mesh-independence analysis was performed. The results show that the total-damage curves under different grid settings exhibit consistent overall evolution trends and similar stabilized values, indicating that the present discretization is sufficient for capturing the main damage evolution features. The corresponding mesh-independence results are provided in Appendix B.

Figure 4.

Schematic diagram of the cross-section of the PBX loading structure. PBX: polymer-bonded explosive.

Figure 5.

Schematic diagram of the PBX computational model and PD bonds. PD: peridynamics; PBX: polymer-bonded explosive.

PD model parameters

In PD theory, the interaction force between two nearby material points is described by paired bond force density functions. Therefore, constructing an appropriate bond force density function is crucial to better capture the true behavior of the material. The previous experimental results indicate that when PBX is under passive confining pressure, it can undergo partial plastic deformation when subjected to external force before failure, rather than breaking directly as in brittle materials. Hence, in this study, the elastoplastic constitutive model is used to capture the plastic behavior of PBX under confining pressure (Silling and Askari, 2005). The constitutive model is shown in Figure 6 and is represented as:

‖ f ‖ = {\begin{matrix} c s & if s < S_{y} \\ c S_{y} & if S_{y} < s < S_{c} \\ 0 & otherwise \end{matrix}

(23)

where

S_{y}

represents the yield relative elongation. Once the elongation of the bond between the material points is greater than the critical relative elongation, the bond fails and the bond force between these material points disappears. From the perspective of the elastoplastic constitutive model, each material has three constitutive model parameters E,

S_{y}

, and

S_{c}

to be determined. Therefore, there are a total of nine unknown PD model parameters because there are three types of materials involved in the current simulations.

Figure 6.

Schematic diagram of the relationship between PD bond force and relative elongation. PD: peridynamics.

The current research (Kostic et al., 2022) focuses on the influence of the fluctuations of these nine model parameters on the numerical simulation results and ranks their sensitivity according to the degree of influence. The literature has verified that the elastic modulus of the same batch of steel generally follows a normal distribution. However, due to the scarcity of literature on PBX in the PD field, the distribution types of other parameters are not yet clear. In this study, a uniform distribution was used to represent the uncertainty of parameters other than the elastic modulus of the steel. It should be noted that the uniform distribution is commonly used when the distribution type is unknown. The range of parameters was determined by referring to previous studies (Deng et al., 2024; Deng and Wang, 2020; Hu et al., 2024; Huang et al., 2022; Li and Hao, 2023; Prakash and Seidel, 2017). It should be noted that the sensitivity analysis results were obtained under the current parameter uncertainty. The distribution and range of the parameters are shown in Table 1.

Table 1.

Range and distribution of model parameters.

Parameters	Distribution	Upper bound/mean	Lower bound/Std
$E^{S t e e l}$	Normal	205 MPa	5 MPa
$S_{y}^{S t e e l}$	Uniform	0.002	0.003
$S_{c}^{S t e e l}$	Uniform	0.3	0.7
$E^{H M X}$	Uniform	10 MPa	32 MPa
$S_{y}^{H M X}$	Uniform	0.0002	0.0004
$S_{c}^{H M X}$	Uniform	0.001	0.005
$E^{B i n d e r}$	Uniform	3 MPa	8 MPa
$S_{y}^{B i n d e r}$	Uniform	0.004	0.006
$S_{c}^{B i n d e r}$	Uniform	0.02	0.08

Results and discussion

Dynamic damage behaviors of PBX

Before analyzing the dynamic damage behavior of PBX under impact loading, the reliability of the present PD model was first examined by comparing the simulated stress–strain curves of PBX with the corresponding experimental results. Good agreement was obtained in both the overall evolution trend and the main mechanical response characteristics, indicating that the current model and parameter settings can reasonably reproduce the basic mechanical behavior of the constituent materials. The detailed comparison results are provided in Appendix A.

To help analyze and understand the mechanism of the influence of model parameters on the dynamic damage behavior of PBX, the damage processes of PBX under impact load are analyzed first. Three different impact speeds of 20, 40, and 60 m/s are selected to simulate the dynamic damage of PBX, in order to analyze the influence of different impact speeds on PBX dynamic response and how the sensitivity of parameters affects the physical process of PBX dynamic response.

Figure 7(a) illustrates the dynamic damage evolution of PBX under the 20 m/s impact velocity, specifically depicting cross-sectional damage profiles of steel-encased PBX specimens. In Figure 7(b), the steel shell components have been removed to enhance the visualization clarity of PBX-specific damage patterns. The first two sub-figures in Figure 7(a) and (b) correspond to the impact loading phase spanning 0–40 $μ s$ , while the subsequent three sub-figures document the damage progression from load termination (50 $μ s$ ) to the stabilized material response at 100 $μ s$ . The colors of the particles in the figure, a graduated color scheme from blue to red quantifies particle-scale damage severity, where red coloration (damage value = 1) indicates complete bond rupture within discrete particles. The trans-granular damage nucleation within the HMX crystals during the initial stage of loading is first observed. This damage pattern is consistent with previous experimental observations (Xiao et al., 2017), where trans-granular fracture was observed to occur first under uniaxial loading. Such consistency suggests that the present PD model can reasonably reproduce the characteristic early-stage damage morphology of PBX.

Figure 7.

Evolution of PBX damage at impact velocity of 20 m/s: (a) a cross-sectional view of a PBX wrapped in a steel shell; and (b) a PBX model with the outer steel shell removed. PBX: polymer-bonded explosive.

Owing to superior interfacial strength between HMX grains and polymer binder compared to HMX's intrinsic fracture strength, crack propagation preferentially occurred within impact-proximal HMX crystallites, with minimal observable damage in the binder phase. As evidenced in Figure 7(a), this interfacial strength hierarchy enabled binder-mediated stress redistribution, resulting in predominant plastic deformation rather than catastrophic fracture during compressive loading. Post 40 $μ s$ , stress wave propagation drove progressive damage accumulation.

A comparative analysis of Figures 8 and 9 demonstrates conserved damage progression trends at elevated velocities (40–60 m/s), albeit with accelerated kinetics. At these higher velocities, both HMX grain damage and interfacial damage developed earlier, and the interaction between the PBX material and the surrounding steel shell became more pronounced during compression.

Figure 8.

Evolution of PBX damage at impact velocity of 40 m/s: (a) a cross-sectional view of a PBX wrapped in a steel shell; and (b) a PBX model with the outer steel shell removed. PBX: polymer-bonded explosive.

Figure 9.

Evolution of PBX damage at impact velocity of 60 m/s: (a) a cross-sectional view of a PBX wrapped in a steel shell; and (b) a PBX model with the outer steel shell removed. PBX: polymer-bonded explosive.

In our simulations, as the impact progresses, the stress wave propagates from the impacted surface toward the opposite surface of the PBX sample. Due to the poly-crystalline nature of PBX and the distinct mechanical properties of HMX and the binder, the damage evolution does not simply follow the direction of stress wave propagation. In order to clearly demonstrate the damage evolution trend of HMX grains, an HMX grain is chosen in Figure 10, and its damage process is analyzed.

Figure 10.

HMX crystal model diagram inside the PBX. PBX: polymer-bonded explosive; HMX: octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine.

Figure 11 presents the temporal evolution of damage in the HMX grain shown in Figure 10 under an impact velocity of 20 m/s. The damage field is visualized using a color gradient from blue (undamaged) to red (fully damaged), indicating the severity of damage. It can be clearly seen that damage initiates in the grain interior and progressively propagates outward the grain boundaries as the impact duration increases, eventually stabilizing as the system approaches a quasi-steady damage state.

Figure 11.

The damage processes of HMX crystal with time at an impact speed of 20 m/s. HMX: octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine.

It can be observed that HMX particles located at the core of the HMX grain are predominantly surrounded by other HMX particles, leading to a local bonding structure primarily composed of HMX–HMX interactions. Toward the grain boundary, the neighborhood of HMX particles contains an increasing number of binder particles, thereby elevating the proportion of HMX–binder bonds. Since the binder exhibits significantly better ductility than that of the HMX grain, it has a greater capacity to accommodate deformation.

To characterize the mechanical property of the HMX–binder interface, an adjustable interface parameter $β_{i n t e r}$ is introduced in the present study. The critical relative elongation of HMX–binder bonds is defined as: $s_{c}^{H M X - B i n d e r} = s_{c, min} + β_{i n t e r} (s_{c, max} - s_{c, min})$ . The $s_{c}^{H M X - B i n d e r}$ denotes the critical relative elongation of HMX–binder bonds, and $s_{c, min}$ and $s_{c, max}$ represent the minimum and maximum values of the critical relative elongations of HMX and binder, respectively.

In the present model, the adjustable interface parameter is taken as $β_{i n t e r} = 0.5$ , which assigns the yield relative elongation, critical relative elongation, and elastic modulus of HMX–binder bonds as the intermediate value between those of HMX and binder. Since the binder phase exhibits better ductility and the HMX–binder bonds are less prone to bond breakage than HMX–HMX bonds under the same local deformation, the interface region has a stronger capacity to accommodate deformation-induced damage. As a result, damage is more likely to initiate within the grain interior, while HMX particles located closer to the grain boundary, which form a larger proportion of HMX–binder bonds, tend to experience less damage at the early stage.

There are three different types of bonds in our simulation: HMX–HMX, HMX–binder, and binder–binder. The damage corresponding to different bonds is calculated by the number of broken bonds divided by its total number of initial bonds. Since the breakage of HMX–HMX bonds corresponds to the formation of cracks within the HMX grains, herein referred to as trans-granular damage, while the breakage of HMX–binder and binder–binder bonds signifies interfacial cracking, that is, inter-granular damage.

To quantitatively assess the dominant damage mechanism within PBX, Deng and Wang (2020) proposed a damage mode index (DMI), defined as:

D M I = \frac{D_{H M X - H M X}}{D_{H M X - B i n d e r} + D_{B i n d e r - B i n d e r}} = \frac{D_{C r y s t a l}}{D_{I n t e r f a c e}}

(24)

where D denotes the number of broken bonds of the corresponding type. We categorize the damage modes as follows: the number of broken HMX–HMX bonds is designated as

D_{C r y s t a l}

to characterize trans-granular damage, while the number of broken HMX–binder and binder–binder bonds is collectively denoted as

D_{I n t e r f a c e}

to characterize inter-granular damage. This classification is based on the physical locations where different bond failures occur in the PBX micro-structure. Since HMX–HMX bonds connect material points within the same crystal, their breakage directly represents crack initiation and propagation inside HMX grains, and can therefore be used to describe trans-granular damage. In contrast, HMX–binder bonds and binder–binder bonds are mainly distributed near grain boundaries and within the binder-related regions. Their breakage reflects interfacial debonding, binder failure, and crack growth along or across the crystal–binder boundaries, which are characteristic features of inter-granular damage. Therefore,

D_{C r y s t a l}

and

D_{I n t e r f a c e}

are not empirical variables introduced only for convenience, but physically meaningful quantities directly associated with two distinct damage modes in PBX. A value of

D M I > 1

indicates that trans-granular damage predominates, whereas

0 < D M I < 1

implies that inter-granular damage is dominant.

To quantify the variability of damage modes induced by parameter uncertainty, N samples were first drawn from the prescribed distributions of the uncertain PD model parameters after the Kriging surrogate model was established. The corresponding DMI values were then predicted by the surrogate model as ${\hat{D M I}}^{(i)} = g_{k} (X^{(i)}) (i = 1, 2, \dots, N)$ , where $X^{(i)}$ denotes the $i$ -th sampled parameter vector and $g_{k}$ represents the trained Kriging model. Based on the resulting DMI sample set, the mean DMI was calculated as $\bar{D M I} = (1 / N) \sum_{i = 1}^{N} {\hat{D M I}}^{(i)}$ , and the 95% interval was estimated using the 2.5th and 97.5th percentiles, that is, $[Q_{0.025} (\hat{D M I}), Q_{0.975} (\hat{D M I})]$ . Figure 12 illustrates the mean values of DMI under three impact velocities (20, 40, and 60 m/s), together with the corresponding percentile-based 95% intervals. The error bars therefore represent the potential variability in the dominant damage mode arising from the uncertainty of material parameters.

Figure 12.

The variation range of the damage mode index at impact velocities of 20, 40 and 60 m/s (the error bars are expressed as 95% confidence intervals).

At an impact velocity of 20 m/s, the mean DMI is approximately 10.257, indicating that trans-granular damage is the dominant failure mode under low-velocity impact. Notably, the error bar at 20 m/s is much longer than those at 40 and 60 m/s, and is also strongly asymmetric, with a particularly large upper bound. This feature indicates that the DMI distribution at low velocity is highly dispersed and right-skewed. Physically, under 20 m/s impact, trans-granular damage is dominant, while inter-granular damage remains relatively limited. Since DMI is defined as the ratio of $D_{C r y s t a l}$ to $D_{I n t e r f a c e}$ , a relatively small $D_{I n t e r f a c e}$ makes the DMI highly sensitive to parameter-induced fluctuations, so even moderate changes in material parameters can produce a large variation in the ratio. As a result, the uncertainty of the dominant damage mode is amplified at low velocity, which is reflected by the very long upper error bar.

As the impact velocity increases to 40 m/s, the mean DMI decreases to 0.96052, which is close to 1. This indicates that the system enters a transition regime between trans-granular and inter-granular damage. Correspondingly, the error bar becomes much shorter than that at 20 m/s, but it is still wider than that at 60 m/s and extends across the critical value of $D M I = 1$ . This suggests that, under 40 m/s impact, the dominant damage mode is no longer strongly biased toward one side, and uncertainty in model parameters can still shift the damage response between the two competing modes. Therefore, 40 m/s can be regarded as a transitional velocity at which the competition between trans-granular and inter-granular damage is strongest.

At 60 m/s, the mean DMI further decreases to 0.81462, and the corresponding error bar is the narrowest among the three cases. More importantly, the entire interval remains below 1, indicating that inter-granular damage becomes the stably dominant failure mode at high impact velocity. This result suggests that, although parameter uncertainty still affects the quantitative value of DMI, it no longer changes the qualitative judgment of the dominant damage mode as easily as in the lower-velocity cases.

Kriging surrogate model

Sensitivity analysis requires comprehensive sampling to capture variations in parameter space and their impacts on model outputs when handling complex systems or non-linear models, as parameters often exhibit multi-dimensional, non-linear, or interactive effects, making sparse sampling insufficient to cover all potential parameter combinations. Therefore, a large number of samples are needed in sensitivity analysis to ensure a uniform distribution in the parameter space. However, practical numerical models may involve high-computational costs and prolonged simulation times, which conflict with the demand for extensive repeated calculations. According to the Section “Global Kriging surrogate model,” Kriging surrogate models corresponding to impact loading conditions at 20, 40, and 60 m/s can be established to address these computational challenges. The critical step in constructing surrogate models lies in defining inputs and outputs: uncertain parameters serve as inputs, while outputs correspond to the model's response to input variations, whose formulation is intrinsically tied to the analysis objectives.

This study focuses on analyzing the influence of constitutive model parameter uncertainties on the dynamic damage outcomes of PBX under impact loading. Since damage is quantified by bond crack counts, the constitutive models of PBX and steel are designated as inputs, with outputs being the crack counts of HMX–HMX, HMX–binder, and binder–binder bonds. As shown in Figure 13, under sustained impact loading, PBX eventually undergoes severe damage and fragmentation, leading to damage values approaching 0.8 across different velocities. The limited variability in output uncertainties results in nearly identical parameter importance rankings at various velocities, failing to adequately represent velocity-dependent parameter sensitivities.

Figure 13.

Damage of PBX at different impact speeds (solid line shows the damage caused by the continuous impact on the PBX particles, dashed line shows the damage caused by the impact applied for a time of 50 $μ s$ ). PBX: polymer-bonded explosive.

Under moderate strain-rate conditions, Pei et al. (2024) investigated trans-granular fracture and interfacial debonding in pressed explosives under multi-axial non-shock loading by adjusting incident wave durations to 50 $μ s$ . Similarly, a 50 $μ s$ impact duration at different velocities induces varying damage levels in PBX while allowing sufficient time for damage evolution to stabilize, thereby effectively capturing velocity-specific influences on PBX damage progression.

Figure 14 displays the relative damage variation rate of PBX over time under fixed 50 $μ s$ impact durations at different velocities. The damage value corresponding to a relative variation rate below $0.1 %$ is selected as the stabilized post-evolution damage of PBX, which serves as a rational output for sensitivity analysis. Consequently, this study defines the dynamic response damage of PBX after 50 $μ s$ impact duration as the output, with PD constitutive model parameters as inputs. Kriging surrogate models for the three impact velocities are trained using a stopping criterion of root mean square error $\leq 0.5 %$ . The training sample sizes and convergence errors are summarized in Table 2. Because an error-controlled adaptive sampling strategy was used, the training sample size was not fixed in advance but depended on the complexity of the response surface under each impact velocity. The 20 m/s case required more training samples than the 40 and 60 m/s cases, indicating a more complex relationship between parameter variation and damage response at the lower impact velocity. In addition, the predictive accuracy of the surrogate models was further evaluated using independent validation samples. The $R^{2}$ values were 0.9946, 0.9514, and 0.9265 for the three impact velocities, respectively, indicating that the Kriging models can accurately reproduce the PD simulation results.

Figure 14.

Relative rate of change of damage during continuous evolution. The red area represents the range of relative rate of change of 0.1% that meets the convergence condition.

Table 2.

Training samples and mean square errors of the Kriging surrogate model at three impact speeds.

Surrogate model (m/s)	Sample size	MSE
20	700	0.0008
40	430	0.0026
60	480	0.0024

MSE: mean square error.

The sensitivity analysis of model parameters

The total damage of PBX is defined as the cumulative number of bond breakages among HMX–HMX, HMX–binder, and binder–binder. This metric reflects the global accumulation of internal defects and their subsequent influence on the mechanical integrity of the material. By evaluating the sensitivity of total damage to various material parameters, one can identify the dominant sources of uncertainty and, consequently, determine the key factors governing the damage evolution in PBX composites.

Figure 15 illustrates the first-order and total Sobol sensitivity indices of the PD constitutive model parameters with respect to total damage under impact velocities of 20, 40, and 60 m/s. At an impact velocity of 20 m/s, the most influential parameter is the critical relative elongation of the binder ( $S_{c}^{B i n d e r}$ ), followed by the elastic modulus of HMX ( $E^{H M X}$ ). Together, the uncertainties of these two parameters (as quantified by their variances) account for approximately 95% of the variance in total damage, indicating their critical role in damage prediction under this loading condition. It should be noted that the uncertain model parameters may interact with each other during the numerical simulation. In the present study, such interaction effects are implicitly taken into account through the total Sobol indices, which include not only the main effect of each parameter but also its coupled effects with all other parameters.

Figure 15.

First-order Sobol index and total index calculation of PD model parameters for total damage under different impact velocities: (a) 20 m/s, (b) 40 m/s, and (c) 60 m/s. PD: peridynamics.

At 20 m/s, trans-granular damage dominates the damage mechanism. Given the relatively limited damage experienced by the binder at this stage, its material parameters exhibit greater regulatory influence over the total damage than those of HMX, and thus emerge as the most significant contributors. As the impact velocity increases, damage initiates within the grains and progressively propagates toward the grain boundaries. Concurrently, the dominant damage mode converts from trans-granular to inter-granular fracture. This shift in failure mechanism is reflected in the changing importance of material parameters: the increasing damage to the binder reduces its ability to regulate total damage, leading to a corresponding decrease in the significance of its associated parameters.

At impact velocities of 40 and 60 m/s, the Sobol index of $E^{H M X}$ surpasses that of the $S_{c}^{B i n d e r}$ , becoming the most influential factor on total damage. Moreover, its contribution continues to rise with increasing velocity. Nevertheless, the combined influence of $E^{H M X}$ and $S_{c}^{B i n d e r}$ still accounts for over 90% of the total damage uncertainty, suggesting that while the dominant failure mechanisms may shift, the overall trend of damage evolution remains largely consistent. Therefore, the changes in total damage response are closely related to the variation in dominant fracture modes at different impact velocities.

The velocity-dependent differences in parameter sensitivity also provide an explanation for the different error-bar lengths observed among the three impact velocities in Figure 12. At 20 m/s, parameter perturbations can induce relatively large fluctuations in $D_{I n t e r f a c e}$ , which amplifies the variability of the ratio $D_{C r y s t a l} / D_{I n t e r f a c e}$ and leads to the widest DMI interval. As the impact velocity increases, interfacial damage develops more sufficiently and the dominant damage mode becomes more stable, so the influence of parameter uncertainty on DMI is reduced, resulting in narrower error bars.

In the simulation, the total damage of PBX is defined as the cumulative number of broken bonds among three types. Each bond type corresponds to damage occurring in different regions of the micro-structure. Due to the poly-crystalline nature of PBX and the ability of HMX and binder to form interfacial bonds, the breakage of a particular bond type is not influenced solely by the material parameters of that specific phase, but may also be affected by the properties of other materials within the composite. Therefore, analyzing the sensitivity of different bond breakages to various material parameters not only identifies key parameters that dominate bond failure, but also sheds light on the inter-dependence among different material phases during damage evolution. In PBX, the initial population of each bond is unequal, and due to their distinct evolution patterns under impact loading, the Sobol sensitivity analysis of the three bonds under related conditions is considered to decompose the uncertainty contributions of each bond to the total damage.

The inner pie chart of Figure 16 depicts the proportion of variance in $D_{T o t a l}$ attributed to the uncertainties in $D_{H M X - H M X}$ , $D_{H M X - B i n d e r}$ , and $D_{B i n d e r - B i n d e r}$ . The outer pie chart distinguishes between the independent contributions (orange) and the correlated contributions (gray) from each bond type to the overall uncertainty. Figure 17 presents the first-order Sobol indices and total Sobol indices of all material parameters for the two damage modes, with the two columns corresponding to $D_{C r y s t a l}$ and $D_{I n t e r f a c e}$ , respectively. An important observation is that the independent contributions of the three bond types to total damage uncertainty sum to only about 5%, indicating that total damage is not governed by any single bond type in isolation, but is controlled predominantly by the coupled evolution among all three bond types. This result reflects the highly interactive nature of damage development in PBX. Because PBX is a heterogeneous composite composed of HMX grains, binder, and their interfaces, the failure of one bond type inevitably alters the local stress redistribution and deformation compatibility of the surrounding micro-structure. For example, the breakage of HMX–HMX bonds reduces the structural integrity of the crystal interior and promotes stress transfer toward grain boundaries and binder-related regions, thereby facilitating the failure of HMX–binder and binder–binder bonds. Conversely, interfacial debonding and binder damage weaken the mechanical constraint around HMX grains, which can enhance local stress concentration inside the grains and further promote trans-granular cracking. Therefore, the three damage modes do not evolve independently, but mutually influence and reinforce each other during impact loading.

Figure 16.

The individual and total contributions of the uncertainty of the three types of bond damage to the total damage uncertainty at different impact speeds.

Figure 17.

First-order Sobol indicators and total indicators of all parameters for the two damage modes at impact speeds, v = 20, 40, and 60 m/s.

The very small sum of independent contributions further suggests that the uncertainty in total damage arises mainly from the correlated fluctuations among bond types rather than from the isolated variation of any one damage component. In other words, total damage is primarily a collective response of the PBX micro-structure to impact, rather than a simple linear superposition of three separate failure processes. This also explains why adjusting the broken-bond number of only one bond type is generally insufficient to significantly change the overall damage evolution, whereas the coupled variation among crystal fracture, interfacial debonding, and binder failure plays a dominant role in determining the final damage state.

At an impact velocity of 20 m/s, the pie chart in panel Figure 16(a) clearly shows that approximately 80.6% of the total damage uncertainty originates from $D_{I n t e r f a c e}$ (sum of the Sobol indices of HMX–binder and binder–binder), while only 19.4% is attributed to $D_{C r y s t a l}$ (Sobol index of HMX–HMX). Despite trans-granular damage being the dominant failure mode at this velocity, the inherently brittle nature of the grain interiors leads to high deformation between HMX particles, greatly exceeding the critical failure strain of HMX–HMX bonds. As a result, nearly all HMX–HMX bonds within the grains undergo fracture following elastic deformation. This explains why the HMX–HMX bond has the highest fracture rate among all types but exhibits the least uncertainty.

In contrast, under lower impact loads, damage at the interfaces is comparatively minimal, with a greater proportion of bonds in a critical (near-failure) state. This makes the evolution of $D_{I n t e r f a c e}$ more sensitive to changes in material parameters, which aligns the Sobol indices for $D_{I n t e r f a c e}$ with those for total damage.

From the first-order Sobol indices of $D_{C r y s t a l}$ in Figure 17(a), it is evident that the most influential parameters are HMX's constitutive properties: 49.6% of the variance originates from the fracture strain of HMX and 33.6% from its elastic modulus. In Figure 17(b), 86.7% of the influence comes from the binder's PD model parameters for $D_{I n t e r f a c e}$ , with 33.6% specifically attributed to $E^{H M X}$ . This is due to the interface factor being set to 0.5, meaning that the modulus uncertainty of HMX–binder bonds is primarily governed by $E^{H M X}$ , while their critical relative elongation uncertainty is dominated by $S_{c}^{B i n d e r}$ . Moreover, the number of HMX–binder bonds is approximately equal to that of binder–binder bonds in PBX, making the contribution of $E^{H M X}$ significantly greater than that of $E^{B i n d e r}$ in $D_{I n t e r f a c e}$ .

As shown in Figure 11, damage in HMX grains initiates from the interior and propagates outward. Under lower impact velocities, damage is concentrated within the HMX grains, where most HMX–HMX bonds are completely broken. Closer to the grain boundary, binder particles become more prevalent, resulting in increased bonding between HMX and binder. Since HMX–binder bonds are more tolerant to deformation-induced failure under the current model setting, edge HMX particles tend to experience less damage, leaving many of their neighboring HMX–HMX bonds in a critical state. Furthermore, the binder undergoes less damage than HMX, and a substantial number of bonds at the interface remain near the failure threshold. This explains why, in the importance analysis of both $D_{C r y s t a l}$ and $D_{I n t e r f a c e}$ , when most bonds are in a critical state, $S_{c}^{H M X}$ and $S_{c}^{B i n d e r}$ exert greater influence than $E^{H M X}$ on their respective bond types.

As the impact velocity increases, the dominant damage mechanism in PBX shifts from trans-granular to inter-granular. Correspondingly, the influence of $D_{I n t e r f a c e}$ on total damage decreases, leading to a divergence between the parameter importance rankings for total damage and $D_{I n t e r f a c e}$ . At 40 m/s, the importance of $S_{c}^{H M X}$ for $D_{C r y s t a l}$ drops below 10%, far less than that of $E^{H M X}$ . Although the ranking of dominant parameters for $D_{I n t e r f a c e}$ remains unchanged, the importance of $S_{c}^{B i n d e r}$ similarly decreases, following the same trend. This suggests that, under higher impact velocities and deeper material damage, particle displacements greatly exceed the critical strain thresholds of bonded pairs, thereby reducing the sensitivity of $S_{c}$ and relatively increasing the importance of E. As PBX experiences more extensive damage, the effects of other material parameters on both damage modes also become more pronounced. For example, the influence of $E^{B i n d e r}$ on $D_{C r y s t a l}$ increases from 3.7% to 8.3%, and the yield relative elongation of steel $S_{y}^{S t e e l}$ shows increased influence on both damage modes.

In summary, the shift in dominant damage mechanisms in PBX leads to the corresponding changes in the ranking of parameter importance. When the damage is relatively mild, $S_{c}$ plays a dominant role. However, as the damage becomes more severe, its influence gradually decreases, whereas the effect of E becomes increasingly pronounced. In addition, under large deformation, the interactions among different material phases are strengthened, which amplifies the contributions of other material parameters to both $D_{C r y s t a l}$ and $D_{I n t e r f a c e}$ . These correlation effects among sensitive parameters also help explain the different widths of the DMI intervals. At low impact velocity, the coupled perturbations of key parameters can induce relatively large fluctuations in both $D_{C r y s t a l}$ and $D_{I n t e r f a c e}$ , making the ratio $D_{C r y s t a l} / D_{I n t e r f a c e}$ more variable. As the impact velocity increases, inter-granular damage becomes more fully developed and the dominant damage mode becomes more stable, thereby weakening the influence of parameter correlations on the DMI interval.

Conclusion

By developing a PD model for steel-encased PBX, this study successfully simulated the dynamic response of PBX under impact loading with confinement effects. Utilizing PD constitutive model parameters as stochastic inputs and three types of broken bond number (HMX–HMX, HMX–binder, binder–binder) as damage outputs, the effects of constitutive parameter uncertainties on dynamic damage of PBX are quantified. Kriging surrogate models were constructed for 20, 40, and 60 m/s impact velocities using these datasets. At a relatively low impact velocity (20 m/s), trans-granular damage dominates within the PBX micro-structure and exhibits pronounced sensitivity to uncertainties in the PD model parameters. This is evidenced by a significantly high DMI of >1 accompanied by a wide confidence interval, indicating substantial variability in damage evolution induced by parameter fluctuations. As the impact velocity increases to 40 and 60 m/s, the dominant failure mechanism transitions to inter-granular damage, with the DMI markedly decreasing to values below 1. Concurrently, the influence of material parameter uncertainty on the damage response becomes less pronounced, as reflected by the narrowing of the DMI confidence intervals.

Sobol global sensitivity analysis was implemented to compute the importance indices of PD model parameters, identify critical factors influencing damage progression, and establish a systematic framework for parameter sensitivity evaluation. The parameter sensitivity analysis based on the simulation results reveals a clear velocity-dependent dominance: at low impact velocity (20 m/s), $S_{c}^{B i n d e r}$ dominates the uncertainty in total damage, with a contribution of 71.8%, highlighting its pivotal role in regulating damage evolution through the deformation accommodation capacity of the binder and the interface region. This dominance can be physically attributed to the relatively mild stress-wave amplitude at low velocity, under which damage initiation and growth remain strongly affected by the failure threshold of the binder phase. In contrast, under high-velocity impact (60 m/s), the sensitivity of $E^{H M X}$ increases to 56.3%, reflecting the reduced protective effect of binder deformation and the increasing role of HMX stiffness in controlling stress concentration and load transfer. From the perspective of stress-wave propagation, the wave speed can be approximately estimated as $c = \sqrt{E / ρ}$ . For HMX and binder, the estimated wave speeds are on the order of several kilometers per second, and the stress wave can traverse the PBX domain several times during the 50 μs loading duration. Therefore, as the impact velocity increases, repeated wave propagation, reflection, and local stress redistribution become more intense, making the stiffness contrast between HMX and binder increasingly important in governing damage evolution.

Moreover, the sensitivity analysis of damage modes indicates that, under 20 m/s impact conditions, the relatively mild damage causes the uncertainties in trans-granular and inter-granular damage to depend mainly on the variations in HMX-related and binder-related parameters, respectively. However, under 40–60 m/s impact conditions, $E^{H M X}$ exerts a critical influence on both trans-granular and inter-granular fracture processes, showing a pronounced coupling effect. This suggests that the dominant control mechanism gradually shifts from binder-governed damage accommodation at low velocity to HMX stiffness-governed stress transfer and fracture development at higher velocities. Furthermore, the sensitivity analysis indicates that uncertainty in total damage primarily originates from the coupled interactions among the three bond types rather than from any individual bond. This emphasizes that the overall damage response of PBX is not a simple superposition of isolated bond failures, but a synergistic outcome of crystal fracture, interfacial debonding, and binder failure within the heterogeneous micro-structure.

Footnotes

ORCID iDs

Yukai Cui

Xiaoliang Deng

Luyi Li

Author contributions

Yukai Cui: conceptualization, methodology, and writing—validation, investigation, and original draft. Xiaoliang Deng: conceptualization, methodology, writing—review and editing, validation, and investigation. Luyi Li: methodology and writing—review, editing, and validation.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors gratefully acknowledge financial support from the National Natural Science Foundation of China (through Grant No. 12272361) and the scientific research project of National Key Laboratory of Shock Wave and Detonation Physics (Grant No. JCKYS2024212104).

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Data availability statement

The data used in this study will be made available upon request to the corresponding author.

Appendix A

To examine the reliability of the present PD model, the simulated stress–strain response was compared with the corresponding experimental result. The same PD constitutive model and parameter settings as those used in the main simulations were adopted.

As shown in Figure A1, the simulated stress–strain curve agrees well with the experimental curve. The PD model captures the main mechanical response of PBX, including the initial loading stage, the peak stress, and the rapid post-peak stress drop caused by damage and fracture. The predicted peak stress and the corresponding strain are also close to the experimental values.

Although slight differences can be observed in the post-peak stage, this is reasonable because the failure process of PBX is affected by local micro-structural heterogeneity, grain distribution, binder deformation, and interfacial debonding. Overall, the comparison indicates that the present PD model can reasonably reproduce the basic mechanical response of PBX, supporting its use for subsequent dynamic damage and sensitivity analyses.

Appendix B

To evaluate the influence of spatial discretization on the simulated damage evolution, a mesh-independence analysis was conducted. Three models with different mesh densities under the 20 m/s conditions were considered, including the baseline mesh used in the main simulations, a refined mesh with twice the number of grid points, and a further refined mesh with five times the number of grid points. Except for the mesh density, the model geometry, material parameters, boundary conditions, and impact-loading conditions were kept the same.

Figure B1 compares the total damage evolution obtained from the three mesh settings. The total-damage curves show consistent overall trends under different mesh densities. In all cases, total damage increases rapidly after impact loading and then gradually approaches a stabilized state. The results obtained from the baseline mesh and the 2× mesh are very close. Although the 5× mesh produces slight differences in the absolute damage value, the main evolution pattern, stage characteristics, and final stabilized state remain consistent with those of the baseline model.

These results indicate that the spatial discretization used in the main simulations is sufficient to capture the primary damage-evolution behavior of PBX. Therefore, the present mesh setting does not change the main conclusions regarding damage evolution, damage-mode transition, or parameter effect analysis.

References

Ahmadi

Sadighi

Hosseini-Toudeshky

(2022) Microstructure-based deformation and fracture modeling of particulate reinforced composites with ordinary state-based peridynamic theory. Composite Structures 279: 114734.

Barua

Horie

Zhou

(2012a) Energy localization in HMX-Estane polymer-bonded explosives during impact loading. Journal of Applied Physics 111(5): 054902.

Barua

Horie

Zhou

(2012b) Microstructural level response of HMX-Estane polymer-bonded explosive under effects of transient stress waves. Proceedings of the Royal Society A-Mathematical Physical and Engineering Sciences 468(2147): 3725–3744.

Barua

Kim

Horie

, et al. (2013) Ignition criterion for heterogeneous energetic materials based on hotspot size-temperature threshold. Journal of Applied Physics 113(6).

Bennett

Haberman

Johnson

, et al. (1998) A constitutive model for the non-shock ignition and mechanical response of high explosives. Journal of the Mechanics and Physics of Solids 46(12): 2303–2322.

Deng

Huang

Zhao

(2024) Numerical investigation of damage and ignition behaviors of PBX under punch loading. Engineering Fracture Mechanics 298: 109903.

Deng

Wang

(2020) Peridynamic modeling of dynamic damage of polymer bonded explosive. Computational Materials Science 173: 109405.

Denimal

Sinou

(2021) Advanced Kriging-based surrogate modelling and sensitivity analysis for rotordynamics with uncertainties. European Journal of Mechanics-A/Solids 90: 104331.

Dienes

Zuo

Kershner

(2006) Impact initiation of explosives and propellants via statistical crack mechanics. Journal of the Mechanics and Physics of Solids 54(10): 2235–2240.

10.

Diyaroglu

Oterkus

Madenci

, et al. (2016) Peridynamic modeling of composite laminates under explosive loading. Composite Structures 144: 14–23.

11.

Field

(1992) Hot spot ignition mechanisms for explosives. Accounts of Chemical Research 25(11): 489–496.

12.

Franzelin

Diehl

Pflüger

(2014)Non-intrusive Uncertainty Quantification with Sparse Grids for Multivariate Peridynamic Simulations. Switzerland: Springer International Publishing.

13.

Gruau

Picart

Belmas

, et al. (2009) Ignition of a confined high explosive under low velocity impact. International Journal of Impact Engineering 36(4): 537–550.

14.

Bobaru

(2011) Characteristics of dynamic brittle fracture captured with peridynamics. Engineering Fracture Mechanics 78(6): 1156–1168.

15.

Deng

Liu

(2024) Anisotropic peridynamic simulation of dynamic response of PBX containing polycrystalline HMX under low velocity impact. Engineering Analysis with Boundary Elements 166: 105861.

16.

Huang

Tao

Wei

, et al. (2024) Microstructure effect analysis of carbon black-filled rubber composites. Métodos Numéricos Para Cálculo y Diseño en Ingeniería: Revista Internacional 40(1): 1–15.

17.

Huang

Deng

Bai

(2022) Peridynamic investigation of dynamic damage behaviors of PBX confined in spherical steel shells. Mechanics of Materials 72: 104389.

18.

Jiang

Chen

German

(2016) Multidisciplinary statistical sensitivity analysis considering both aleatory and epistemic uncertainties. AIAA Journal 54(4): 1326–1338.

19.

Kim

Barua

Horie

, et al. (2014) Ignition probability of polymer-bonded explosives accounting for multiple sources of material stochasticity. Journal of Applied Physics 115(17).

20.

Kim

Wei

Horie

, et al. (2018) Prediction of shock initiation thresholds and ignition probability of polymer-bonded explosives using mesoscale simulations. Journal of the Mechanics and Physics of Solids 114: 97–116.

21.

Kostic

Miljojkovic

Simunovic

, et al. (2022) Uncertainty in the determination of elastic modulus by tensile testing. Engineering Science and Technology-An International Journal-Jestech 25: 100998.

22.

Miao

, et al. (2025) A fully coupled mechano-chemical digital model for environment-induced damage and damage propagation. Métodos Numéricos Para Cálculo y Diseño en Ingeniería: Revista Internacional 41(2): 1–21.

23.

Hao

(2023) On deformation and fracture of PBX simulant employing modified three-body potential peridynamics model with deformation-based failure criteria. Applied Mathematical Modelling 115: 100–126.

24.

Liu

Cai

Ong

(2017) An adaptive sampling approach for Kriging metamodeling by maximizing expected prediction error. Computers & Chemical Engineering 106: 171–182.

25.

Liu

Kong

Deng

, et al. (2025) Prediction and uncertainty quantification of the fatigue life of corroded cable steel wires using a Bayesian physics-informed neural network. Journal of Bridge Engineering 30(5): 04025018.

26.

Madenci

Oterkus

(2013) Peridynamic Theory and Its Applications. Springer.

27.

Napadensky

Stresau

Savitt

(1960) The behavior of explosives at impulsively induced high rates of strain.Red Hook, New York: Proceedings of the 3rd international detonation symposium.

28.

Pei

Xiao

, et al. (2024) Mesoscopic damage model and constitutive model of explosive transmission under impact load. Journal of Ordnance Equipment Engineering 45(2): 150–157.

29.

Prakash

Seidel

(2017) Computational electromechanical peridynamics modeling of strain and damage sensing in nanocomposite bonded explosive materials (NCBX). Engineering Fracture Mechanics 177: 180–202.

30.

Prakash

Seidel

(2018) Effects of microscale damage evolution on piezoresistive sensing in nanocomposite bonded explosives under dynamic loading via electromechanical peridynamics. Modelling and Simulation in Materials Science and Engineering 26(1): 015003.

31.

Silling

(2000) Reformulation of elasticity theory for discontinuities and long-range forces. Journal of the Mechanics and Physics of Solids 48(1): 175–209.

32.

Silling

Askari

(2005) A meshfree method based on the peridynamic model of solid mechanics. Computers & Structures 83(17–18): 1526–1535.

33.

Silling

Epton

Weckner

, et al. (2007) Peridynamic states and constitutive modeling. Journal of Elasticity 88(2): 151–184.

34.

Sobol

(2001) Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Mathematics and Computers in Simulation 55(1–3): 271–280.

35.

Sundararajan

Keerthi

(2001) Predictive approaches for choosing hyperparameters in Gaussian processes. Neural Computation 13(5): 1103–1118.

36.

Wyart

Coulon

Pardoen

, et al. (2009) Application of the substructured finite element/extended finite element method (S-FE/XFE) to the analysis of cracks in aircraft thin walled structures. Engineering Fracture Mechanics 76(1): 44–58.

37.

Xiao

Pei

Zou

, et al. (2024) Dynamic damage and non-shock ignition behavior of polymer bonded explosive under lower velocity impact. Propellants, Explosives, Pyrotechnics 49(7): e202300016.

38.

Xiao

Sun

Zhen

, et al. (2017) Characterization, modeling and simulation of the impact damage for polymer bonded explosives. International Journal of Impact Engineering 103: 149–158.

39.

Gertner

(2008) Uncertainty and sensitivity analysis for models with correlated parameters. Reliability Engineering & System Safety 93(10): 1563–1573.

40.

Zhang

Bobaru

(2016) Modeling the evolution of fatigue failure with peridynamics. The Romanian Journal of Technical Sciences Applied Mechanics 61(1): 22–40.