Failure Modeling and Sensitivity Analysis of Ceramics Under Impact

Ceramics are brittle or quasi-brittle materials used in a vast range of applications from body and vehicle armors, semi-conductors, scratch-resistant shields, kitchenware, and everyday appliances. While the properties of ceramics can be as varied as their chemical composition and structure, armor ceramics in particular are characterized by high energy absorption, high hardness, high compressive strength and preferably low density. All these properties make them suitable for impact and blast resistance during combat. Boron carbide, with its relatively low density and high hardness, is of special interest for research in the field of armor ceramics [1].

A ceramic constitutive model suitable for a given size and/or application might not be able to capture the relevant physics for another application. Particle-based models [2] attempting to capture the very large number of microscale defects in ceramics are computationally intensive. Numerical models of discrete cracks [3,4] or crack-like features [5–10] are not suitable for modeling simultaneous propagation of millions of micro-cracks, as is often the case for high rate simulations of armor ceramics. Continuum damage models [11–14] are applicable so long as the continuum assumption holds true. The applicability of such models in penetration or fragmentation problems are dependent on the numerical solver and the output property of interest. Ensemble averaging of atomistic properties of single crystal ceramics provides a more meaningful representation of micro-mechanical properties such as fracture toughness, moduli, and Poisson’s ratio. However, these models also have some non-physical parameters that are difficult to calibrate experimentally. Often, calibration of these parameters based on macroscopic response might be appropriate for modeling the physics and mechanisms involved in the experiments they are used to calibrate against but may not be suitable for other loading scenarios. Despite these limitations, many of these models can be used for extreme environment simulations, some of which are often infeasible or expensive to replicate in laboratories. They can also help us understand key trends and guide material processing toward improving material performance.

Like most other materials, scale separation becomes quintessential in modeling armor ceramics. From the point of simulations, this might not only mean differences in the physics at different scales but also the calibration of parameters. As an example, fragmentation in ceramics is observed typically at two scales: a microstructure-dependent scale leading to smaller fragments and a geometry-dependent macroscale [15]. Depending on the problem, fragmentation can either lead to granular phase transition and granular flow as observed in the Mescall region [16,17], or it can also lead to disintegration and limit the peak strength, as observed in unconfined Kolsky bar experiments [18]. From a modeling viewpoint, this can lead to a different calibration of the fragmentation or granular transition criterion for different experiments, depending on the dominant fragmentation mechanism. In addition to this, the resolution of the simulations limits the fragment size that can be captured in a continuum granular mechanics model. The fragment statistics are used to calibrate such granular mechanics models [19] and thereby influence the calibration at a given resolution.

The key micro-mechanisms in boron carbide under impact are amorphization [20], crack growth, crack interaction, and crack coalescence [21] leading to fragmentation followed by granular flow [22]. The current work exercises an integrated ceramics model [23–25] that combines these key features with a modified granular transition criterion developed by Ref. [19] and tries to simulate key failure mechanisms observed in Sphere Indentation experiments [26] and Edge On Impact experiments [27–29].

In this study, the integrated model has been used to assess the sensitivity of material behavior to selected model parameters for sphere indentation simulations. Broadly, there are two types of sensitivity analysis methods in the existing literature [30]: local and global. Some of the popular local methods include the one-at-a-time (OAT) method and the derivative-based methods. OAT method [31–33] involves studying the effect of the output by locally perturbing one input variable at a time while keeping the other variables at their baseline values. Local derivative-based methods [34–36] involve estimating the partial derivative of the output with respect to an input variable by small perturbations of that variable around a fixed point in the input space. The partial derivatives act as natural sensitivity measures. Although fast and transparent, the main drawback of local methods is that they do not thoroughly explore the entire input space. The global sensitivity methods try to alleviate this problem and also try to study the effect of large changes in the input variables. Among the global sensitivity analysis method, variance-based methods [31,37,38] and linear regression analysis [39,40] have been most extensively used. Variance-based methods are based on a decomposition of the variance of the output into terms corresponding to the different input variables as well as their interactions. The sensitivity measure of the output for a particular input variable is the amount of variance in the output attributable to that input variable. Since variance-based methods allow the examination of sensitivities across the whole input space, they often make use of emulators/surrogates [41–44] to reduce the computational expense of too many model runs. Linear regression involves fitting of the input-output data assuming that the output is linearly related to the data. The standardized regression coefficients obtained from the regression fit can then be used as measures of parameter sensitivity.

