Abstract

A micromechanical multi-physics model for ceramics has been recalibrated and used to simulate impact experiments with boron carbide in abaqus. The dominant physical mechanisms in boron carbide have been identified and simulated in the framework of an integrated constitutive model that combines crack growth, amorphization, and granular flow. The integrative model is able to accurately reproduce some of the key cracking patterns of Sphere Indentation experiments and Edge On Impact experiments. Based on this integrative model, linear regression has been used to study the sensitivity of sphere indentation model predictions to the input parameters. The sensitivities are connected to physical mechanisms, and trends in model outputs have been intuitively explored. These results help suggest material modifications that might improve material performance, prioritize calibration experiments for materials-by-design iterations, and identify model parameters that require more in-depth understanding.

1 Introduction

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 [510] 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 [1114] 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 [2325] 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 [2729].

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 [3133] 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 [3436] 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 [4144] 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.

2 Overview—Integrative Model

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

2.1 Kinematics and Equation of State.

Multiplicative split of the deformation gradient tensor (F) [49] into deformation associated with micro-crack induced damage (Fed), amorphization (Fa), and granular plasticity (Fgp) is used to model the kinematics as
(1)
The temperature (θ), Hugoniot pressure (pH), Grüneisen coefficient (Γ0), reference density (ρ0), the cold energy(ec), and the specific heat at constant entropy (cη) is used in the Mie Grüneisen Equation of State to compute the pressure of the intact solid (pS) as
(2)

2.2 Amorphization.

Amorphization is modeled as parallel bands along which sliding occurs [45]. There are three primary parts that describe the phenomenon. This includes a criterion for initiation of amorphization bands followed by sliding along these bands and damage induced due to these bands which in turn affects the transition to granular flow and the rate of crack growth due to degradation of the critical stress intensity factor. The damage induced due to amorphization is calculated as
(3)
where γs is the shear deformation, ka and na are material parameters associated with the number density of failed bands, and d is the band spacing.
It has been assumed that damage induced due to amorphization linearly degrades the fracture toughness (KIC) of the material as shown in Eq. (4).
(4)
where KICeq is the equivalent effective fracture toughness and Dc is the critical failure damage.

2.3 Fracture and Fragmentation.

The model, developed by Ref. [23], incorporates defect distribution in the material as micro-cracks and is based on classical wing-crack growth from sliding flaw models of Refs. [50,51]. Macroscale material variability has been addressed in Ref. [23] by generating microstructural realizations of local flaw distribution. Crack interactions are modeled using an effective medium approach incorporating dynamic crack growth [46]. A non-dimensional damage parameter D [23,46] is used to calculate the properties of the effective medium. The damage parameter is a summation of the damage induced due to amorphization Da and micro-cracking Dm (Eq. (5)).
(5)
where Nbins represent the number of flaw families with ωk being the number of flaws in each flaw family having a representative flaw size sk. The degradation of the elastic properties of the medium with growth of damage is calculated as
(6a)
(6b)
with Zr=16(1ν02)/3E0(1(ν0/2)), Zn=16(1ν02)/3E0, Zc = (−Zn)/8.

K0, G0, and ν0 are the bulk modulus, shear modulus, and the Poisson’s ratio of the undamaged material.

The local compliance matrix is used to calculate the local stress state of the material. A dynamic crack growth criterion based on Ref. [52] is used to calculate the crack velocity (l˙) as
(7)
where Vmax is the maximum crack velocity, KI is the stress intensity factor, KIC is the fracture toughness, and γc is the crack growth exponent which is a fitting parameter. The crack velocity is then used to compute incremental crack growth and the rate of damage.

2.4 Transition to Granular Mechanics.

The damage parameter serves as a switch between a damaged continuum to a granular continuum determined by the threshold damage (Ωf). This threshold damage is computed using the transition equations proposed by Ref. [19]. Reference [19] estimates the extent of fragmentation at a damage threshold using a numerical crack coalescence model. For the onset of granular mechanics, sufficient fragmentation has to occur. Based on a parametric study, empirical transition equations were proposed, that correspond to a threshold degree of fragmentation. The empirical equations suggest the transition to granular mechanics occurs when the effective wing crack length reaches a certain threshold (lf) determined by the initial flaw size (li), effective initial flaw density (η), and initial flaw orientation distribution. For a fixed defect orientation, the transition wing crack length (lfd) criterion is expressed as
(8)
Equation (8) can be reformulated in terms of a damage-based transition criterion, where the critical damage (Ωfd) is a function of the initial damage (Ωi) as
(9)
Similarly, for a random flaw orientation distribution, the transitional criterion can be expressed in terms of a transition wing crack length (lfr) as Eq. (10) or a transition damage (Ωfr) as Eq. (11).
(10)
(11)

2.5 Continuum Breakage Mechanics Model.

A rate-dependent constitutive model developed by Refs. [47,48] based on breakage theory in geomechanics [5355] and previously implemented in an integrated ceramics model in Ref. [25], is used to model the continuum deformation of granular media. Overstress theory of viscoplasticity has been used to introduce rate dependency in the model. The yield surface is defined as
(12)
where EB is breakage energy and is physically related to the energy dissipation due to particle breakage, B is the breakage index which denotes how close the fragment distribution is to the ultimate distribution, τRP is the relative porosity, and p and q are the pressure and deviatoric stress, respectively. EC, γd, Md, and MBM are material and model parameters. EC is the critical breakage energy density and is related to the strain energy density at the onset of comminution, γd is related to the behavior associated with dilation, Md is a dilatancy parameter, and MBM is the friction parameter. Energy dissipation due to particle breakage, reorganization of particles and friction dissipation is accounted for in the model.
The model captures the refragmentation of the initial fragments during granular flow via the evolution of breakage (B) as
(13)
The evolution of porosity (ϕ) and plastic strain (ɛp) is given by
(14)
(15)
where λ˙BM is a non-negative multiplier defined by Eq. (16) where NBM is the strain rate sensitivity coefficient and ηBM is the viscosity parameter. γB and κBM are material parameters that control the initiation and evolution of breakage. H(FBM) is the Heaviside step function of FBM defined in Eq. (17). H(FBM) = 1 if FBM ≥ 0, otherwise H(FBM) = 0.
(16)
(17)
The constitutive relationship is defined by
(18)
where D is the stiffness and ϑBM is the grading index. ϑBM denotes how far the initial distribution (d0) is from the ultimate fragment size distribution (du). It is calculated from the ratio of the second moments of the initial and final fragment size distributions as,
(19)