In this work, sensitivity analysis is performed based on linear regression fit to the ceramics model data to study the importance of parameters with respect to the selected quantities of interest. The most sensitive parameters have been identified, and their trends have been introspectively speculated and compared against simulation results. Finally, suggestions have been made toward material processing modifications and prioritization of calibration experiments and/or more in-depth modeling to support-specific parameters.

The integrated model builds on Ref. [23] and combines multiple mechanisms that are deemed dominant in boron carbide. It incorporates amorphization induced damage from Ref. [45], fracture dominated by the growth and interaction of wing cracks from sliding flaws [23,46], and fragmentation and transition to granular mechanics [19] and granular plasticity. A Continuum breakage mechanics model (CBM) [47,48] is used to simulate granular flow-induced plasticity. The current implementation in the integrative model is an improvement over previous implementations by Refs. [23–25] as it combines all the mechanisms independently developed and/or improved by the aforementioned authors with a modified physics-based granular transition criterion based on Ref. [19] with the recalibration of some model parameters.

K₀, G₀, and ν₀ are the bulk modulus, shear modulus, and the Poisson’s ratio of the undamaged material.

Calibration of mechanical, equation of state, microstructural, micromechanical, amorphization and granular flow parameters follow [24,25]. Critical transition damage for initiation of granular mechanics [23–25] is newly recalibrated in this work as per [19]. Critical breakage energy density (E_C) in Ref. [47] is recalibrated using Ref. [56]. Grading index (⁠$ϑ$⁠) is re-estimated and some of the arguments of fragment size evolution for granular media in Ref. [47] are revisited in the context of a low porosity fragmenting ceramic. These updated calibrations are discussed in the sections that follow. A summary of the calibrated integrative model parameters is presented in Table 1.

For the chosen material properties (Table 1), $Ωfd=3.5$ and $Ωfr=2.9$⁠.

Reference [56] explores the particle size effect in the behavior of granular boron carbide under quasi-static compression. The evolution of porosity, bulk modulus with changing hydrostatic pressure is used to compute the critical breakage energy density, E_C [47,55]. Equation (5.7) in Ref. [55] is used to compute the critical breakage energy density. The grading index (⁠$ϑBM$⁠) is calculated from the initial and final fragment size distribution for granular boron carbide from Ref. [56]. The pressure beyond which the bulk modulus of granular boron carbide starts increasing can be thought of as the critical comminution pressure, p_CR (Eq. (5.7) in Ref. [55]) at which the breakage of fragments initiates. The grading index, critical comminution pressure, and the corresponding bulk modulus (K_g) are used to calculate the value of E_C. The results of the calculation are listed in Table 2.

In general, it can be expected that the strength of an individual fragment might decrease with an increase in fragment size due to the higher statistical probability of defects. However, the critical breakage energy density in a granular medium might have a more complicated response than what can be explained by simple Weibull size scaling. This might arise for two reasons. The coordination number of particles varies with the size and shape distribution of particles and the volume fraction of the granular media. This has an influence on the average stress a particle experiences. On the other hand, although larger particles are more prone to failure due to the presence of more defects, the interaction gets more complicated when the defects are comparable to the particle size. This, in addition to other uncertainties, might help account for the initial increase in E_C (from particle size 170 μm to 190 μm in Table 2) with an increase in particle size followed by a general decreasing trend with increasing particle size (from particle size 190 μm to 470 μm in Table 2). In this study, given the uncertainties, an average value of 1.15 MPa has been used for E_C. The sensitivity of the output to this parameter will be evaluated in the latter part of this paper.

Grading index is a metric related to the evolution of fragment size distribution that signifies the proximity of the initial fragment size distribution to the final or critical distribution [47,53–55]. It is calculated using Eq. (19). For most granular media, it has been observed in simulations and experiments that the largest grains do not undergo fragmentation as they are surrounded by many smaller grains leading to a more hydrostatic state of stress often referred to a “cushioning effect” [60,61]. Therefore, it is assumed in breakage mechanics models that while the smallest grains fragment, the largest ones often remain intact. This is contrary to the expectation of larger fragments being weaker due to the statistical size effect [62]. In reality, the interplay between size effect due to the number of defects and coordination number, both of which can be expected to have opposite effects on the probability of fracture with size, control the evolution of grain fragmentation until a steady-state is reached. In case of extremely low porosity systems like the comminuted zone in ceramics, all fragments (or grains) are expected to be sufficiently supported. Hence, the effect of coordination number should not play a role initially, and statistical size effect should dominate the failure response. For such systems, the assumption of the largest fragments (grains) not re-fragmenting might not hold true. Most experiments for granular media are not meant for such low porosity high-pressure systems. Under high-pressure dynamic loading conditions, this can lead to a very high grading index if there is significant fragmentation of the largest grain. In this work, a grading index value of 0.95 has been selected. This is close to the values calculated from Ref. [56] in Table 2. As with all the parameters, the sensitivity of model outputs to grading index will be evaluated in the latter part of this paper.

Sphere indentation experiments of boron carbide were simulated in abaqus. The details of the experimental technique can be found in Ref. [26]. The geometry consists of a 1/4 in. diameter tungsten carbide sphere impacting a boron carbide cylinder, 1 in. diameter and 1 in. height (Fig. 1). Reduced integration eight-noded brick elements (C3D8R) were used to model both the tungsten carbide sphere and the boron carbide cylinder. Similar to Ref. [25], the cylinder was discretized using a mesh size of approximately 0.55 mm, with 99,682 elements. A kinematic contact algorithm for frictionless surfaces was used. The integrated ceramics model was used for the boron carbide cylinder, with properties as listed in Table 1. The tungsten carbide sphere was modeled using the Johnson–Cook material model parameters determined by Ref. [63]. The equation of state values from Ref. [64] was used (see Table 3). The simulations were conducted for 10 different cases of 100 to 1000 m/s impact velocity.

Edge-On Impact experiments were performed by Strassburger on boron carbide and other ceramics. The experimental technique was developed in EMI, and the details of the experiments can be found in Refs. [27–29]. The simulation geometry consists of a steel projectile 30 mm in diameter and 23 mm in height impacting a 100 mm² and 10 mm thick boron carbide plate along the edge (Fig. 2). Reduced integration eight-noded brick elements were used in abaqus for both the plate and the projectile. A kinematic contact algorithm for frictionless surfaces was used. Johnson–Cook model calibration for hardened steel [65] was used for the steel projectile (see Table 4). The integrated ceramics model was used for modeling the boron carbide plate. K_IC and ρ₀ were modified in this model to $2.9MPam$ and 2530 kg/m³, respectively, to be consistent with the values reported in Ref. [29]. All the other parameters have been set to the same values as Table 1. The simulations were run until 7.5 μs after impact for eight different impact velocities from 50 to 1010 m/s.

Typical features observed in Sphere Indentation experiments are a comminuted zone under the impactor with visible radial cracking on the surface. Immediately under the comminuted zone, cone cracks are also observed. In our simulations damage (Fig. 3(a)) and density (Fig. 3(b)), localizations are observed that are presumed to be related to larger-scale radial cracks. Figure 4(a) shows the damage contour along a section of the ceramic cylinder 15 μs after a spherical indenter impact at 600 m/s. Figures 4(b) and 4(c) show the corresponding density contour along two orthogonal planes. High damage and low-density regions under the indenter signify the presence of a comminuted zone. In addition to that, we observe slanted damage and density localizations that have the appearance of cone cracks. Figures 4(b) and 4(c) demonstrate that the slanted low-density regions are repeatable across different orthogonal planes. This confirms that these low-density regions are 3D cone crack-like features, not simply an artifact of one particular cross section. Some of the other features observed in experiments like crack branching and lateral cracking are observed at certain impact velocities but they are not repeated across all scenarios.