3 Numerical Simulations—Integrative Model

3.1 Calibration of Model Parameters.

Calibration of mechanical, equation of state, microstructural, micromechanical, amorphization and granular flow parameters follow [24,25]. Critical transition damage for initiation of granular mechanics [2325] is newly recalibrated in this work as per [19]. Critical breakage energy density (EC) 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.

3.1.1 Calibration of Transition Damage.

Damage parameter is a summation of damage due to wing crack growth (Eq. (5)) and amorphization-induced damage. When the damage exceeds the threshold set in Ref. [19], granular physics is activated in the model. The constant critical threshold damage used in Refs. [2325] is able to capture some structural fragmentation (macroscale fragments), but does not capture microscale fragments sufficiently, at or smaller than the resolution of integrative model simulations. This does not properly represent the initiation of granular mechanics in the Mescall region as explained in Ref. [19]. The current model uses the threshold defined in Ref. [19]. Flaw sizes are assumed to follow a Pareto distribution with the minimum half flaw size (smin), maximum half flaw size (smax), and the distribution exponent (αflaw) defining the distribution. The mean-squared average of initial flaw size can be expressed as
(20)
From the transitional criteria (Eqs. (9) and (11)), the transitional wing crack damage for fixed defect orientation distribution (Ωfd) and random defect orientation distribution (Ωfr) is estimated.

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