Figure 5 shows the damage contour at the top of the cylinder, 15 μs after impact at different velocities. Figure 6 shows the corresponding density contour. Density contour shows more radial crack-like localization than the corresponding damage contour. In either case, the number of radial crack-like features increases with an increase in impact velocity. However, determining the number of radial cracks is a subjective assessment. The number of observed radial cracks is however slightly less than the number of radial cracks observed in Ref. [26] for boron carbide. Radial cracks were not quantified in Refs. [66,67].

Figure 7 shows that both the percentage of material that is granular and the percentage of material that is amorphized, 15 μs after impact, increase with an increase in impact velocity. The percentage of material that undergoes amorphization is, however, insignificant even at an impact velocity of 1000 m/s (less than 0.1$%$⁠). Figure 8 shows that the maximum depth of penetration of the indenter almost linearly increases with an increase in impact velocity. Amorphization does not seem to affect the depth of penetration of the indenter significantly for the impact velocity range studied here. Reference [24] argues that significant amorphization is not observed in sphere indentation experiments until an impact velocity of around 2 km/s. We can expect the effect of amorphization to be more pronounced at higher impact velocities or for different shapes of indenter.

Figure 9 shows the evolution of the surface damage contour with time for 1010 m/s impact velocity in Edge On Impact simulation. The damaged region around the indenter grows with time. Around 2.625 μs after impact, damage starts localizing into cone crack-like features, which mostly become visible around 4 μs. There is a distributed damaged region with a darker or more fragmented granular region inside.

Figures 10 and 11 show a comparison of the experiments performed by Strassburger [27–29] at two different impact velocities of 469 m/s and 1010 m/s. The experimental image (Figs. 10(a) and 11(a)) is based on the intensity of light reflected from the surface of the ceramic plate viewed via a high-speed camera. It is not clear exactly which model output would best represent the changes in the intensity of reflected light. The intuitive options are damage, density, and out of plane displacement at the surface. Figures 10(b) and 11(b) show the damage pattern for 469 m/s and 1010 m/s impact velocity, respectively. The predicted damage pattern differs slightly in the two cases, and there is a longer cone crack-like damage localized feature observed for 1010 m/s impact, along with a wider fragmented or granular region. In this paper, the numerical damage front represents the boundary capturing Ω ≥ 0, and the numerical granular front or the edge of the granular region represents the boundary capturing Ω ≥ Ω_f. The numerical granular front in Fig. 10(b) seems to be around the same location as the experimental damage front in Fig. 10(a). However, in Fig. 11(b), the numerical damage front seems to be in the same location as the experimental damage front in Fig. 11(a). In both Figs. 10(c) and 11(c), the density pattern has been adjusted to show even the slightest amount of dilation. The density pattern in either figures mimics the corresponding numerical damage pattern. Predictions of out of plane strain (Figs. 10(d) and 11(d)) show a marked difference between the two impact velocities. The region experiencing higher out of plane strain grows with increasing impact velocity. However, it is still not clear if out of plane strain best correlates with the change in light intensity observed in experiments.

Figure 12 shows a comparison of the numerical damage contour (bottom row) with the experimental image [27] (top row) at three different time instants. It appears as though the numerical damage front (dotted line) always exceeds the experimental damage front (dashed line), while the numerical granular front (dashed-dotted line) is at or behind the experimental damage front.

Another interesting feature is that the velocity of the damage front in the middle of the plate always exceeds the velocity of the damage front at the surface of the plate (Fig. 13). For an impact velocity of 469 m/s, 4.5 μs after impact, the distance to which granular and damage front at the middle section of the plate (Fig. 13(b)) has progressed exceeds the granular and damage front at the surface of the plate (Fig. 13(a)). Also mid-section damage, orthogonal to the plane of the plate, shows much higher damage at the middle than the surface (Fig. 13(c)).

The numerical damage front velocity and granular front velocity has been calculated and compared against the experimentally reported damage front velocity in Fig. 14.

In the experiments, the damage front velocity initially rises sharply as a function of impact velocity and then plateaus until an impact velocity of around 469 m/s. Beyond this impact velocity, the damage front velocity gradually rises and reaches to around 12 km/s at an impact velocity of 1010 m/s. However, in the simulations, both the numerical damage front velocity and the granular front velocity sharply increase initially, followed by a more gradual increase with further increase in impact velocity. The experimental damage front velocity appears somewhat bounded between the granular and numerical damage front velocities, 6 μs after impact. This is not necessarily true at an early stage when the numerical damage front and granular front almost coincide. The inset image in Fig. 14 shows the numerical granular and damage fronts. Although the simulations do not capture the sudden rise in damage front velocity beyond 469 m/s impact velocity, they feature onset of amorphization around 742 m/s impact velocity, with the volume of amorphized region gradually increasing with further increase in impact velocity (see Fig. 15). Figure 16 shows a comparison of the growth in granular material percentage (GMP) versus amorphized material percentage, 7.5 μs after impact with change in impact velocity. Despite the coincidence of the onset and increase in amorphization with the sudden rise in damage front velocity reported in experiments, it is not clear as to how amorphization can induce a sudden change in the damage front velocity. Reference [45] does not account for new crack nucleation as a consequence of amorphization, but it is not obvious if that would lead to a rise in front velocity. Another possible explanation of the experimentally observed rise in damage front velocity, beyond 469 m/s impact velocity could be some sort of phase transformation that changes material properties. The current integrative model does not capture such phase transformation. To summarize, experimentally observed damage pattern in Edge On Impact experiments is correlated with the growth of cracks, change in density and out of plane displacements observed in simulations using the integrative model. Typical features such as cone cracks, a distributed crack front, and even secondary crack zones (for low impact velocity) observed in experiments [28] can be reproduced in the simulations. However, unlike the simulations, in the experiments for boron carbide, these patterns are not always symmetric and consistently discernible. The numerical damage front in most cases exceeds the experimental damage front, except 1010 m/s impact velocity. The change in experimental damage front velocity with impact velocity is not accurately captured by the numerical damage front.

Ballistic performance of boron carbide can be improved by modifying the mechanical properties via doping and grain boundary engineering. Reducing defect population by controlling the free carbon content and densification techniques [68] can help improve hardness. Significant research toward enhancing fracture toughness, strength, and/or modulus through different sintering aids [69–76] has been conducted. Amorphization mitigation via silicon doping [77,78] and boron enrichment [79,80] has also been explored. However, performance enhancement via controlling granular flow through material modifications is not well understood, and the focus is driven on characterization of granular flow [81]. In addition to this plethora of research serving as a guide for material modifications, ranking the properties with the most pronounced influence over ballistic performance is desirable.

The integrative model currently has a huge parameter space with significant uncertainty around most of the parameters. While many of these parameters are flags to control the onset or suppression of mechanisms, a significant number of these are model parameters which are directly or indirectly related to physical properties. As a result, there is a significant challenge in designing new materials due to several reasons. First, the influence of model parameters and related material properties on model predictions is not well understood. Second, the design of new materials by modifying individual properties poses a fiscal challenge. Performing a set of experiments to calibrate multiple iterations of a material is expensive and laborious. A work-around that addresses both these issues is to understand the sensitivities of model predictions to model parameters. This will not only help us in understanding the trends that guide us towards improving material performance, but also in the selection of the most significant model parameters. The following subsections will describe the setup, results, and conclusions from a sensitivity analysis study of sphere indentation simulations in boron carbide.