3.1.2 Calibration of Critical Breakage Energy Density.

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, EC [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, pCR (Eq. (5.7) in Ref. [55]) at which the breakage of fragments initiates. The grading index, critical comminution pressure, and the corresponding bulk modulus (Kg) are used to calculate the value of EC. 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 EC (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 EC. The sensitivity of the output to this parameter will be evaluated in the latter part of this paper.

3.1.3 Calibration of Grading Index.

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

3.2 Simulation of Sphere Indentation Experiments.

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.

3.3 Simulation of Edge On Impact Experiments.

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. [2729]. The simulation geometry consists of a steel projectile 30 mm in diameter and 23 mm in height impacting a 100 mm2 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. KIC and ρ0 were modified in this model to 2.9MPam and 2530 kg/m3, 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.

4 Results—Integrative Model

4.1 Sphere Indentation Simulations.

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.

4.2 Edge On Impact Simulations.

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

4.3 Necessity of Parameter Prioritization.

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

5 Problem Setup—Sensitivity Analysis

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 this study, a linear regression model is fit to the selected QoIs, and statistical inference is performed to understand the relative importance between each of the parameters and the QoIs. If x~=[x1,x2,,xn] denotes the set of n input parameters and y denotes the QoI, then the linear regression model assumes a linear relationship between x~ and y given by
(21)
where b~=[b0,b1,b2,,bn] is the vector of unknown n + 1 coefficients which needs to be determined. After an initial analysis with the original dataset, leave-one-out-cross-validation (LOOCV) is performed and the individual cross-validated errors are used to detect potential outliers based on the residual interquartile range (IQR) given by
(22)
where Q3res is the third quartile and Q1res is the first quartile of the residual values of each data point from the LOOCV analysis. After removing the outliers, linear regression is again performed on the remaining dataset. A parameter reduction is then performed by eliminating a group of parameters which does not reduce the coefficient of determination (R2) significantly compared to the one obtained with all the parameters. Linear regression is then performed on the remaining dataset with the reduced parameters.

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

6 Results—Sensitivity Analysis

6.1 Original Parameter Set: Correlation Study.

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.

6.2 Derived Parameter Set: Correlation Study.

The list of 19 derived parameters is shown in Table 8. Minimum half flaw size (smin), maximum half flaw size (smax), flaw distribution exponent (αflaw), flaw density (η), density (ρ0), and Poisson’s ratio (ν0) are removed and replaced by mean half size (smean), flaw spacing (sspacing), range half flaw size (srange), initial damage (Ωi), and longitudinal wave speed (VL).

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.

6.3 Granular Material Percentage.

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 570%.

6.3.1 Analysis With Original Parameter Set.

The linear correlation coefficients between GMP and each of the input parameters are calculated and shown in Table 10. Fracture toughness (KIC) has the strongest correlation with GMP and the correlation is negative, which means the GMP decreases with an increase in KIC. 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 (R2), the adjusted coefficient of determination (adj. R2) 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 R2 value is 0.827, adj. R2 value is 0.820, and the RMSE value is 0.4159. For the LOOCV analysis, the R2, adj. R2, 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 (R2, adj. R2, RMSE) for both the full data and LOOCV analysis as shown in Table 11. For example, the R2 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 (ν0, η, αflaw, κBM, γB, γd, u, l) are removed such that the full data R2 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 R2 value reduces from 0.841 to 0.839 and the RMSE value increases from 0.392 to 0.395 but the adj. R2 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 (KIC) is the most sensitive parameter followed by minimum half flaw size (smin) and friction parameter (MBM). On the other hand, the strain rate sensitivity coefficient (NBM) and maximum porosity without damage (ϕu) are among the least sensitive parameters, but still significant.

6.3.2 Analysis With Derived Parameter Set.

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 (KIC) has the strongest correlation with GMP (similar to the case in Sec. 6.3.1) while crushability parameter (κBM) and range half flaw size (srange) 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, NBM, u, l) are removed such that the full data R2 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 (KIC) is the most sensitive parameter followed by mean half flaw size (smean) and friction parameter (MBM). 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 smin replaced by the derived parameter smean.

6.4 Depth of Penetration.

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.

6.4.1 Analysis With Original Parameter Set.

The linear correlation coefficients between MDR and each of the input parameters are shown in Table 14 which suggests friction parameter (MBM) 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 R2, adj. R2, 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 (MBM) is the most sensitive parameter followed by grading index (ϑBM) and fracture toughness (KIC). On the other hand, flaw orientation (θflaw) and maximum porosity without damage (ϕu) are among the least sensitive parameters, but still significant.

6.4.2 Analysis With Derived Parameter Set.

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 (MBM) is the most sensitive parameter followed by grading index (ϑBM) and fracture toughness (KIC). On the other hand, crushability parameter (κBM) and maximum porosity without damage (ϕu) are among the least sensitive parameters, but still significant.

7 Discussion

7.1 Physical Mechanisms and Trends.

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 ρ0, ν0, α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 ρ0, ν0, and κBM.

7.2 Implications Toward Designing Materials.

Table 19 lists the top 10 most sensitive parameters from regression studies for percent granular material and depth of penetration using original and derived model parameters. In either case, fracture toughness (KIC), granular friction (MBM), shear modulus (G0), grading index (ϑBM), minimum flaw size (smin) seem to be important. The strain rate sensitivity coefficient of granular flow (NBM) is important for the depth of penetration. This study provides us insights that could guide material processing modifications. It is difficult if not impossible to control some of these parameters, particularly the granular mechanics parameters other than MBM and EC. Fortunately, these parameters are not the most significant ones. As mentioned earlier, researchers have explored different techniques to improve KIC, G0 [7476]. Although these parameters have been assumed to be independent, some of them might be related to one another through available processing routes, and they may be challenging to independently control. For example, it might not be possible to vary the defect population without affecting the polycrystalline fracture toughness or the modulus.

The initial damage (Ωi) is physically related to the volume fraction of defects. While processing a ceramic, Ωi might be optimized to meet a certain performance goal. Further optimization to improve performance would likely involve controlling the defect size and the defect spacing. Our study suggests that the flaw density (η) or flaw spacing (sspacing) are not as significant as the minimum flaw size (smin) or the mean flaw size (smean). This might mean that ensuring smaller closely spaced defects might be more desirable than larger agglomerates. Flaw distribution parameters are related to the secondary phases in boron carbide matrix. The most abundant of these secondary phases is free carbon. Location of free carbon along boundaries of fragments suggests crack growth from these sites [88,89]. Therefore, controlling the size and volume fraction of these graphitic inclusions [90] can help address fracture mechanisms. Smaller defect spacing might also lead to smaller initial fragments. In addition to this, smaller fragments might also increase EC, due to the size effect, which might improve impact performance. Higher granular friction (MBM) is one of the most desirable traits and more angular particles can achieve that. However, it is not clear how to control angularity. Reference [15] suggests that larger fragments have lower circularity. However, it might not be fair to compare a bulk-averaged estimate of granular friction with individual particle shape. Reference [91] investigated the influence of particle morphology on frictional behavior of sand. Further research toward microstructural features that control the angularity of subsequent fragments and therefore granular friction will be useful.

7.3 Implications Towards Modeling and Calibration.

The results from the sensitivity analysis suggest that certain parameters, some of which control the evolution of state variables in the CBM model (Eqs. (1315)), are not significant for the output quantities of interest studied in this work. If the same can be replicated for simulation of other impact experiments (summarized in Ref. [92]), it might be assumed that either the model can be simplified to ignore those parameters or that those parameters do not need recalibration with slight changes in the material. For example, many of the granular mechanics parameters are calibrated through multiple drained and undrained triaxial compression tests, oedometric compression tests on granular solids [48]. Often the classical geomechanics experimental setups cannot be employed at the high-pressure conditions of impact experiments. So, for a new material, while it will be ideal to recalibrate granular friction (MBM) using pressure shear impact experiments [81], one can rely on past data for γB, γd, κBM. Similarly, accurate estimation of polycrystalline fracture toughness (KIC) is essential [93], and should be prioritized over flaw friction (μflaw) or relative density tests to calibrate l, u and dry density test to calibrate ϕu. This can save both computational and experimental effort and expenditure.

8 Summary

A ceramics model that integrates multiple physical mechanisms has been recalibrated and used to simulate sphere indentation and Edge On Impact experiments in ABAQUS using boron carbide as a model material. The simulations are able to replicate key cracking patterns observed in experiments through damage and density localizations. Two simulation outputs from sphere indentation experiments have been identified as quantities of interest for a sensitivity analysis study: % granular material and indentation depth, 15 μs after impact. 20 micro-mechanical and granular flow model parameters have been varied to generate 500 space-filled samples. Linear regression analysis for the two quantities of interest has been conducted to identify the most significant model parameters. Connections to physical mechanisms for these model parameters have been argued and the implications towards material design from the sensitivity analysis have been explored. The results from the sensitivity study assists in prioritizing calibration experiments for new materials.

Footnotes

Acknowledgment

Research was sponsored by the Army Research Laboratory and was accomplished under Cooperative Agreement Number W911NF-12-2-0022. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the Army Research Laboratory or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation herein.

The HPC at the Maryland Advanced Research Computing Center (MARCC)2 and the currently decommisioned COPPER cluster at the DoD HPC3 were used for conducting the simulations.

The authors acknowledge the contributions of Prof. J. D. Hogan, Prof. K. T. Ramesh, Prof. Nilanjan Mitra, Prof. Mark Robbins, and others in the CMEDE4 Ceramics Modelling group through insightful discussions.The authors would also like to thank Dr. Elmar Strassburger for permission to reproduce the experimental images. Copyright for Fig. 10(a) was obtained from the publisher of Ref. [29], John Wiley and Sons, under license number 4996030533359.

Data Availability Statement

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

References

1.
Crouch
,
I.
,
Franks
,
G.
,
Tallon
,
C.
,
Thomas
,
S.
, and
Naebe
,
M.
,
2017
, “7 – Glasses and Ceramics,”
The Science of Armour Materials
,
I. G.
Crouch
, ed.,
Woodhead Publishing
,
Duxford, UK
, pp.
331
393
.
2.
Clemmer
,
J. T.
,
2019
, “
An Improved Computational Constitutive Model for Brittle Materials
,” Ph.D. Thesis, Physics, Johns Hopkins University, Baltimore, MD. https://jscholarship.library.jhu.edu/handle/1774.2/62208
3.
Daux
,
C.
,
Moës
,
N.
,
Dolbow
,
J.
,
Sukumar
,
N.
, and
Belytschko
,
T.
,
2000
, “
Arbitrary Branched and Intersecting Cracks With the Extended Finite Element Method
,”
Int. J. Numer. Methods Eng.
,
48
(
12
), pp.
1741
1760
. https://doi.org/10.1002/1097-0207
4.
Falk
,
M. L.
,
Needleman
,
A.
, and
Rice
,
J. R.
,
2001
, “
A Critical Evaluation of Cohesive Zone Models of Dynamic Fracture
,”
J. Phys. IV France
,
11
(
PR5
), pp.
Pr5–43
Pr5–50
. https://doi.org/10.1051/jp4:2001506
5.
Bažant
,
Z. P.
, and
Oh
,
B.
,
1983
, “
Crack Band Theory for Fracture of Concrete
,”
JMatériaux et Construction
,
16
(
3
), pp.
155
177
. 10.1007/BF02486267
6.
Le
,
J.-L.
, and
Eliáš
,
J.
,
2016
, “
A Probabilistic Crack Band Model for Quasibrittle Fracture
,”
ASME J. Appl. Mech.
,
83
(
5
), p.
051005
. 10.1115/1.4032692
7.
Spatschek
,
R.
,
Brener
,
E.
, and
Karma
,
A.
,
2011
, “
Phase Field Modeling of Crack Propagation
,”
Philos. Mag.
,
91
(
1
), pp.
75
95
. 10.1080/14786431003773015
8.
Hofacker
,
M.
, and
Miehe
,
C.
,
2012
, “
Continuum Phase Field Modeling of Dynamic Fracture: Variational Principles and Staggered FE Implementation
,”
Int. J. Fracture
,
178
(
1–2
), pp.
113
129
. 10.1007/s10704-012-9753-8
9.
Borden
,
M. J.
,
Verhoosel
,
C. V.
,
Scott
,
M. A.
,
Hughes
,
T. J.
, and
Landis
,
C. M.
,
2012
, “
A Phase-Eield Description of Dynamic Brittle Fracture
,”
Comput. Methods Appl. Mech. Eng.
,
217–220
, pp.
77
95
. 10.1016/j.cma.2012.01.008
10.
Schlüter
,
A.
,
Willenbücher
,
A.
,
Kuhn
,
C.
, and
Müller
,
R.
,
2014
, “
Phase Field Approximation of Dynamic Brittle Fracture
,”
Comput. Mech.
,
54
(
5
), pp.
1141
1161
. 10.1007/s00466-014-1045-x
11.
Johnson
,
G. R.
, and
Holmquist
,
T. J.
,
1994
, “
An Improved Computational Constitutive Model for Brittle Materials
,”
AIP. Conf. Proc.
,
309
(
1
), pp.
981
984
. 10.1063/1.46199
12.
Ravichandran
,
G.
, and
Subhash
,
G.
,
1995
, “
A Micromechanical Model for High Strain Rate Behavior of Ceramics
,”
Int. J. Solids. Struct.
,
32
(
17
), pp.
2627
2646
. 10.1016/0020-7683(94)00286-6
13.
Deshpande
,
V.
, and
Evans
,
A.
,
2008
, “
Inelastic Deformation and Energy Dissipation in Ceramics: A Mechanism-Based Constitutive Model
,”
J. Mech. Phys. Solids.
,
56
(
10
), pp.
3077
3100
. 10.1016/j.jmps.2008.05.002
14.
Fernández-Fdz
,
D.
,
Zaera
,
R.
, and
Fernández-Sáez
,
J.
,
2011
, “
A Constitutive Equation for Ceramic Materials Used in Lightweight Armors
,”
Comput. Struct.
,
89
(
23
), pp.
2316
2324
. 10.1016/j.compstruc.2011.08.003
15.
Hogan
,
J. D.
,
Farbaniec
,
L.
,
Daphalapurkar
,
N.
, and
Ramesh
,
K.
,
2016
, “
On Compressive Brittle Fragmentation
,”
J. Am. Ceram. Soc.
,
99
(
6
), pp.
2159
2169
. 10.1111/jace.14171
16.
Mescall
,
J.
, and
Weiss
,
V.
,
1983
,
Material Behavior Under High Stress and Ultrahigh Loading Rates
, Vol.
29
,
Springer
,
New York, United States
. https://www.springer.com/gp/book/9780306414749
17.
Curran
,
D.
,
Seaman
,
L.
,
Cooper
,
T.
, and
Shockey
,
D.
,
1993
, “
Micromechanical Model for Comminution and Granular Flow of Brittle Material Under High Strain Rate Application to Penetration of Ceramic Targets
,”
Int. J. Impact Eng.
,
13
(
1
), pp.
53
83
. 10.1016/0734-743X(93)90108-J
18.
Farbaniec
,
L.
,
Hogan
,
J.
,
Xie
,
K.
,
Shaeffer
,
M.
,
Hemker
,
K.
, and
Ramesh
,
K.
,
2017
, “
Damage Evolution of Hot-Pressed Boron Carbide Under Confined Dynamic Compression
,”
Int. J. Impact Eng.
,
99
, pp.
75
84
. 10.1016/j.ijimpeng.2016.09.008
19.
Bhattacharjee
,
A.
,
Hurley
,
R.
, and
Graham-Brady
,
L.
,
2020
, Predicting High Rate Granular Transition and Fragment Statistics at the Onset of Granular Flow for Brittle Ceramics, Nov. arXiv:2011.08331.
20.
Chen
,
M.
,
McCauley
,
J. W.
, and
Hemker
,
K. J.
,
2003
, “
Shock-Induced Localized Amorphization in Boron Carbide
,”
Science
,
299
(
5612
), pp.
1563
1566
. 10.1126/science.1080819
21.
Huq
,
F.
,
Liu
,
J.
,
Tonge
,
A.
, and
Graham-Brady
,
L.
,
2019
, “
A Micromechanics Based Model to Predict Micro-Crack Coalescence in Brittle Materials Under Dynamic Compression
,”
Eng. Fract. Mech.
,
217
, p.
106515
. 10.1016/j.engfracmech.2019.106515
22.
Shockey
,
D. A.
,
Marchand
,
A.
,
Skaggs
,
S.
,
Cort
,
G.
,
Burkett
,
M.
, and
Parker
,
R.
,
1990
, “
Failure Phenomenology of Confined Ceramic Targets and Impacting Rods
,”
Int. J. Impact Eng.
,
9
(
3
), pp.
263
275
. 10.1016/0734-743X(90)90002-D
23.
Tonge
,
A. L.
, and
Ramesh
,
K.
,
2016
, “
Multi-scale Defect Interactions in High-Rate Brittle Material Failure. Part I: Model Formulation and Application to ALON
,”
J. Mech. Phys. Solids.
,
86
, pp.
117
149
. 10.1016/j.jmps.2015.10.007
24.
Zeng
,
Q.
,
Tonge
,
A. L.
, and
Ramesh
,
K.
,
2019
, “
A Multi-Mechanism Constitutive Model for the Dynamic Failure of Quasi-Brittle Materials. Part II: Integrative Model
,”
J. Mech. Phys. Solids.
,
131
, pp.
20
42
. 10.1016/j.jmps.2019.06.015
25.
Cil
,
M. B.
,
Zeng
,
Q.
,
Hurley
,
R. C.
, and
Graham-Brady
,
L.
,
2020
, “
An Integrative Model for the Dynamic Behavior of Brittle Materials Based on Microcracking and Breakage Mechanics
,”
J. Dyn. Behav. Mater.
,
6
(
4
), pp.
472
88
. 10.1007/s40870-020-00251-x
26.
Leavy
,
R. B.
,
Brannon
,
R. M.
, and
Strack
,
O. E.
,
2010
, “
The Use of Sphere Indentation Experiments to Characterize Ceramic Damage Models
,”
Int. J. Appl. Ceram. Technol.
,
7
(
5
), pp.
606
615
. 10.1111/j.1744-7402.2010.02487.x
27.
Strassburger
,
E.
,
2002
, Investigation of Fracture Propagation During Impact in Boron Carbide. Technical Report, United States Army, European Research Office of the U.S. Army, 01.
28.
Strassburger
,
E.
,
2004
, “
Visualization of Impact Damage in Ceramics Using the Edge-On Impact Technique
,”
Int. J. Appl. Ceram. Technol.
,
1
(
3
), pp.
235
242
. 10.1111/j.1744-7402.2004.tb00175.x
29.
Strassburger
,
E.
,
2013
, “Edge‐On Impact Investigation of Fracture Propagation in Boron Carbide,”
Advances in Ceramic Armor
, Vol.
IX
,
J. C.
LaSalvia
,
S.
Kirihara
, and
S.
Widjaja
, eds.,
John Wiley & Sons, Ltd
,
New Jersey
.
30.
Borgonovo
,
E.
, and
Plischke
,
E.
,
2016
, “
Sensitivity Analysis: A Review of Recent Advances
,”
Eur. J. Oper. Res.
,
248
(
3
), pp.
869
887
. 10.1016/j.ejor.2015.06.032
31.
Saltelli
,
A.
,
Ratto
,
M.
,
Andres
,
T.
,
Campolongo
,
F.
,
Cariboni
,
J.
,
Gatelli
,
D.
,
Saisana
,
M.
, and
Tarantola
,
S.
,
2008
,
Global Sensitivity Analysis: The Primer
,
John Wiley & Sons
,
West Sussex, UK
.
32.
Saltelli
,
A.
,
Annoni
,
P.
,
Azzini
,
I.
,
Campolongo
,
F.
,
Ratto
,
M.
, and
Tarantola
,
S.
,
2010
, “
Variance Based Sensitivity Analysis of Model Output. Design and Estimator for the Total Sensitivity Index
,”
Comput. Phys. Commun.
,
181
(
2
), pp.
259
270
. 10.1016/j.cpc.2009.09.018
33.
Campbell
,
J. E.
,
Carmichael
,
G. R.
,
Chai
,
T.
,
Mena-Carrasco
,
M.
,
Tang
,
Y.
,
Blake
,
D. R.
,
Blake
,
N. J.
,
Vay
,
S. A.
,
Collatz
,
G. J.
,
Baker
,
I.
, and
Berry
,
J. A.
,
2008
, “
Photosynthetic Control of Atmospheric Carbonyl Sulfide During the Growing Season
,”
Science
,
322
(
5904
), pp.
1085
1088
. 10.1126/science.1164015
34.
Helton
,
J. C.
,
1993
, “
Uncertainty and Sensitivity Analysis Techniques for Use in Performance Assessment for Radioactive Waste Disposal
,”
Reliab. Eng. Syst. Saf.
,
42
(
2–3
), pp.
327
367
. 10.1016/0951-8320(93)90097-I
35.
Kleijnen
,
J. P.
,
2010
, “Sensitivity Analysis of Simulation Models.” Wiley Encyclopedia of Operations Research and Management Science.
36.
Errico
,
R. M.
,
1997
, “
What is An Adjoint Model?
Bull. Am. Meteorol. Soc.
,
78
(
11
), pp.
2577
2592
. 10.1175/1520-0477(1997)078<2577:WIAAM>2.0.CO;2
37.
IM
,
S.
,
1993
, “
Sensitivity Estimates for Nonlinear Mathematical Models
,”
Math. Model. Comput. Exp
,
1
(
4
), pp.
407
414
.
38.
Homma
,
T.
, and
Saltelli
,
A.
,
1996
, “
Importance Measures in Global Sensitivity Analysis of Nonlinear Models
,”
Reliab. Eng. Syst. Saf.
,
52
(
1
), pp.
1
17
. 10.1016/0951-8320(96)00002-6
39.
Chatterjee
,
S.
, and
Hadi
,
A. S.
,
1988
,
Sensitivity Analysis in Linear Regression
(
Wiley Series in Probability and Statistics
),
Wiley
,
Hoboken, NJ
.
40.
Yan
,
X.
, and
Su
,
X. G.
,
2009
,
Linear Regression Analysis
,
World Scientific
,
Singapore
.
41.
Buhmann
,
M. D.
,
2003
,
Radial Basis Functions: Theory and Implementations
, Vol.
12
,
Cambridge University Press
,
Cambridge, MA
.
42.
Bhaduri
,
A.
, and
Graham-Brady
,
L.
,
2018
, “
An Efficient Adaptive Sparse Grid Collocation Method Through Derivative Estimation
,”
Probab. Eng. Mech.
,
51
, pp.
11
22
. 10.1016/j.probengmech.2017.11.002
43.
Bhaduri
,
A.
,
He
,
Y.
,
Shields
,
M. D.
,
Graham-Brady
,
L.
, and
Kirby
,
R. M.
,
2018
, “
Stochastic Collocation Approach With Adaptive Mesh Refinement for Parametric Uncertainty Analysis
,”
J. Comput. Phys.
,
371
, pp.
732
750
. 10.1016/j.jcp.2018.06.003
44.
Bhaduri
,
A.
,
Brandyberry
,
D.
,
Shields
,
M. D.
,
Geubelle
,
P.
, and
Graham-Brady
,
L.
,
2020
, “
On the Usefulness of Gradient Information in Surrogate Modeling: Application to Uncertainty Propagation in Composite Material Models
,”
Probab. Eng. Mech.
,
60
, p.
103024
. 10.1016/j.probengmech.2020.103024
45.
Zeng
,
Q.
,
Tonge
,
A. L.
, and
Ramesh
,
K.
,
2019
, “
A Multi-Mechanism Constitutive Model for the Dynamic Failure of Quasi-Brittle Materials. Part I: Amorphization As a Failure Mode
,”
J. Mech. Phys. Solids.
,
130
, pp.
370
392
. 10.1016/j.jmps.2019.06.012
46.
Paliwal
,
B.
, and
Ramesh
,
K.
,
2008
, “
An Interacting Micro-Crack Damage Model for Failure of Brittle Materials Under Compression
,”
J. Mech. Phys. Solids.
,
56
(
3
), pp.
896
923
. 10.1016/j.jmps.2007.06.012
47.
Cil
,
M. B.
,
Hurley
,
R. C.
, and
Graham-Brady
,
L.
,
2019
, “
A Rate-Dependent Constitutive Model for Brittle Granular Materials Based on Breakage Mechanics
,”
J. Am. Ceram. Soc.
,
102
(
9
), pp.
5524
5534
. 10.1111/jace.16376
48.
Cil
,
M. B.
,
Hurley
,
R. C.
, and
Graham-Brady
,
L.
,
2020
, “
Constitutive Model for Brittle Granular Materials Considering Competition Between Breakage and Dilation
,”
J. Eng. Mech.
,
146
(
1
), p.
04019110
. 10.1061/(ASCE)EM.1943-7889.0001690
49.
Simo
,
J.
, and
Ortiz
,
M.
,
1985
, “
A Unified Approach to Finite Deformation Elastoplastic Analysis Based on the Use of Hyperelastic Constitutive Equations
,”
Comput. Methods Appl. Mech. Eng.
,
49
(
2
), pp.
221
245
. 10.1016/0045-7825(85)90061-1
50.
Ashby
,
M.
, and
Hallam
,
S.
,
1986
, “
The Failure of Brittle Solids Containing Small Cracks Under Compressive Stress States
,”
Acta Metall.
,
34
(
3
), pp.
497
510
. 10.1016/0001-6160(86)90086-6
51.
Nemat-Nasser
,
S.
, and
Horii
,
H.
,
1982
, “
Compression-Induced Nonplanar Crack Extension With Application to Splitting, Exfoliation, and Rockburst
,”
J. Geophys. Res.: Solid Earth
,
87
(
B8
), pp.
6805
6821
. 10.1029/JB087iB08p06805
52.
Freund
,
L. B.
,
1990
,
Dynamic Fracture Mechanics
(
Cambridge Monographs on Mechanics
),
Cambridge University Press
.
53.
Einav
,
I.
,
2007
, “
Breakage Mechanics–Part I: Theory
,”
J. Mech. Phys. Solids.
,
55
(
6
), pp.
1274
1297
. 10.1016/j.jmps.2006.11.003
54.
Einav
,
I.
,
2007
, “
Breakage Mechanics–Part II: Modelling Granular Materials
,”
J. Mech. Phys. Solids.
,
55
(
6
), pp.
1298
1320
. 10.1016/j.jmps.2006.11.004
55.
Einav
,
I.
,
2007
, “
Fracture Propagation in Brittle Granular Matter
,”
Proc. R. Soc. A: Math., Phys. Eng. Sci.
,
463
(
2087
), pp.
3021
3035
. 10.1098/rspa.2007.1898
56.
Nicewicz
,
P.
,
Peciar
,
P.
,
Macho
,
O.
,
Sano
,
T.
, and
Hogan
,
J. D.
,
2020
, “
Quasi-Static Confined Uniaxial Compaction of Granular Alumina and Boron Carbide Observing the Particle Size Effects
,”
J. Am. Ceram. Soc.
,
103
(
3
), pp.
2193
2209
. 10.1111/jace.16871
57.
Thévenot
,
F.
,
1990
, “
Boron Carbide–A Comprehensive Review
,”
J. Eur. Ceram. Soc.
,
6
(
4
), pp.
205
225
. 10.1016/0955-2219(90)90048-K
58.
Dandekar
,
D. P.
,
2001
, Shock Response of Boron Carbide. Technical Report, Army Research Lab, Aberdeen Proving Ground, MD, USA, 04.
59.
Chocron
,
S.
,
Nicholls
,
A. E.
, and
King
,
N. L.
,
2012
, “
Intact and Predamaged Boron Carbide Strength Under Moderate Confinement Pressures
,”
J. Am. Ceram. Soc.
,
95
(
1
), pp.
350
357
. 10.1111/j.1551-2916.2011.04931.x
60.
Sammis
,
C.
,
King
,
G.
, and
Biegel
,
R.
,
1987
, “
The Kinematics of Gouge Deformation
,”
Pure Appl. Geophys.
,
125
(
5
), pp.
777
812
. 10.1007/BF00878033
61.
Åström
,
J.
, and
Herrmann
,
H.
,
1998
, “
Fragmentation of Grains in a Two-Dimensional Packing
,”
Eur. Phys. J. B - Condens. Matter Complex Syst.
,
5
(
3
), pp.
551
554
. 10.1007/s100510050476
62.
McDowell
,
G.
, and
Amon
,
A.
,
2000
, “
The Application of Weibull Statistics to the Fracture of Soil Particles
,”
Soils and Found.
,
40
(
5
), pp.
133
141
. 10.3208/sandf.40.5_133
63.
Holmquist
,
T.
, and
Johnson
,
G.
,
2005
, “
Modeling the 14.5 Mm BS41 Projectile for Ballistic Impact Computations
,”
Comput. Ballistics II, WIT Trans. on Modelling and Simul.
,
4
(
2087
), pp.
61
75
.
64.
Tonge
,
A. L.
, and
Ramesh
,
K.
,
2016
, “
Multi-Scale Defect Interactions in High-Rate Failure of Brittle Materials, Part II: Application to Design of Protection Materials
,”
J. Mech. Phys. Solids.
,
86
, pp.
237
258
. 10.1016/j.jmps.2015.10.006
65.
Iqbal
,
M.A.
,
Senthil
,
K.
,
Bhargava
,
P.
, and
Gupta
,
N.K.
,
2015
, “
The characterization and ballistic evaluation of mild steel
,”
International Journal of Impact Engineering
,
78
, pp.
98
113
. 10.1016/j.ijimpeng.2014.12.006
66.
LaSalvia
,
J. C.
,
Normandia
,
M. J.
,
Miller
,
H. T.
, and
MacKenzie
,
D. E.
,
2005
, “Sphere Impact Induced Damage in Ceramics: II. Armor‐Grade B4C and WC,”
Advances in Ceramic Armor
, Vol.
26
,
J. J.
Swab
, ed.,
John Wiley & Sons, Ltd.
,
Ohio
.
67.
Aydelotte
,
B.
, and
Schuster
,
B.
,
2016
, “Observation and Modeling of Cone Cracks in Ceramics,”
Dynamic Behavior of Materials
, Vol.
1
,
B.
Song
,
L.
Lamberson
,
D.
Casem
, and
J.
Kimberley
, eds.,
Springer International Publishing
,
Cham, Switzerland
, pp.
19
23
.
68.
Toksoy
,
M. F.
,
Rafaniello
,
W.
,
Xie
,
K. Y.
,
Ma
,
L.
,
Hemker
,
K. J.
, and
Haber
,
R. A.
,
2017
, “
Densification and Characterization of Rapid Carbothermal Synthesized Boron Carbide
,”
Int. J. Appl. Ceram. Technol.
,
14
(
3
), pp.
443
453
. 10.1111/ijac.12654
69.
Kim
,
H.-W.
,
Koh
,
Y.-H.
, and
Kim
,
H.-E.
,
2000
, “
Densification and Mechanical Properties of B4C With Al2O3 As a Sintering Aid
,”
J. Am. Ceram. Soc.
,
83
(
11
), pp.
2863
2865
. 10.1111/j.1151-2916.2000.tb01647.x
70.
Yamada
,
S.
,
Hirao
,
K.
,
Yamauchi
,
Y.
, and
Kanzaki
,
S.
,
2003
, “
High Strength B4C–TiB2 Composites Fabricated by Reaction Hot-Pressing
,”
J. Eur. Ceram. Soc.
,
23
(
7
), pp.
1123
1130
. 10.1016/S0955-2219(02)00274-1
71.
Frage
,
N.
,
Hayun
,
S.
,
Kalabukhov
,
S.
, and
Dariel
,
M. P.
,
2007
, “
The Effect of Fe Addition on the Densification of B4C Powder by Spark Plasma Sintering
,”
Powder. Metall. Met. Ceram.
,
46
(
11–12
), pp.
533
538
. 10.1007/s11106-007-0082-9
72.
Subramanian
,
C.
,
Roy
,
T.
,
Murthy
,
T.
,
Sengupta
,
P.
,
Kale
,
G.
,
Krishnaiah
,
M.
, and
Suri
,
A.
,
2008
, “
Effect of Zirconia Addition on Pressureless Sintering of Boron Carbide
,”
Ceram. Int.
,
34
(
6
), pp.
1543
1549
. 10.1016/j.ceramint.2007.04.017
73.
Sun
,
J.
,
Liu
,
C.
, and
Duan
,
C.
,
2009
, “
Effect of Al and TiO2 on Sinterability and Mechanical Properties of Boron Carbide
,”
Mater. Sci. Eng. A.
,
509
(
1
), pp.
89
93
. 10.1016/j.msea.2009.01.067
74.
Madhav Reddy
,
K.
,
Guo
,
J.
,
Shinoda
,
Y.
,
Fujita
,
T.
,
Hirata
,
A.
,
Singh
,
J.
,
McCauley
,
J.
, and
Chen
,
M.
,
2012
, “
Enhanced Mechanical Properties of Nanocrystalline Boron Carbide by Nanoporosity and Interface Phases
,”
Nat. Commun.
,
3
(
1
), p.
1052
. 10.1038/ncomms2047
75.
Shen
,
Y.
,
Li
,
G.
, and
An
,
Q.
,
2019
, “
Enhanced Fracture Toughness of Boron Carbide From Microalloying and Nanotwinning
,”
Scr. Mater.
,
162
, pp.
306
310
. 10.1016/j.scriptamat.2018.11.035
76.
He
,
Y.
,
Shen
,
Y.
,
Tang
,
B.
, and
An
,
Q.
,
2020
, “
Strengthening Boron Carbide Through Lithium Dopant
,”
J. Am. Ceram. Soc.
,
103
(
3
), pp.
2012
2023
. 10.1111/jace.16889
77.
Proctor
,
J. E.
,
Bhakhri
,
V.
,
Hao
,
R.
,
Prior
,
T. J.
,
Gregoryanz
,
T. S. E.
,
Chhowalla
,
M.
, and
Giulani
,
F.
,
2015
, “
Stabilization of Boron Carbide Via Silicon Doping
,”
J. Phys.: Condens. Matter.
,
27
(
1
), p.
015401
. 10.1088/0953-8984/27/1/015401
78.
Khan
,
A. U.
,
Etzold
,
A. M.
,
Yang
,
X.
,
Domnich
,
V.
,
Xie
,
K. Y.
,
Hwang
,
C.
,
Behler
,
K. D.
,
Chen
,
M.
,
An
,
Q.
,
LaSalvia
,
J. C.
,
Hemker
,
K. J.
,
Goddard
,
W. A.
, and
Haber
,
R. A.
,
2018
, “
Locating Si Atoms in Si-doped Boron Carbide: A Route to Understand Amorphization Mitigation Mechanism
,”
Acta Mater.
,
157
, pp.
106
113
. 10.1016/j.actamat.2018.07.021
79.
Chauhan
,
A.
,
Schaefer
,
M. C.
,
Haber
,
R. A.
, and
Hemker
,
K. J.
,
2019
, “
Experimental Observations of Amorphization in Stoichiometric and Boron-Rich Boron Carbide
,”
Acta Mater.
,
181
, pp.
207
215
. 10.1016/j.actamat.2019.09.052
80.
Schaefer
,
M. C.
, and
Haber
,
R. A.
,
2020
, “
Amorphization Mitigation in Boron-Rich Boron Carbides Quantified by Raman Spectroscopy
,”
Ceramics
,
3
(
3
), pp.
297
305
. 10.3390/ceramics3030027
81.
Sun
,
X.
,
Chauhan
,
A.
,
Mallick
,
D. D.
,
Tonge
,
A. L.
,
McCauley
,
J. W.
,
Hemker
,
K. J.
,
LaSalvia
,
J. C.
, and
Ramesh
,
K.
,
2020
, “
Granular Flow of An Advanced Ceramic Under Ultra-high Strain Rates and High Pressures
,”
J. Mech. Phys. Solids.
,
143
, p.
104031
. 10.1016/j.jmps.2020.104031
82.
Batra
,
R.
, and
Pydah
,
A.
,
2020
, “
Impact Analysis of PEEK/ceramic/gelatin Composite for Finding Behind the Armor Trauma
,”
Composite Struct.
,
237
, p.
111863
. 10.1016/j.compstruct.2020.111863
83.
Shen
,
Z.
,
Hu
,
D.
,
Yang
,
G.
, and
Han
,
X.
,
2019
, “
Ballistic Reliability Study on SiC/UHMWPE Composite Armor Against Armor-Piercing Bullet
,”
Composite Struct.
,
213
, pp.
209
219
. 10.1016/j.compstruct.2019.01.078
84.
Matoušek
,
J.
,
1998
, “
On the L2-discrepancy for Anchored Boxes
,”
J. Complex.
,
14
(
4
), pp.
527
556
. 10.1006/jcom.1998.0489
85.
Joe
,
S.
, and
Kuo
,
F. Y.
,
2003
, “
Remark on Algorithm 659: Implementing Sobol’s Quasirandom Sequence Generator
,”
ACM Trans. Math. Softw. (TOMS)
,
29
(
1
), pp.
49
57
. 10.1145/641876.641879
86.
James
,
G.
,
Witten
,
D.
,
Hastie
,
T.
, and
Tibshirani
,
R.
,
2013
,
An Introduction to Statistical Learning
, Vol.
112
,
Springer
,
New York
.
87.
Sheather
,
S.
,
2009
,
A Modern Approach to Regression With R
,
Springer Science & Business Media
,
New York
.
88.
Farbaniec
,
L.
,
Hogan
,
J.
, and
Ramesh
,
K.
,
2015
, “
Micromechanisms Associated with the Dynamic Compressive Failure of Hot-pressed Boron Carbide
,”
Scr. Mater.
,
106
, pp.
52
56
. 10.1016/j.scriptamat.2015.05.004
89.
Xie
,
K. Y.
,
Kuwelkar
,
K.
,
Haber
,
R. A.
,
LaSalvia
,
J. C.
, and
Hemker
,
K. J.
,
2016
, “
Microstructural Characterization of a Commercial Hot-Pressed Boron Carbide Armor Plate
,”
J. Am. Ceram. Soc.
,
99
(
8
), pp.
2834
2841
. 10.1111/jace.14295
90.
Bakhshi
,
M.
,
Souri
,
A.
, and
Amin
,
M. K.
,
2019
, “
Carbothermic Synthesis of Boron Carbide With Low Free Carbon Using Catalytic Amount of Magnesium Chloride
,”
J. Iranian Chem. Soc.
,
16
(
6
), pp.
1265
1272
. 10.1007/s13738-019-01602-9
91.
Alshibli
,
K. A.
, and
Cil
,
M. B.
,
2018
, “
Influence of Particle Morphology on the Friction and Dilatancy of Sand
,”
J. Geotech. Geoenviron. Eng.
,
144
(
3
), p.
04017118
. 10.1061/(ASCE)GT.1943-5606.0001841
92.
Li
,
H.
,
Lo
,
C.
,
Toussaint
,
G.
,
Sano
,
T.
, and
Hogan
,
J. D.
,
2020
, “Dynamic Fracture and Fragmentation of Boron Carbide,”
Boron Carbide: Structure, Processing, Properties and Applications
,
Kolan Madhav
Reddy
, ed.,
Nova Science Publishers
,
New York
.
93.
Salem
,
J.
,
Quinn
,
G.
, and
Jenkins
,
M.
,
2005
, “Measuring the Real Fracture Toughness of Ceramics: ASTM C 1421,”
Fracture Mechanics of Ceramics
, Vol.
6
,
R. C.
Bradt
,
D.
Munz
,
M.
Sakai
, and
K. W.
White
, eds.,
Springer
,
New York
, pp.
531
553
.