Reference [82] performed sensitivity analysis on ballistic impact of silicon carbide (SiC) ceramic plate with poly-ether-ether-ketone (PEEK) layer and selected peak normal contact force, plastic dissipation in ceramic and PEEK, transmitted impulse to the ceramic back face as output quantities of interest (QoIs). Penetration state function defined using the residual bullet velocity and the depth of penetration is used as a QoI in ballistic impact of SiC/ultra-high molecular weight polyethylene composite plate in Ref. [83]. Similarly crater size, number of radial cracks, ejecta velocity are useful QoIs for further studies of ballistic performance of ceramics.

In this paper, sensitivity analysis of sphere indentation simulations has been performed using the same geometry and setup as highlighted in Sec. 3.2. The depth of penetration of the spherical indenter and the granular material percentage, 15 μs after impact, are selected as the two output QoIs.

As seen before in Fig. 7, amorphization does not play a significant role in the range of impact velocities studied. To simplify our problem, amorphization is deactivated and the CBM model for granular mechanics is employed. We have selected 20 parameters which we suspect play an important role. These include four mechanical, six microstructural, and 10 granular flow parameters, shown in Table 5.

The range of these parameters are chosen in order to bound commonly observed micro-mechanical values for different ceramics and granular solids as much as possible. In some cases when the parameter is not well understood, the range was selected based on engineering intuition. Scrambled Sobol sequences [84,85] are used to generate 500 space-filled samples in the 20-dimensional parameter space bounded by the ranges given in Table 5. The integrative model was then run for each of the 500 parameter combinations. From each simulation, the percentage granular material as well as the evolution of the depth of penetration is obtained.

In addition, some of the 20 integrative model parameters are combined to obtain a set of derived parameters which are either non-dimensional or more physically representative. These parameters are summarized in Table 6. After replacing some of the independent parameters,they are associated with a separate study is also performed with a new set of parameters (derived and original).

For an accurate interpretation of the linear regression coefficients, multicollinearity [86] of parameters should be avoided. The linear correlation coefficient heat map is shown in Fig. 17. It shows the negligible correlation among the input parameters, which suggests the desired lack of multicollinearity. This is supported by the variance inflation factor (VIF) [86] values of each parameter shown in Table 7. VIF for a parameter is obtained by regressing that parameter with respect to the other remaining parameters and has a lower bound of 1. The closer the VIF value for a parameter is to 1, the less the dependence of that parameter on the others. VIFs between 1 and 5 [87] suggest there is a moderate correlation but represent an acceptable level of multicollinearity.

The list of 19 derived parameters is shown in Table 8. Minimum half flaw size (s_min), maximum half flaw size (s_max), flaw distribution exponent (α_flaw), flaw density (η), density (ρ₀), and Poisson’s ratio (ν₀) are removed and replaced by mean half size (s_mean), flaw spacing (s_spacing), range half flaw size (s_range), initial damage (Ω_i), and longitudinal wave speed (V_L).

A correlation study, similar to the one with the original parameter set, is performed here. The correlation coefficient heat map is shown in Fig. 18, and the VIF values for each parameter are shown in Table 9. All the VIF values are within the acceptable range of 1–5. This suggests that the parameter combinations have an acceptable level of multicollinearity and linear regression analysis can be reliably performed.

In this section, we study how the different mechanical, microstructural, and granular parameters influence the GMP at 15 μs after impact. The histogram plot for GMP is shown in Fig. 19. The results show that this QoI varies considerably in the sampled parameter space from $∼5–70%$⁠.

The linear correlation coefficients between GMP and each of the input parameters are calculated and shown in Table 10. Fracture toughness (K_IC) has the strongest correlation with GMP and the correlation is negative, which means the GMP decreases with an increase in K_IC. Minimum half flaw size and friction parameters are also strongly correlated with GMP. On the other hand, Poisson’s ratio and the flaw distribution coefficient have the weakest correlation. These correlation coefficients, however, only provide a preliminary crude idea about the sensitivities of the parameters. For a more detailed sensitivity analysis, we resort to linear regression analysis.

Three iterations of linear regression analysis are performed and the corresponding prediction results are shown in Table 11. As part of the prediction results, the coefficient of determination (R²), the adjusted coefficient of determination (adj. R²) and the root-mean-squared error (RMSE) are reported. In the first iteration, analysis is done using the original dataset with sample size 500 and all 20 parameters. For the full data analysis, the R² value is 0.827, adj. R² value is 0.820, and the RMSE value is 0.4159. For the LOOCV analysis, the R², adj. R², and RMSE values are 0.811, 0.803, and 0.435, respectively. In the second iteration, the cross-validated errors for each individual data were obtained from the LOOCV analysis results in the first iteration and used to detect potential outliers based on the residual IQR (Eq. (22)). Eight outliers are detected and then linear regression was performed with a dataset of sample size 492 and the same set of 20 parameters. The elimination of outliers helps improve all the three metrics (R², adj. R², RMSE) for both the full data and LOOCV analysis as shown in Table 11. For example, the R² improved from 0.827 to 0.841 for the full data analysis and from 0.811 to 0.826 for the LOOCV analysis. In the third iteration, a group of eight parameters (ν₀, η, α_flaw, κ_BM, γ_B, γ_d, u, l) are removed such that the full data R² without these parameters is not significantly reduced. Linear regression analysis is then performed with the dataset of sample size 492 and a set of 12 reduced parameters. It is found from the results in Table 11 that for the full data analysis, the R² value reduces from 0.841 to 0.839 and the RMSE value increases from 0.392 to 0.395 but the adj. R² increases slightly from 0.834 to 0.835 which is desirable. For the LOOCV analysis, however, there is an improvement in all the three metrics compared to those in the second iteration.

The standardized t-test statistic of each parameter obtained from the regression analysis in the third iteration is shown in Fig. 20. From the p-values (not shown), it is found that all 12 parameters are statistically significant for a significance level of 0.05. The t-test statistic bar plot in Fig. 20 gives an idea of the sign of correlation between each of the parameters and the output GMP. In the figure, the parameters are arranged such that the importance of parameters increases from top to bottom. It is thus seen that fracture toughness (K_IC) is the most sensitive parameter followed by minimum half flaw size (s_min) and friction parameter (M_BM). On the other hand, the strain rate sensitivity coefficient (N_BM) and maximum porosity without damage (ϕ_u) are among the least sensitive parameters, but still significant.

The linear correlation coefficients between GMP and each of the 19 derived input parameters shown in Table 8 are calculated and reported in Table 12. Fracture toughness (K_IC) has the strongest correlation with GMP (similar to the case in Sec. 6.3.1) while crushability parameter (κ_BM) and range half flaw size (s_range) have the weakest correlation.

Next, linear regression analysis is performed and the corresponding prediction results are shown in Table 13. It is noted that in the third and final iteration of the analysis, a group of five parameters (κ_BM, γ_B, N_BM, u, l) are removed such that the full data R² without these parameters is not significantly reduced. All these five parameters happen to be the granular parameters of the CBM model.

The standardized t-test statistic of each parameter obtained from the regression analysis in the third iteration are shown in Fig. 21. From the p-values, it is found that all the reduced 14 parameters are statistically significant for a significance level of 0.05. Fracture toughness (K_IC) is the most sensitive parameter followed by mean half flaw size (s_mean) and friction parameter (M_BM). On the other hand, dilative behavior parameter (γ_d) and maximum porosity without damage (ϕ_u) are among the least sensitive parameters, but still significant. To sum up the sensitivity analysis of the granular material percentage, similar trends are found in the original and the derived parameter cases with respect to the order of importance of the parameters. For example, in both cases, the three most sensitive parameters and their importance order is very similar, the only difference being the original parameter s_min replaced by the derived parameter s_mean.

One of the outputs of the integrative model is the evolution of the depth of penetration values with time. Two distinct cases are observed in the 500 simulation outcomes. In one case, the sphere rebounds after impacting the cylinder, and in the other case, the sphere continues to penetrate into the cylinder. The maximum depth of penetration for the rebound case (MDR) is selected as the output QoI for sensitivity analysis. To be exact, 342 simulations lead to rebound of the sphere while the rest lead to sphere penetration. Thus, 342 samples are used for the sensitivity analysis of MDR. The histogram for the MDR output is shown in Fig. 22.

The linear correlation coefficients between MDR and each of the input parameters are shown in Table 14 which suggests friction parameter (M_BM) has the strongest correlation with MDR while maximum half flaw size has the weakest correlation with MDR.

Next, linear regression analysis is performed and the corresponding prediction results are shown in Table 15. After removing the outliers and reducing the number of parameters, the R², adj. R², and RMSE values are reported to be 0.918, 0.915, and 0.270, respectively, for the LOOCV analysis.

The standardized t-test statistic of the 11 significant parameters obtained from the regression analysis in the third iteration is shown in Fig. 23. Friction parameter (M_BM) is the most sensitive parameter followed by grading index (⁠$ϑBM$⁠) and fracture toughness (K_IC). On the other hand, flaw orientation (θ_flaw) and maximum porosity without damage (ϕ_u) are among the least sensitive parameters, but still significant.

The linear correlation coefficient values between MDR and each of the derived input parameters (not reported here) are very similar to that shown in Sec. 6.4.1. The prediction results from the linear regression analysis shown in Table 16 indicate a slight improvement in all the three metrics compared to that reported in Sec. 6.4.1. Figure 24 shows the standardized t-test statistic of the 11 significant parameters obtained from the regression analysis in the third iteration. Friction parameter (M_BM) is the most sensitive parameter followed by grading index (⁠$ϑBM$⁠) and fracture toughness (K_IC). On the other hand, crushability parameter (κ_BM) and maximum porosity without damage (ϕ_u) are among the least sensitive parameters, but still significant.

The influence of model parameters on a given model output is a complex interplay of the various mechanisms with which they are associated. Although there are multiple physical mechanisms active at any instant, only a few of them can be intuitively expected to play a crucial role in influencing the output. For example, as discussed before, amorphization may not play a critical role in these particular experiments because of the small amorphization volume observed in the present simulations, but wave propagation, crack growth, and energy dissipation due to granular flow might. We have tried to isolate the parameters by the mechanism they might be responsible for, and determine the possible correlations with model output. Some parameters might be responsible for multiple mechanisms. The suspected influence of model parameters on physical mechanisms and the corresponding correlation with percentage granular region has been highlighted in Table 17. Most of the correlations reported in Tables 10 and 12 match the physical intuition in Table 17. However, the suspected correlations from ρ₀, ν₀, α_flaw, μ_flaw, γ_d do not match the observed correlations. This might be due to complicated parameter interactions. At any rate, the magnitude of correlation corresponding to each of these parameters is very small, and they are not parameters that have a significant influence on model QoIs.

As mentioned in Sec. 6.4, only 342 out of the 500 simulations led to rebound within 15 μs after impact. The simulations in which rebound did not occur within that duration may feature rebound at a later time, or may not feature rebound because of cylinder fragmentation. Hence, the 342 rebound cases were used to conduct the sensitivity study for the depth of penetration. The depth of penetration is a more localized measure, intricately related to deformation mechanisms in the Mescall zone. This region is highly comminuted and almost certainly under granular flow. Unsurprisingly, the model parameters associated with granular mechanics seem to play a more significant role. In addition, the percentage of the region under granular flow is expected to have a weak influence as well. When more of the region is granular, the region right under the indenter might exhibit less granular flow. Thus, parameters associated with crack growth might have competing influences. The influence of the modulus, fracture toughness, critical breakage energy density is physically related to overall stiffness and can be expected to be negatively correlated with the depth of penetration. For other granular mechanics model parameters, there is an indirect influence via porosity change and volumetric deformation. Table 18 highlights the suspected correlation of model parameters with the instantaneous depth of penetration. The influence of crack growth parameters on the maximum depth of penetration is complicated and unclear. Once again physical intuition can be successfully used to justify correlations corresponding to wave speed and granular flow parameters except ρ₀, ν₀, and κ_BM.