Multi-scale defect interactions in high-rate brittle material failure. Part I: Model formulation and application to ALON

Journal of the Mechanics and Physics of Solids 86 (2016) 117–149

Contents lists available at ScienceDirect

Journal of the Mechanics and Physics of Solids journal homepage: www.elsevier.com/locate/jmps

Multi-scale defect interactions in high-rate brittle material failure. Part I: Model formulation and application to ALON Andrew L. Tonge, K.T. Ramesh n Department of Mechanical Engineering, The Johns Hopkins University, 223 Latrobe Hall, 3400 N. Charles St., Baltimore, MD 21218, USA

a r t i c l e i n f o

abstract

Article history: Received 30 April 2014 Received in revised form 5 October 2015 Accepted 12 October 2015 Available online 23 October 2015

Within this two part series we develop a new material model for ceramic protection materials to provide an interface between microstructural parameters and bulk continuum behavior to provide guidance for materials design activities. Part I of this series focuses on the model formulation that captures the strength variability and strain rate sensitivity of brittle materials and presents a statistical approach to assigning the local flaw distribution within a specimen. The material model incorporates a Mie–Grüneisen equation of state, micromechanics based damage growth, granular flow and dilatation of the highly damaged material, and pore compaction for the porosity introduced by granular flow. To provide initial qualitative validation and illustrate the usefulness of the model, we use the model to investigate Edge on Impact experiments (Strassburger, 2004) on Aluminum Oxynitride (AlON), and discuss the interactions of multiple mechanisms during such an impact event. Part II of this series is focused on additional qualitative validation and using the model to suggest material design directions for boron carbide. & 2015 Elsevier Ltd. All rights reserved.

Keywords: A. Microcracking A. Dynamic fracture B. Ceramic material B. Constitutive behavior

1. Introduction Failure processes in brittle materials are fundamentally linked to the nature and distribution of defects within the materials. In many brittle materials, the controlling defects are crack-like flaws. By understanding the interactions of a population of crack-like defects, researchers have developed models for the evolution of damage parameters (e.g. Nemat-Nasser and Horii, 1982), which capture the history dependence of the strength and failure of brittle materials. The more specific problem of the failure of brittle materials subjected to impact loading is an important problem in geology and geophysics, planetary science, and defense. In order to address these types of problems, the failure mechanisms that occur during an impact event must be captured. Capturing these failure mechanisms requires mechanism-based models. Most of the micromechanics models for failure under compressive loading (Paliwal and Ramesh, 2008; Grechka and Kachanov, 2006; Huang and Subhash, 2003; Basista and Gross, 1998; Nemat-Nasser and Obata, 1988; Ashby and Hallam, 1986; Deshpande et al., 2011; Deshpande and Evans, 2008) provide for the growth of a population of cracks based on some measure of the crack tip driving force and a crack kinetics law. The models differ in how crack interactions are handled. This is typically done either through the use of a crack array (e.g. Deng and Nemat-Nasser, 1992) or through an effective medium approach (e.g. Budiansky and O'Connell, 1976). The distribution of cracks is important for capturing the correct scaling response of strength with strain rate (Kimberley et al., 2013). n

Corresponding author. E-mail addresses: [email protected] (A.L. Tonge), [email protected] (K.T. Ramesh).

http://dx.doi.org/10.1016/j.jmps.2015.10.007 0022-5096/& 2015 Elsevier Ltd. All rights reserved.

118

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

As the microcracks grow, they will eventually intersect creating numerous fragments; once this occurs, the material is reasonably described as a granular material (Brannon et al., 2009). The behavior of granular materials has been extensively studied within the context of geologic materials and the civil engineering community (Terzaghi et al., 1996). Recent work has focused on developing models for granular materials based on explicit modeling of grain interactions using the discrete element method (Andrade et al., 2012). The granular flow contributes two key aspects to an impact event. First, as the granular material is sheared, it “bulks” and increases the pressure on the surrounding material (Bolton, 1986). Second it provides a mechanism for dissipating additional energy (Shockey et al., 1990). Some of the largest microcracks grow and become unstable macro-scale cracks. A number of methods have been developed to explicitly track crack fronts and crack paths, including the cohesive element method (Falk et al., 2001), and the extended finite element method (Daux et al., 2000). These methods are very effective for tracking a few dynamically interacting cracks; however, they become computationally expensive when many dynamically interacting cracks are considered, and some have poor scaling performance in parallel environments (Radovitzky et al., 2011). Techniques such as element deletion (Pandolfi and Ortiz, 2012) can provide successful representations of regions that look like, and behave similar to, explicit cracks at lower computational cost. Alternatively, one could model cracks as damage zones where the material behavior is described by the behavior of the highly damaged material as in the Kayenta model (Brannon et al., 2009) and the brittle material models by Holmquist and Johnson (2011), Johnson and Holmquist (1993) and Johnson et al. (2003). In this work we develop a model for brittle material failure during impact events, such as the Edge On Impact experiments (EOI) conducted on AlON (Strassburger, 2004; Strassburger et al., 2006). In these events the structural dimensions are typically much larger than any of the microscale features. This separation of scales necessitates the use of a homogenized damage approach because of the prohibitive computational cost required to resolve the cohesive zone for the microcracks while still capturing the macroscopic loading. The key physical processes that the model presented in this work captures are summarized in Fig. 1. Over the years a number of models have been developed that incorporate internal variables that are directly related to microstructural features (e.g. Clayton, 2010; Deshpande and Evans, 2008; Deshpande et al., 2011; Gailly and Espinosa, 2002; Paliwal and Ramesh, 2008). In all of these models, the key microstructural feature is a distribution of cracks. Although these models acknowledge the statistical nature of brittle materials, they do not explicitly incorporate the statistical variation of these microstructural quantities as a function of position within a simulation volume. Other authors (Meyer and Brannon, 2012; Brannon and Gowen, 2014) have incorporated a statistical variation of material parameters (such as strength) but did not relate the material parameters to microstructural features (although the idea and justification for the variation of the material parameters was based on the existence of a heterogeneous microstructure at some fine unresolved scale). In this work we use the volume of particle domains within the Convected Particle Domaine Interpolation (CPDI) technique (Sadeghirad et al., 2011) and a Poisson process within each of these particle domains to create local realizations of the flaw distribution, which changes the local constitutive response. Variability is a key feature of brittle materials. Introducing this variability has the added benefit of stabilizing the damage problem and providing natural initiation sites for damage (Molinari et al., 2007; Meyer and Brannon, 2012; Leavy et al., 2009; Brannon and Gowen, 2014). Micromechanics-based damage models have the potential to introduce physically based variability into the constitutive response. Initial work has linked the physical variability observed in brittle materials to the process of sampling flaws in a material and the interactions of these flaws (Graham-Brady, 2010). In this work, we present a mechanism based material model for use in high rate simulations that provides a physically based approach for producing macroscale (computationally resolved) variability based on microstructural sampling. Additionally, we look at the implications of the variability for structural problems. This paper is organized as follows. We first discuss flaw distributions and an approach to generating realizations of flaw distributions in the context of numerical simulations. This discussion is followed by a detailed development of the new material model. We then look at some results from the model, in two parts. We first discuss results for simple uniaxial compression simulations and determine the influence of variability and specimen size on the rate dependent strength. Finally we use the model to look at an edge on impact experiment, providing initial qualitative validation, and then discuss

Fig. 1. Physical processes captured in the material model.

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

119

the interaction of multiple mechanisms. Finally, we consider the implications of our results.

2. Incorporating flaw distributions into simulations Material variability is particularly important for simulations of brittle materials. In such materials, the defects that constitute the defect/flaw distribution include: pores, micro-cracks, and inclusions. These flaws may have both size and orientation distributions; for this work, we focus on the size distribution. Since the computational discretization is determined by the geometry and macroscopic loading conditions, it is important to consider the effects of sampling the defect population. We assume that microstructural features such as grain boundaries and grain orientations are homogenized into effective properties, which include the stiffness, fracture toughness, and crack growth laws. These effective properties are constant throughout the material. For each material (defined by specific composition and processing, e.g. PAD Boron Carbide), there is a distribution of defects that characterize the material. Within the micromechanics model discussed in Section 3.4, the flaws are approximated as disk shaped cracks and the flaw size corresponds to the disk radius. Some approaches for mapping complex crack geometries to an effective disk crack are discussed in Sevostianov and Kachanov (2002). One could also derive an effective crack size by computing the disk radius required to achieve the same stress intensity factor from the disk as the defect produces. This work seeks to determine the response of a material after the flaw population has been converted to an equivalent population of disk flaws. For this population of disk flaws, there is an associated parent flaw size probability distribution g(s) (defined below) and parent flaw density η; these characterize the material. The parent flaw size probability density function (g(s), where s is the radius of an equivalent disk shaped flaw) determines the probability (P) of finding a flaw of size between a lower limit sl and an upper limit sh through the relation:

P=

∫s

sh

g (s ) ds.

(1)

l

The cumulative distribution function (CDF) for flaw size is defined as G (s ) =

s

∫0 g (a) da . Using the CDF we can rewrite the

probability that a given flaw has size between flaw sizes sl and sh as P (sl ≤ s < sh ) = G (sh ) − G (sl ). Knowing the CDF simplifies computing the fraction of a flaw distribution that lies between two limiting flaw sizes. ^ The parent flaw density (η) describes the expected number of flaws per unit reference volume. In a finite volume (V ) of ^ material we can define a local flaw density η^ as follows. If the number of flaws in V is Nf, then we define the local flaw Nf density η^ = ^ . In general, the local flaw density is not equal to the parent flaw density (η^ ≠ η ). By creating a normalized V histogram of the flaw sizes in the local flaw population, we can also calculate the local flaw distribution g^ (s ). For sufficiently large volumes (or numbers of flaws), η^ and g^ (s ) converge to the parent flaw density and the parent flaw distribution function. ^ Consider a computational discretization of the sampled volume V . This discretization introduces another volume associated with the size of each discretization volume. We give this discretization volume the label V0. In general, the finite size of V0 introduces variability in the number of flaws as discussed above. Within each discretization volume (V0), we capture the effect of microcracks using a homogenized damage model. Our scalar damage definition is the sum of the cubed crack sizes (si) over all Nf cracks in V0 divided by V0:

D=

1 V0

Nf

∑ si3. i=1

(2)

When there is a large number of flaws, computing this quantity directly becomes unrealistic. To reduce the computational cost of the damage parameter, we group similar sized flaws in Nbins flaw families (or “bins”). Each of the Nbins represents the h Nk flaws with sizes between slk and sk in the volume V0 using the representative flaw size sk and the family flaw density Nk ωk = V . Using this binning approach, the damage parameter in Eq. (2) can be approximated using 0

N bins

D=

∑ k=1

ωk sk3.

(3)

Each discretization volume within the simulation has its own local flaw distribution determined by the collection of flaw families at that location. Different discretization volumes will in general have different values for the family densities in that discretization volume. We treat both the flaw family density and the representative flaw size for a family as statistical quantities. To simplify the discussion we introduce the notation where an over tilde (s˜k ) denotes a statistical quantity; an overbar (s¯k ) denotes the expected value of the quantity, and the quantity alone (sk) denotes a realization of that quantity. As an approximation for the spatial heterogeneity of the flaw distribution, we assume that the locations of flaws are independent of each other and of the flaw size, and model the spatial distribution of flaws using a Poisson process. The expected value of the family flaw density ω¯ k is the parent flaw density (η) multiplied by the probability that a given flaw size h is between sk and slk (ω¯ k = η (G (skh ) − G (skl ))). From the assumption of a Poisson process, the statistical number of flaws in each

120

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

Fig. 2. Algorithm for generating a simulated microstructural realization of flaws in a material. The material is defined by the parent flaw size distribution (g (s)) and parent flaw density (η) and the geometry is discretized into particles with volume V0 for the Material Point Method. The microstructure is simulated by generating the local flaw distribution at every particle.

family in V0 becomes a Poisson distributed random variable defined by

N˜ k = Pois ⎡⎣ V0 ω¯k ⎤⎦.

(4)

Since the particle volume is fixed, the family flaw density ω˜ k has a distribution that depends on only the distribution of N˜ k . N 1 We define sk as the sample mean of the flaw sizes within the family resulting in sk = N ∑i =k1 ai , where ai are the sizes in k the family (for an extended discussion of this choice see Tonge, 2014, chapter 2). Each flaw size ai within the family is a statistical quantity described by a probability function gk(a), which is related to the parent flaw distribution g(s) through

gk (a) =

G (skh )

⎧ g (a) s l ≤ a < s h 1 k k ⎨ . l − G (sk ) ⎩ 0 otherwise

(5)

For large values of Nk the distribution of s˜k becomes Gaussian with the mean given by the first moment of gk:

s¯k =

∞

∫−∞ agk (a) da,

(6)

and the variance is given by the variance of gk divided by Nk: ∞

σ s2k =

∫−∞ a2gk (a) da − s¯k2 Nk

.

(7)

After discretizing the physical problem of interest, we create a realization of the flaw distribution by assigning sk and ωk for each family in every particle. The algorithm that we use to generate these realizations is shown in Fig. 2. By simulating the flaw family density and representative flaw family size for each family in each particle in the simulation, we generate a ^ realization of the size and spatial distribution of flaws within the simulation volume V . ˜ ˜ To ensure a reasonable computational cost associated with simulating Nk and sk , and to avoid numerical errors due to ˜ finite precision arithmetic when evaluating e−Nk for large values of N˜ k (numerical underflow) which may occur when generating a realization from a Poisson distribution, we introduce a cutoff number of flaws Ncutoff . When the expected

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

121

number of flaws in family k exceeds Ncutoff , we use a Gaussian approximation to the Poisson distribution to generate a realization of the number of flaws in the family Nk; otherwise, we directly sample a Poisson distribution using the procedure in Knuth (1970, p. 132). After creating a realization of Nk for a flaw family in a specific subvolume, we test Nk and determine if it exceeds Ncutoff ; if so, we create a realization of the representative flaw size (sk) for that flaw family in that subvolume using a Gaussian distribution with a mean given by Eq. (6) and variance given by equation (7); otherwise, we explicitly simulate Nk realizations from the flaw size distribution for the k flaw family (gk(a)) and set sk equal to the sample mean of these values. Using Ncutoff = 20 provides a balance between accuracy and computational cost. In Tonge (2014, chapter 2), Tonge shows that changing Ncutoff does not noticeably change the simulated flaw distribution. Additionally, very large values of Ncutoff can cause numerical underflow issues in the algorithm for creating a realization from a Poisson distribution. Introducing Ncutoff limits the computational cost (CPU time) of assigning a family density and family flaw size to each flaw family; however, there is still an incremental cost for each flaw family. The computational cost for each time step scales with the number of flaw families with a non-zero family flaw density. For most flaw distributions of interest, there is a small number of large flaws and a large number of small flaws. We bias the size of the flaw families so that they are larger for large flaws and smaller for small flaws so that for a given number of bins, the biased flaw families should give us a better representation of the flaw distribution (for a discussion of this approach see Tonge, 2014, chapter 2). This algorithm (Fig. 2) takes the computational discretization and the parent flaw distribution as inputs and computes a ^ specific realization of the distribution of flaws within the sample volume V . Since each discretization volume V0 is associated ^ with a position X , the local flaw distribution at specific locations within the sampled volume V is a function of position. We next discuss a constitutive model accounting for the initial flaw distribution.

3. A constitutive model that includes flaw distributions We develop a model that incorporates finite deformation kinematics, thermodynamic effects through an equation of state, micromechanics based damage growth, and granular flow of the highly damaged material. 3.1. Kinematics Impact events produce large deformations in the material. We choose to use the framework outlined in Simo and Ortiz (1985). We use the multiplicative split of the deformation gradient into an elastic part F e and a plastic part F vp such that F = F eF vp , where F vp represents the granular viscoplastic flow of the highly damaged material. We introduce the plastic deformation tensor C vp = F vpT F vp , which is a material measure of the plastic deformation. The elastic finger tensor (be−1) is a spatial measure of the elastic deformation, it follows that we can write be in terms of material tensors as be = FC vp−1FT . This identity is important when we define an objective integration procedure for the evolution of be with time following Simo and Ortiz (1985). We also introduce a multiplicative split of the volumetric and volume preserving deformation. The volume change ratio is the Jacobian of the deformation gradient J = det (F ) =

ρ0 . ρ

We similarly define the volume change ratio due to granular flow

and the elastic volume change ratio = det (F e ). The volume preserving deformation measures (F¯ , F¯ e , and JGP = vp vp −1/3 vp F ). We can then F¯ ) are obtained by dividing by the cube root of the associated Jacobian (F¯ = J−1/3 F , F¯ e = Je−1/3 F e , F¯ = JGP e e vp vp −2/3 T − 2/3 T ¯ ¯ ¯ ¯ ¯ define the corresponding symmetric deformation measures C , C , and b as C = J F Fvp , and b¯ = J−2/3 Fe FT . F F, C = J

det (F vp)

Je

GP

The total rate of deformation tensor d is the symmetric part of the spatial velocity gradient l = 1

vp

vp

∂v . ∂x

e

e

We define the rate of

vp

deformation associated with the viscoplastic granular flow as d vp = 2 F e [(F ̇ F vp−1)T + F ̇ F vp−1] F eT , and the elastic rate of deformation is then given by de = d − d vp . These kinematic definitions are similar to those in Simo and Ortiz (1985), except that our material model is not plastically incompressible so we account for the volume change associated with the plastic flow. 3.2. Deviatoric elastic response Experimental observations suggest that materials support much larger elastic deformations resulting in a compressive volume change than elastic deformations resulting in a shape change. This difference in response leads us to adopt a similar decomposition in our constitutive relation. We divide the Kirchhoff stress (τ ) into a hydrostatic portion (−ps J e ) and a deviatoric portion (τdev )

τ = τdev − ps Je I .

(8)

The Kirchhoff stress is related to the Cauchy stress (σ ) through the determinant of the deformation gradient τ = Jσ . Com1 puting the pressure as − 3 tr (σ ) = p results in p = ps /JGP . In this material model the granular flow may have some dilatation associated with it, and therefore JGP is in general not equal to 1 and Je is in general not equal to J. The volume change associated with granular plasticity JGP is identical to the distension (usually given the label α) discussed in Carroll and Holt

122

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

(1972). The subscript in ps represents the pressure in the solid material if it was subjected to the volume change described by Je. Since there is limited information regarding the deviatoric elastic response of ceramic materials subjected to large elastic deformations, we assume a simple linear form that uses a finite deformation strain measure. This form is motivated by a generalized Neo-Hookean model where the volumetric and isochoric deformations are decoupled and reduces to the infinitesimal strain linear elastic response in the limit of infinitesimal deformation (Simo and Ortiz, 1985). The deviatoric part of the Kirchhoff stress tensor is assumed to be a linear function of the shear modulus G and the isochoric measure of elastic e deformation b¯ given by

τdev = G (b¯ e −

1 3

tr (b¯ e ) I ).

(9)

The coupling with damage is discussed after defining the volumetric response. 3.3. Equation of state (EOS) We use a Mie–Grüneisen equation of state with the assumption of a linear shock speed-particle velocity relationship and artificial viscosity to stabilize the shock front. Our implementation is based on the discussion in Drumheller (1998). From the Rankine–Hugoniot shock jump conditions and a linear relation between the shock speed and the particle velocity one can derive the pressure on the principal Hugoniot as a function of the elastic deformation (Je), the bulk wave speed at room temperature and pressure (C0), the density at room temperature and pressure (ρ0), and the slope of the shock speed curve (S):

⎧ ρ C 2 (1 − J e ) 0 0 ⎪ J e < 1.0 ⎪ 2 pH (J e ) = ⎨ ( 1 − S (1 − J e ) ) . ⎪ 2 ⎪ e otherwise ⎩ ρ0 C0 (1 − J )

(10)

The Mie–Grüneisen equation of state introduces the Grüneisen parameter (Γ) to quantify the energy coupling between pressure and temperature through the atomic vibrations. We make the simplifying assumption that the quantity Γρ = Γ0 ρ0 is a constant (see Drumheller, 1998 for a discussion of the implications of this choice). The total energy is decomposed into a cold energy ec, which is a function of only the volumetric elastic deformation, and a thermal energy, which depends on only the specific heat cv and the temperature θ. A detailed discussion of the derivation of the equation of state (along with the definition of the cold energy) is contained in Appendix A. The resulting equation of state is

⎡ ⎞⎤ Γ ⎛ ps (J e , θ ) = pH (J e ) ⎢ 1 − 0 ⎜ 1 − J e ⎟ ⎥ + ρ0 Γ0 ⎡⎣ ec (J e ) + c η (θ − θ 0 ) ⎤⎦. ⎣ ⎠⎦ 2⎝

(11)

We use the equation of state parameters reported in a recent review on AlON (McCauley et al., 2009). The reported bulk wave speed was 7670 m/s and the slope S was 1.3 with a density of 3595 kg/m3. We assume a Grüneisen parameter Γ0 of 1.6 based on the formula Γ0 ≈ 2S − 1. The specific heat (Cv) and thermal conductivity (αθ ) are 800 J/(kg K) and 9 W/(m K) respectively (Warner et al., 2005). In this coupled thermo-mechanical problem, we split the solution procedure into an isentropic mechanical step followed by a thermal step under fixed mechanical conditions (Armero and Simo, 1993). Since we are interested in short time evolution, we assume adiabatic conditions and skip the heat transfer calculation. Since the pressure and temperature are ̇ ) coupled through the Grüneisen parameter (Eq. (11)), isentropic changes in volume result in a change in temperature (θent given by

̇ = θent

−θΓ0 ρ0 tr (d e ). ρ

(12)

Although we do not track shock fronts in the computational scheme, we capture the dissipative nature of shocks by introducing an “artificial” bulk viscosity based on the form proposed in Wilkins (1980). We implement the artificial viscosity by introducing an additional viscous pressure pvisc that is activated under volumetric compression and provides an additional dissipation mechanism given by

⎧ A C | tr (d )| dx + A tr (d )2dx2 ρ tr (d ) ≤ 0 ( 10 ) e 2 e e pvisc = ⎨ . ⎪ ⎩0 otherwise ⎪

(13)

Here dx is the average edge length for the cell, and the parameters A1 ¼0.4 and A2 ¼4.0 were chosen to smooth the shock front over 9 cells and eliminate oscillations behind the shock. Computational implementation of the equation of state and artificial viscosity are verified by comparison to the analytic solution of a plate impact scenario for a material with negligible shear strength (Tonge, 2014, chapter 2). We assume that all of the energy dissipated by the artificial viscosity is converted into heat resulting in a temperature rise given by

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

123

Fig. 3. Multi-scale approach including the wing crack geometry for the self consistent calculation. The boundary loads (s1 and s3) bridge from the macroscale to the microscale. The crack length l bridges back to the macroscale through the damage parameter D. The initial boundary value problem is discretized at the macroscale. At this scale we assume an isotropic response and apply the granular flow and pore compaction models. The effective medium at the microscale is treated as an anisotropic material to develop tension along e^2 in response to compression along e^1 within the elliptical inclusion.

̇ = θav

−Jpvisc tr (d ) ρ0 cv

(14)

where pvisc is defined in Eq. (13). We next discuss the coupling between the elastic behavior and damage through a micromechanics approach. 3.4. Micromechanics of dynamic fracture We begin by extending (Grechka and Kachanov, 2006) to account for a flaw orientation distribution function and the effect of interacting microcracks. We then discuss the model for dynamic damage growth, which takes into account crack dynamics. The elastic response discussed in the previous section assumes an isotropic material behavior. However, at the micromechanics scale we must address the induced anisotropy associated with the cooperative growth of microcracks. Schematically, this is shown in Fig. 3. The macroscale potato is treated as an isotropic medium while the microscale calculation utilizes the local stress state to define the orientation of the micromechanics calculation. The computational discretization occurs at the macroscale in a three dimensional continuum. At each material point, we compute the damage evolution using a micromechanics based damage model. This damage model computes the effective interaction of the surrounding cracks using an ellipse imbedded within an effective matrix. In three dimensions the crack would be embedded in an ellipsoid instead of an ellipse. However, we only have an analytic solution for the two dimensional micromechanics problem, and so we do not address the complicated three dimensional wing cracking process. For the micromechanics calculation we assume plane strain conditions and pick the two extreme (maximum and minimum) principal stresses as boundary conditions, thus providing a local two dimensional approximation. We use the effective medium approach to account for crack interactions following Paliwal and Ramesh (2008). This provides some interaction between the subscale cracks, but the interaction is not as strong as the interactions in formulations based on an array of cracks like (Deshpande et al., 2011). The presence of cracks increases the elastic compliance of the system by introducing an additional strain, which is proportional to the applied stress. Grechka and Kachanov (2006) demonstrate that the strain energy density associated with a population of Ncracks cracks within a representative volume V0 can be written as

f=

1 1 τ: 0 : τ + τ: 2 2V0

Ncracks

∑ ( n ⊗ Z ⊗ nA)(i) : τ .

(15)

i

Here 0 is the fourth order compliance tensor of the undamaged material, and n i , Zi , and Ai are the unit normal, crack compliance and area of each individual crack. For penny shaped cracks the radial and transverse directions are the same, so the crack compliance can be written in terms of the crack normal and the identity tensor I as

Z=

a ( Zr I − (Zr − Zn ) n ⊗ n). π

(16)

Here Zr and Zn are the compliances in the normal and radial directions normalized by

Zr =

ν02 )

16 (1 − , ⎛ ν ⎞ 3E 0 ⎜ 1 − 0 ⎟ ⎝ 2⎠

Zn =

16 (1 − 3E 0

ν02 )

a π

(Grechka and Kachanov, 2006):

. (17)

The solution to the inclusion problem (see Appendix D) using the strain energy described in Eq. (15) with uniaxial compressive stress boundary conditions (σ1 = 0 and σ 3 < 0) results in a compressive stress state normal to the wing crack plane

124

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

((σ ·e^2 )·e^2 < 0). This is contrary to the results initially reported in Paliwal and Ramesh (2008), where an erroneous factor of 2 in the planar compliance matrix terms related to the shear modulus caused an amplification of the stress intensity factory ((σ ·e^2 )·e^2 > 0) under uniaxial compression as damage increases. Corrections to Paliwal and Ramesh (2008) are discussed in Hu et al. (2015, Appendix B). The compressive stress normal to the crack plane causes the driving force on the crack tip to decrease with increasing crack length. However, experiments (Paliwal et al., 2006) show that as the level of damage grows, interactions between the cracks lead to an acceleration of the damage rate. In addition to an accelerated damage growth rate, brittle materials exhibit a bulking behavior when loaded in compression. We account for both effects by introducing an additional term into the strain energy density function of Eq. (15): f=

1 1 τ: 0 : τ + τ: 2 2V0

Ncracks

∑ ( n ⊗ Z ⊗ nA)(i) : τ + τ: i

1 2V0

Ncracks

∑ i

Z c a i Ai ( n ⊗ n ⊗ I + I ⊗ n ⊗ n − 2n ⊗ n ⊗ n ⊗ n)(i) : τ . π

(18)

The form for this term is based on the representation for a transversely isotropic fourth order tensor (Macon et al., 2014). This additional fourth order tensor provides a way to introduce a Poisson-like effect. The parameter Zc controls the strength of this coupling between loading normal to the crack face and transverse to it. It may be possible to determine a value for Zc based on the stiffness observed in meso-scale simulations such as (Grechka and Kachanov, 2006); however, in this work we choose Zc from the range of thermodynamically admissible values to enhance crack interaction and maintain a decreasing bulk modulus (derived in following paragraphs) with increasing damage. Remembering the scalar damage parameter D, we note that for a delta distribution of flaw sizes arranged in a periodic array when the flaws have grown so that on average a is equal to half of the average spacing between flaws D ¼0.125. When the cracks have grown so that a is equal to the average spacing between flaws, D¼1. In this model, all positive values of D are mathematically permitted, however at some value of D between 0.125 and 1.0 the cracks will strongly interact resulting in fragmentation of the material and a transition to granular flow. We choose D ¼0.125 as the critical condition for the onset of granular flow, and when damage equals 1.0 we stop computing additional damage growth. Physically this choice means that we allow continued fracture within the grains in the granular plasticity model for damage values between 0.125 and 1.0. However when damage reaches 1.0 we assume that there is no further cracking of the grains (and thus no further softening of the elastic moduli). The details of the granular flow model are discussed in the next section. In general, there is a distribution of flaw orientations that can be characterized by the flaw orientation distribution function ρ (n ). Using the damage and the orientation distribution function we can rewrite the strain energy density as an integral over all possible flaw orientations (n ) as f=

1 D 1 τ : 0 : τ + 2 2 4π

∬ω ρ (n) ( Zr τ ·τ: n ⊗ n − (Zr − Zn ) τ: n ⊗ n ⊗ n ⊗ n: τ + Zc τ: ( n ⊗ n ⊗ I + I ⊗ n ⊗ n − 2n ⊗ n ⊗ n ⊗ n): τ ) dω.

(19)

This integral places all of the orientation information in the orientation distribution function ( ρ (n )), and all of the flaw size and density information in D. We now define an equivalent isotropic strain energy function fiso by evaluating Eq. (19) with ρ (n ) = ρiso (n ) = 1. The details of the integration are presented in Appendix B with the final result:

⎞⎞ ⎛ ν ⎞⎞ ⎛ 1 + ν0 D ⎛ D ⎛ 2 fiso = ⎜ + ⎜ Z r − Z n − 8Zc ⎟ ⎟ ( tr (τ ) ) . ⎜ 3Z r + 2Z n − 4Z c ⎟ ⎟ τ : τ − ⎜ 0 + ⎠⎠ ⎝ 2E 0 ⎠⎠ ⎝ 2E 0 30 ⎝ 30 ⎝

(20)

This equivalent strain energy density function is now used at the discretization scale. Although we need to use the anisotropic effective medium in order to compute the crack growth in the micromechanics model, the only information that we use at the discretization scale is the resulting scalar damage parameter. Thus the elastic response at the discretization scale only needs to be isotropic and we use fiso at this scale. In this framework, therefore the compliance at the micromechanics scale and at the computational discretization scale are derived from the same form of the strain energy function, but make different assumptions about the effective local distribution of flaw orientations. Using these two different assumptions can introduce errors when the elastic anisotropy of the problem is important, but is an acceptable approximation in other cases. The loading conditions where the elastic anisotropy may be important are conditions where:

 there are many aligned families of cracks (gas and oil exploration),  there is a single or only a few dominant cracks that grow through cleavage,  or cases where the highly damaged but not fully damaged state under proportional loading is of interest (e.g. interrupted quasi-static compression tests). Under conditions more applicable to impact problems where the failure process occurs rapidly, and the principal stresses will tend to rotate (complex loading histories) we expect a very wide distribution of crack orientations and therefore the approximation of an isotropic distribution of flaws at the computational discretization scale may be reasonable. Accepting the limitations introduced by this approximation, we leave the issue of handling the fully anisotropic problem to future work.

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

125

From this isotropic strain energy density function (Eq. (20)) we compute an effective bulk modulus K(D) and an effective shear modulus G(D):

(

(

K (D) = K0−1 + D Z n + 4Zc

−1

))

,

⎛ ⎞ ⎞−1 2D ⎛ ⎜ 3Z r + 2Z n − 4Z c ⎟ ⎟ . G (D) = ⎜ G0−1 + ⎠⎠ ⎝ 15 ⎝

(21)

Since we are using a Mie–Grüneisen equation of state for the hydrostatic response of the material, we scale the pressure computed using Eq. (11) by the ratio of the damaged bulk modulus (K(D)) to the undamaged bulk modulus (K0). This gives us a pressure in the solid material defined by

ps (J e , θ , D) =

⎡ ⎡ ⎤⎞ ⎞⎤ K (D) ⎛ Γ ⎛ ⎜ p (J e ) ⎢ 1 − 0 ⎜ 1 − J e ⎟ ⎥ + ρ0 Γ0 ⎢ ec (J e ) + c η θ ⎥ ⎟. ⎠⎦ ⎣ ⎣ ⎦⎠ K0 ⎝ H 2⎝

(22)

A sufficient condition to insure that damage growth Ḋ causes non-negative energy dissipation is that the rate of change of ∂K ∂G the shear and bulk moduli with damage are less than or equal to zero ( ∂D ≤ 0 and ∂D ≤ 0). Applying this condition results in the following constraint on the coefficient of the interaction term in Eq. (19):

(23)

− Z n ≤ 4Z c ≤ 3Z r + 2Z n .

For Zc < 0 the interactions between flaws cause a transverse tension as a result of an applied axial compressive stress. This transverse tension accelerates the damage growth rate and is consistent with experimental observations (Paliwal et al., Z 2006). However, at the limit where the tensile interactions are maximized (Zc = − 4n ), the bulk modulus is independent of damage. In tension, it does not make physical sense for the bulk modulus of a damaged material to be independent of damage; as a compromise, we choose an interaction term that gives the expected interactions, while maintaining some softening of the bulk modulus and set Zc equal to one half of its minimum thermodynamically admissible value:

Zc =

−Z n . 8

(24)

This choice of Zc results in an effective Poisson's ratio that decreases only slightly with damage. If measurements or simulations are performed that quantify the anisotropic stiffness of a population of cracks as a function of the crack growth becomes available, then the definition of Zc can be refined. For now it represents a place holder that was not included in the compliance terms defined in Dienes et al. (2006) or Grechka and Kachanov (2006). At the microscale, we are solving a self consistent problem using an effective medium approach as shown in Fig. 3. We select the local coordinate system such that the most compressive principal stress is aligned with e^1 and the most tensile principal stress is aligned with e^2. For wing crack growth, the crack face normals are given by n = e^2. The equivalent orientation density function is a delta function centered at n = e^2 multiplied by 4π. The resulting strain energy density function in the effective medium ( faniso ) is given by

faniso =

1 D DZc τ: 0 : τ + Z r τ ·τ: e^2 ⊗ e^2 − (Z r − Z n ) τ: e^2 ⊗ e^2 ⊗ e^2 ⊗ e^2: τ + τ 2 2 2 : e^2 ⊗ e^2 ⊗ I + I ⊗ e^2 ⊗ e^2 − 2e^2 ⊗ e^2 ⊗ e^2 ⊗ e^2 : τ .

(

(

)

)

(25)

The anisotropic compliance matrix associated with the strain energy function faniso is given by

aniso = 0 + ΔNIA (D) + ΔINT (D).

(26) ΔNIA

Here the first term is the isotropic compliance of the base material, is the additional compliance computed from the non-interacting assumption in Grechka and Kachanov (2006), and ΔINT is the additional compliance term associated with the crack interactions. The isotropic compliance tensor written in terms of the Poisson's ratio and the Young's modulus is 1+ν ν 0 = E 0 s − E0 I ⊗ I , where s is the fourth order symmetric identity tensor such that for any tensor A , the operation s : A 0

0

returns the symmetric part of A . From Grechka and Kachanov (2006) the additional compliance associated with the noninteracting assumption is NIA Δijlm =

8 (1 − ν02 ) p p p . α p δjm + αim δjl + αjlp δim + αjm δil + 4βijml 3E0 (2 − ν0 ) il

(

)

(27)

Here the second order tensor α and the fourth order tensor β are given by α = De^2 ⊗ e^2 and β = − The interaction term is

ΔINT =

−2 (1 − ν02 ) D e^2 ⊗ e^2 ⊗ I + I ⊗ e^2 ⊗ e^2 − 2e^2 ⊗ e^2 ⊗ e^2 ⊗ e^2 . 3E 0

(

)

ν0 ^ De 2 2

⊗ e^2 ⊗ e^2 ⊗ e^2.

(28)

The original self consistent solution for the damage model (Paliwal and Ramesh, 2008) contained an error. We use the corrected solution provided in Hu et al. (2015, Appendix B) to compute the evolution of the damage variable D. For completeness we summarize the important features of the micromechanics growth model in Appendices C and D.

126

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

From Paliwal et al. (2008), we identify the fracture toughness, quasi-static and dynamic strength of AlON, the material used in the Edge-On-Impact experiments (Strassburger et al., 2005) that we simulate later in the paper. We select a fracture toughness of 2.9 MPa m1/2 (Paliwal et al., 2008). We assume that the crack face friction coefficient is similar to contact friction in well rounded granular materials which have a friction angle of close to 30°. This translates into a coefficient of friction (μ) of 0.57 which has a most damaging crack orientation of ϕ = 90° − arctan (1/μ) /2 = 60° (Ashby and Hallam, 1986). We assume that the flaw distribution (g(s)) is a bounded Pareto distribution which has the form

g (s ) =

α s −(α + 1) αsmin . ⎛ smin ⎞α 1−⎜ ⎟ ⎝ smax ⎠

(29)

An exponent of α ¼3 gives self-similar flaw scaling where, in three dimensions, the average distance between flaws longer than a specified size (a) scales linearly with a (Poncelet, 1946). Since we are looking at dynamic failure, we calibrated our material model using a dynamic strength under uniaxial compression of 3.5 GPa at a strain rate of 103 s  1. A flaw size range from smin = 2 μm to smax = 40 μm with a flaw density of 4  1012 m  3 fits the dynamic compressive strength of AlON. This flaw distribution is equivalent to an average spacing of 63 μm between flaws that are at least 2 μm in size. Paliwal et al. (2006) observed some carbonaceous defects on the surfaces of fragments postmortem. These defects had sizes that are consistent with the 2–40 μm flaw size distribution used in this work. The pores and processing defects can occur during any phase of processing. Defects within material grains could result from either grain growth that consumed a foreign particle or precipitation of contaminants during the powder synthesis process. In our model, the flaw density introduces a lower bound on the mesh size given our homogenization approach. A reasonable lower bound is the mesh size corresponding to an average of 8 flaws per discretization volume. For this flaw density that lower bound is close to 125 μm. The structural problems that we are interested in solving have a minimum spatial dimension of 10 mm. It is reasonable to expect that meshes in the sub-millimeter size range are required to capture the dynamics in the problem. Considering these constraints, we recognize that these parameters represent an AlON-like material and may not reflect the behavior of commercial AlON with high purity levels and very low defect densities. In such materials, other deformation and failure mechanisms are likely active. 3.5. Granular plasticity A number of researchers have investigated the dynamic strength of fragmented and granular materials under high confining pressure and dynamic loading (Chen and Luo, 2004; Chocron et al., 2012, 2013; Luo et al., 2006; Nie and Chen, 2011; Martin et al., 2013). In these experiments a linear relationship between pressure and shear strength was observed up to pressures in excess of 1 GPa. Deviations from linearity are often attributed to the activation of additional deformation mechanisms (e.g. plasticity within the grains). In this work we use a two surface model consisting of a Drucker–Prager yield cone representing the frictional response, and a hydrostatic cap representing the effect of pore collapse resulting from grain fracture or plasticity within a grain. Once the pressure is sufficiently large that the hydrostatic compression cap is reached, the observed behavior of the material will no longer be a linear relationship between pressure and shear strength. Single surface models such as the teardrop shape developed by Foster et al. (2005) may be better suited for applications where the response is dominated by the behavior of the damaged material under finite compression and finite shear strains. Future work will examine the use of a single surface three invariant type model for the granular flow and pore compaction behavior of the material; however, the problems of interest for this work avoid high pressures during granular flow and the pore compaction parameters chosen suppress pore compaction in the edge on impact simulations. In this work, we are not driving the material into the extreme compression regime and therefore assume that the deviatoric strength of the granular material is linearly dependent on the hydrostatic pressure. We assume the existence of a yield surface defined by f (τ ) = 0, which defines the onset of plastic flow. We use the common Drucker–Prager yield surface where the deviatoric strength increases linearly with pressure resulting in a yield function defined by

f (τ ) =

⎛ tr (τ ) ⎞ τdev: τdev − Y + A ⎜ − B⎟ . ⎝ 3 ⎠

(30)

Here the parameter B represents the cohesive forces between grains and provides a hydrostatic tensile limit, the parameter A represents the magnitude of the coupling between the hydrostatic and deviatoric components, and the deviatoric yield stress Y represents the deviatoric strength when the hydrostatic term is equal to 0. The parameter A is related to the effective coefficient of friction between grains of material, and for a Mohr–Coulomb type yield surface would be the coefficient of friction. However, because Eq. (30) does not depend on the third stress invariant (either J3 or the Lode angle), A is a scalar multiple of the friction coefficient that is independent of the pressure, but depends on the Lode angle. We use the parameters A¼0.6, B¼0.1 MPa, and Y¼0 based on confined dynamic compression experiments on dry sand (Martin et al., 2013). Granular flow involves the rearrangement of a collection of grains within a representative material volume. Since this rearrangement generally takes some amount of time, we model the granular flow as a viscoplastic process. We decompose the tensorial viscoplastic rate of deformation (d vp ) into a scalar viscoplastic flow rate (λ )̇ and a flow direction (m ):

̇ d vp = λm

(31)

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

127

Assuming that the granular flow is associative to the yield surface f ¼0 results in a granular flow direction given by m = Evaluating the derivative

∂f (τ , ∂τ

J vp )

=

∂f (τ ) ∂τ

∂f . ∂τ

gives

A τdev I. + τdev: τdev 3

(32)

Introducing the deviatoric unit tensor n =

τdev , ∥ τdev ∥

^ and the hydrostatic unit tensor I =

1 3

I , the direction of plastic flow can be

written as

^ m = n + AI .

(33)

The ratio of the bulking deformation to the shear deformation is A, consistent with our assumption of associative flow. We use an associative flow rule for simplicity, but the approach can be extended to allow for non-associative flow. However, the non-associative material model may require additional considerations to avoid unstable growth of “Sandler–Rubin–Pučik” instabilities (Brannon, 2007, Section 6.5). For computational simplicity and to ensure that we maintain positive plastic dissipation in the presence of non-linear effective moduli, we define an effective strain direction meff using

⎧ ⎪ n + ζAI^ tr ( τ ) < 0 m eff = ⎨ . ⎪ ⎩ n + AI^ otherwise

(34)

The original return direction m is used for all of the projections and to compute the effective strain rate discussed below. The parameter ζ starts as 1.0 and is iteratively reduced by multiplication with 0.9 if the computed plastic dissipation is less than e zero. The plastic dissipation is computed by subtracting the strain energy in the final state (i.e., b¯n + 1, JnGP+ 1, Jn + 1) from the strain e GP energy in the trial state (i.e., b¯ , JGP ) plus the thermal energy required to move along an isentropic path from JGP n , J n to J tr

n+1

n+1

at constant J. This thermal term is only non-zero when using an equation of state like the Mie–Grüneisen that couples the mechanical deformation to the temperature. We require this extra step in the solution procedure to ensure non-negative plastic work because most algorithms for pressure dependent plasticity assume an additive decomposition of the strain rate and assume that the stiffness is constant during a time step. Since our problem may violate either of these assumptions, we use the return algorithm for Drucker–Prager plasticity discussed in Brannon (2007) to provide an initial estimate of the return state, then use the iterative procedure described above to ensure that there is non-negative plastic dissipation. This numerical adjustment will reduce effective pressure during granular flow, which could result in more deviatoric granular flow and less bulking than if the actual return direction was used. However, this approach is still more consistent than the approach used in Johnson and Holmquist (1993) where the authors use the pressure from the previous time step to set the deviatoric strength in the current time step. Future work will include updating to a fully non-linear solver for the stress projection that can account for the non-linear material behavior during a time step. The influence of the rate of deformation on granular flow can be examined in a variety of ways. One appealing concept is based on the argument that the arrangement of particles within a granular medium requires some time to reach the lowest energy configuration corresponding to a particular deformation. A representative volume of the granular material will then require a specific stress τ^ to be deformed at an extremely low “quasi-static” rate of deformation so that the particles are able to achieve the lowest energy configuration. However, when the same original representative volume element is taken to the same deformation state at a higher rate of deformation, the particles are unable to reach their lowest energy configuration in the time available, and so the stress τ required to achieve that deformation state is higher. The difference between the stress required at high strain rates and that for quasi-static deformations is defined as the overstress τ¯ = τ − τ^ . We expect that the rate of viscoplastic flow λ ̇ should be a scalar function of the overstress τ¯ and a material timescale τGP . Since the overstress contains both a deviatoric and hydrostatic component, the three standard invariants of τ¯ are not good scalar measures of the energetic driving force for rearranging the particles. Instead, we define a strain-like scalar quantity μ¯ , which measures the magnitude of the elastic strain associated with τ¯ along the direction m

μ¯ m = −1: τ¯.

(35)

We now write the flow rate using the material timescale τGP and μ¯ as

λ̇ =

μ¯ . τGP

(36)

The timescale τGP is essentially a material characteristic timescale defining the time needed for the particles in the RVE of the granular material to relax from the configuration developed at the high rate of loading to the ground state configuration corresponding to quasi-static deformations. Such a material timescale has been defined in the past for granular materials within dynamic loading by Curran et al. (1993), who defined this time from a microscale viewpoint in terms of the average particle size, the packing density, and the shear wave speed (i.e., the timescale arises because of the time needed to communicate across a particle or block within the granular medium). Based on the flaw distribution that we selected the average flaw spacing (η−1/3) is 63 μm. A shear wave traverses this distance in just over 10 ns while a longitudinal wave takes just over 6 ns. In our simulations, we use a timescale for granular flow that is 7 ns (τGP = 7 ns).

128

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

From the definition of the overstress τ¯ it follows that τ^ = τ − τ¯ . Further, at yield, we must have f (τ^ ) = 0. From this condition and Eq. (35), it follows that the overstress τ¯ must satisfy the relation f (τ − τ¯) = 0. We note that the magnitude of the viscoplastic flow rate λ ̇ can now be computed from the condition that f (τ − τ¯) = f (τ −  : (μ¯ m )) = 0 once the characteristic timescale τGP has been defined. For simplicity we treat this characteristic timescale as a constant which is related to the initial flaw density in the material and the shear wave speed, but it is a straight forward extension of the model to make τGP a function of the dilatation and total equivalent granular plastic strain. 3.6. Pore collapse model Our granular plasticity formulation introduces dilatation through the interaction of grains as a result of shear deformation. However, under high confining pressures, we expect a reduction in the porosity. Experimentally, the relationship 1 between the distension JGP and the applied pressure (P = − 3 tr (σ )) is expressed through the crush curve for a given maP −P terial. For simplicity we choose a primitive crush curve defined by J GP = 1 + (J0GP − 1) J2 ( P c− P )2. The material parameter P0 c 0 represents the pressure when inelastic compaction of the porous medium begins for an initial distension of JGP 0 , and Pc represents the pressure required for full densification. We select a reference crush pressure (P0) of 1 GPa with a reference distension (JGP 0 ) of 2.0 and assume that full densification (Pc) occurs at a pressure of 10 GPa. Although these numbers are not based on experimental data, they are reasonable because typical porosities of ceramics prior to sintering are about 50 percent and typical manufacturing pressures (up to a couple hundred MPa) do not significantly reduce this porosity (R. Haber personal communication). Since our material model can develop a distension greater than JGP 0 prior to loading in compression, we require that the crush curve is defined for all values of distension. We extend the crush curve beyond the inelastic compaction limit by assuming an exponential relationship between the distension and the crush pressure (Pα = P0 exp ( − κ (J GP − J0GP ))). The Pc − P 0 crush pressure should be a smooth function of the distension requiring κ = . Writing the crush curve as a yield 2P 0 (J0GP − 1) surface results in:

⎧ ⎛ ⎞ P0 Pc − P0 ⎪ P exp ⎜⎜ − − (J GP − J0GP ) ⎟⎟, P < P0 GP ⎪ Pc − P0 Pc − P0 2P0 (J0 − 1) ⎝ ⎠ ⎪ 2 . fϕ (P , J GP , J ) = ⎨ ⎞ ⎛ ⎪ (J GP − 1) − (J GP − 1) J 2 ⎜ Pc − P ⎟ , P0 ≤ P < Pc 0 ⎝ Pc − P0 ⎠ ⎪ ⎪ GP Pc < P ⎩ J − 1,

(37)

This yield surface represents the quasi-static crush curve. Since we are incorporating dynamic effects into the granular flow portion of the model, we also incorporate dynamic effects into the pore collapse portion of the model. We do this following GP the same procedure as in Section 3.5. The equilibrium distension (JGP eq ) satisfies fϕ (Jeq , J ) = 0. We define a strain like measure GP . We assume the same simple model as Eq. (36) for the rate dependence for the rate of excess distension as ΔJ GP = J GP − Jeq

GP GP −ΔJ GP ̇ ̇ of change of the distension: Jpore when ΔJ GP > 0 and Jpore = τ = 0 otherwise. GP The total rate of change of the distension is the sum of the distension change due to granular flow and the distension change due to pore compaction:

GP

J̇

GP ̇ = J GP tr (d vp) + Jpore

(38)

This pore collapse model, combined with the granular flow model that allows dilatation, produces a steady state porosity (after large shear deformations) that is pressure and rate dependent. The equilibrium porosity may be important for the excavation flow in impact cratering problems and in the flow of the comminuted material during penetration events. Both the pore compaction process and the granular flow process are dissipative processes that can generate a significant amount of heat. Granular flow is primarily a frictional process so we assume that all of the energy dissipated through these processes is converted into heat resulting in a temperature change given by GP

̇ = θGP

̇ Je ( − p ) / ( − tr (τ )) d vp: τ + Jpore s ρ0 cv

.

(39)

Here / is the Heaviside step function that only allows heating due to granular flow if the mean stress is compressive in nature. If it is tensile, then the sub-scale fragments are being disassembled and the only resistance to deformation is their momentum and the requirement that one piece must move before neighboring pieces can move. One can think of this process as transferring the macroscale energy into the kinetic energy of the fragments. Since we do not track this form of energy, we allow it to be dissipated. The rate of plastic work due to granular flow (τ: d vp ) is computed using the difference between the strain energy in the trial (reached through elastic deformation) state and the state at the end of the time step. We have found that this is more accurate than computing the work rate from the plastic strain rate. The total rate of temperature change is the sum of the contributions from granular flow, isentropic heating, and viscous heating ̇ + θav ̇ + θGP ̇ ). (θ ̇ = θent

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

129

Table 1 Material model parameters for AlON. Thermal parameters

Equation of state parameters

Elastic parameters

Flaw distribution

Micromechanics parameters

Granular flow parameters

Pore compaction

Density (ρ) Specific heat capacity (Cv) Thermal conductivity (αθ )

3595 kg/m3 800 J/(kg K)

C0 S Γ0

7670 m/s 1.3 1.6

Shear modulus (G) Bulk modulus (κ)

125 GPa 211 GPa

Minimum flaw size (smin ) Maximum flaw size (smax ) Distribution exponent (α) Flaw density (η)

2 μm

9 W/(m K)

40 μm 3.0 4  1012 m  3

Fracture toughness (KIC) Rayleigh wave speed (Cr) Maximum crack velocity (Vm) Crack growth exponent (γc) Crack face coefficient of friction (μ) Crack orientation (ϕ)

2.5 MPa m 5.41 km/s 0.2 Cr

A Y Damage cohesive strength (B) Relaxation time (τGP ) Damage for granular flow (Dc) Maximum damage (Dmax )

0.6 0 MPa 0.1 MPa

Reference crush pressure (P0) Reference distension (JGP 0 ) Consolidation pressure (Pc)

1 GPa

1.0 0.57 60°

7  10  9 s 0.125 1.0

2.0 10 GPa

Pore compaction is the final piece of the material model description. In addition to pore compaction, the model also incorporates a Mie–Grüneisen equation of state, micromechanics based damage growth with a dynamically interacting distribution of flaws, degradation of the elastic behavior of the material with damage, and granular flow once a sufficient level of damage is reached. All of the material parameters are summarized in Table 1.

4. Flaw sampling and the coupling between specimen size and strength Coupling the flaw sampling discussion in Section 2 with the material model discussed in Section 3 results in one of the major motivations for retaining local information about the flaw distribution. When testing ceramic materials, significant variability in the strength is observed from one specimen to the next. This variability depends on the specimen size, type of loading, and loading rate. While Weibull based arguments reproduce these observed trends, the Weibull distribution is an extreme value distribution that assumes the outcome (strength) is controlled by a single event (the first activated flaw). These arguments have been extended to treat multiple flaws based on an “event horizon” approach (e.g Hild et al., 2003; Benz and Asphaug, 1994; Das and Cleary, 2010). Damage models based on Weibull arguments may be valid for the over driven systems discussed in Hild et al. (2003) (tensile fragmentation of glass within a fractal regime), but it is unclear how these models should extend to other loading regimes. For example, under dynamic compressive loading Paliwal et al. (2006)

130

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

Fig. 4. Weibull plot of the simulated compressive strength of a 0.5 mm cube at a strain rate of 103 1/s using 1000 samples. A Weibull distribution seems to fit the majority of the distribution of simulated strengths, but provides a poor fit for the low strengths.

observed multiple bright spots prior to failure indicating the activation of multiple flaws that initially grow slowly then accelerate as they interact. The approach discussed in this work provides a reasonable approximation for extension from two dimensions to three dimensions while computing the strength distribution and rate dependence based on the input flaw distribution. The statistical variation of specimen strength with specimen size for different spatial distributions of flaws using a similar micromechanics model and a Gaussian flaw size distribution was investigated in Graham-Brady (2010). The author concluded that, while a standard two parameter Weibull distribution reproduces the expected trends in the relationship between strength and specimen size, it does not provide a good fit to a simulated distribution of strengths in dynamic compression. We reach the same conclusion using a bounded Pareto distribution of flaws instead of a Gaussian distribution as shown in the Weibull probability plot given in Fig. 4. The Weibull distribution over-predicts the probability of low compressive strengths indicating that it is a conservative engineering tool when attempting to define a safe load for a structure. In this section we expand the study using our material model and sampling procedure to illustrate the coupling between flaw sampling, strain rate sensitivity, and the distribution of compressive strengths. Since we expect the strength of a specimen to depend on both size and loading rate, we simulate the uniaxial compression problem using a variety of specimen sizes from 4.0 mm on a side to 0.125 mm on a side and a strain rate of 103 1/s. To ensure a homogeneous stress state, these simulations were performed using a single particle with free boundaries in the e1 and e 2 directions, a constant velocity on the positive e3 surface and zero displacement on the negative e3 surface. The results are summarized in the box plots shown in Fig. 5. In general, large specimen sizes have lower median strengths, as indicated by the red horizontal bar, and a narrower distribution of strengths, as indicated by the box and whisker sizes. Larger specimens are expected to be weaker because larger specimens have a higher probability of containing a large flaw. However, one also expects larger specimens to have less variability than smaller specimens since the large specimens should contain more of the flaw distribution and, as a result, the sampling effect of the flaw distribution becomes less important.

Fig. 5. A box plot showing the change in the distribution of uniaxial compressive strengths as the simulated specimen edge length is decreased from 4 mm to 0.125 mm. The median strength and the variability in the strength increase as the specimen size decreases. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

131

Fig. 6. Box plot of the compressive strength variation with strain rate for a sample size of 0.5 mm per side. The median strength increases with strain rate while the spread in the strength distribution as indicated by the whisker length decreases with strain rate.

These trends are consistent with Weibull models, but we expect a different quantitative distribution than a Weibull model would predict. In general, for a fixed specimen size, as the strain rate increases the strength increases and the spread in the strength distribution decreases. This trend is shown in Fig. 6 for a specimen size of 0.5 mm per side. Note that the nature of the rate sensitivity depends on the local flaw distribution, because it is a competition between the stress required to drive the activated cracks faster and the activation of the next set of available flaws by the cracks. This complex interplay between local sampling, strain rate sensitivity, and the specimen size is one of the reasons that we elect to retain the local flaw distribution data despite the associated additional computational cost (discussed in Section 6). The combined influence of specimen size and strain rate on the strength distributions is shown in Fig. 7, which presents the empirical CDFs of the peak strength for two specimen sizes (0.5 mm/cell and 0.125 mm/cell) and three loading rates from ϵ̇ = 103 1/s to ϵ̇ = 105 1/s. In this plot the blue lines correspond to 0.5 mm specimens and the red lines correspond to the 0.125 mm specimens. The different strain rates are denoted by the line patterns (solid lines correspond to strain rates of 103 s  1, dashed lines to rates of 104 s  1, and dash-dotted lines to rates of 105 s  1). For all specimen sizes, the strength increases when the strain rate is increased. This effect is a result of the micromechanics damage model, and specifically the flaw distribution and limiting crack growth velocity. For all of the strain rates, as the specimen size decreases, the median strength increases and the variability in the strength also increases. The physical reason for these trends relates directly to the flaw sampling process discussed in Section 2. As the specimen size increases, the local flaw distribution approaches the parent flaw distribution. Since the local flaw distributions in larger specimens are in general closer to the parent flaw distribution, they are in general also closer to each other. This explains the trend towards greater variability in smaller specimens. The mean strength increases with decreasing specimen size because smaller specimens are less likely to have large flaws, which results in an increase in strength. Even though a large portion of the flaw distribution participates in the damage growth process, the initial damage growth is determined by the largest flaw.

Fig. 7. The distribution of material point strengths is a function of both material point size and the loading rate. As the specimen size decreases from 0.5 mm/cell (blue) to 0.125 mm/cell (red) the median strength and the variability in the strength increases (Fig. 5). The strain rate sensitivity, shown by the difference between the solid and dash-dotted lines, decreases with smaller specimen sizes. These changes are a result of the flaw distribution sampling and are a key benefit of using the micromechanics based damage model. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

132

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

The final observation relating to these strength distributions is that the rate sensitivity (strength increase between a strain rate of 103 s  1 and 105 s  1) decreases as the sample size decreases. This is caused by both the sampling and the micromechanics damage model as previously noted. The high rates require more active flaws to relax the stresses, and the additional flaws that are available are smaller and require high stresses for activation. When the specimen size is reduced, most of the realizations will be missing the large flaws in the flaw distribution. Since the large flaws are missing, higher stresses are required to initiate damage growth, resulting in an increased strength. The flaw sampling procedure discussed in Section 2 enforces integer numbers of flaws within each discretization volume, which results in a quantization of the flaw densities. In smaller specimens adding or removing a single flaw has a large effect on the local flaw density. This results in a reduced strain rate sensitivity because the strengths at lower strain rates are increased more than the strengths at higher strain rates. Although the mechanisms behind these size and rate effects are general, the degree of rate sensitivity and size effect depend on the specific defect distribution. We expect that the trends would be qualitatively similar but quantitatively different if a different flaw size distribution was used. The power of this approach is that the rate and size effects are natural outcomes of the input distribution, and do not need to be specified a priori. In the following section we use this framework and model to perform simulations of Edge on Impact (EOI) experiments (Strassburger et al., 2006) conducted on AlON tiles.

5. Simulations of edge on impact experiments The EOI experiment is interesting because it provides real time information about the dynamic propagation of damage as a result of an impact event. In the experiment (Fig. 8), a projectile strikes the edge of an armor ceramic tile. As a result of the impact, stress waves travel through the target and cause the formation of a network of cracks. Since the tile is thin, and in the case of AlON transparent, the development and propagation of the crack network can be recorded using high speed cameras. In the particular experiment that we simulate (labeled experiment number 14897 in Strassburger et al. (2005)), a 10 mm thick, 100 mm square AlON tile is impacted at a velocity of 381 m/s by a 23 mm long steel cylinder with a diameter of 30 mm. We model the steel projectile as a simple elastic plastic material with linear strain hardening, and whose properties are summarized in Table 2. These material properties represent a high strength maraging steel; however, we have preformed simulations using a yield strength of 180 MPa that represents SAE 1010 steel in a hot rolled condition Norton, 2000, Table C-9, and observed a similar results for the damage pattern. The impactor strikes the edge of the tile causing a damage front to develop and propagate through the tile. The impact event and subsequent damage front propagation were imaged using high speed photography (Strassburger et al., 2006). Images were captured every 0.5–1.0 μs (Strassburger et al., 2005).

Fig. 8. Comparison of experimental results with simulation results for 3 and 6 μs after impact with a resolution of 0.25 mm/cell. The upper row of images shows a normal view of the of simulation with the top half of the view showing the center plane of the plate and the bottom half showing the surface of the plate. The color scale indicates the damage level. Blue is undamaged, red indicates damage that has reached the granular flow threshold, and black indicates the maximum allowable damage. The bottom row of images show perspective views of the simulations and the corresponding experimental images. Damage develops first on the interior of the plate and then arrives at the surface. Non-homogeneous damage patterns develop at later times. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

133

Table 2 Material properties for projectile. Bulk modulus Shear modulus Yield stress Hardening modulus Density

173 GPa 80 GPa 2000 MPa 750 MPa 7.83  103 kg/m3

Using these images (Strassburger, 2004) computed the velocity of the damage front resulting from the impact event. The authors also computed the velocity of discrete cracks that nucleated in front of the primary damage front. Since this material model is not designed to track the evolution of individual cracks, we compare our simulation results to the reported damage front location. 5.1. Simulation setup These simulations were performed using convected particle domain interpolation (CPDI) (Sadeghirad et al., 2011) and explicit time integration (Wallstedt and Guilkey, 2008) implemented in the Uintah (Parker et al., 2011) computational framework. CPDI is based on the Material Point Method (Sulsky et al., 1994) in which the constitutive description of the material is carried in a set of Lagrangian material points. A background Eulerian grid is used to compute gradients and solve the equations of motion. The stability of some implementations of the Material Point Method are discussed in Love and Sulsky (2006). In these simulations we use 1 particle per cell in the grid. This choice focuses the computational effort on resolving high stress gradients while sacrificing the ability to resolve sub grid variations in the constitutive response. We compare our simulations to experimental shadowgraphs that show the damage pattern at discrete instants in time. Examples of the experimental images are the bottom two images in Fig. 8, taken 3.0 μs and 6.0 μs after impact. These figures show light and dark regions. In the light regions, the light is fully transmitted, and we interpret this as a state of no damage. In the dark regions, the light is blocked by a large amount of damaged material through the thickness of the tile. We interpret the gray regions in the experimental images as regions where there is some cracking or damage in some location through the thickness of the tile, but the damage is less extensive than in the fully black regions. In the simulation results we present an isometric view of the damage pattern where a quarter of the specimen is removed. This view allows us to see damage on both the surface and on the mid plane of the tile. At both 3.0 and 6.0 μs, the computed damage extent is larger at the mid-plane of the tile than at the surface. The damaged region is concentrated in a zone downrange of the projectile. The material is fully damaged near the projectile. As one moves further from the projectile and closer to the surface the damage zone becomes more heterogeneous. This is especially evident in the 6.0 μs image. Although the simulations do not show well developed fingers of damaged regions as seen in the experiments, they do show a heterogeneous damage pattern near the leading edge of the damage front. It is always important to understand the effect of mesh resolution on the results and ensure that the mesh has sufficient resolution. We conducted a mesh refinement study (Fig. 9) using four different meshes (1.0–0.125 mm/cell) and examining the time history of the damage extent at the specimen surface and the center plane of the specimen measured along the impact direction (Fig. 8). We compute the damage extent by extracting the damage value for all points along the specified path then finding the maximum distance along the path where there is a damage value grater than 0.124. The threshold of

Fig. 9. Mesh refinement study results. The damage extent at the center plane of the target shows good agreement at all resolutions. Damage reaches the top surface of the target at earlier times with a finer resolutions, but the change in the arrival times is decreasing with increased resolution. (a) Damage extent at the center plane of the tile along the impact direction. (b) Damage extent at the top surface of the tile along the impact direction.

134

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

0.124 is used to identify locations where the damage is sufficiently large that granular flow is possible. The damage growth rate is very rapid when damage reaches 0.124 so it is unlikely that damage will reach 0.124 but not continue growing to reach 0.125 and enable granular flow. The location history of the damage extent at the mid plane of the tile agrees (Fig. 9(a)) for all times for resolutions larger than 0.5 mm/cell. The damage extent at the top surface is more sensitive to the cell size (Fig. 9(b)), but the 0.25 mm/cell and the 0.125 mm/cell results are similar. Damage growth on the surface of the tile is sensitive to interactions through the thickness of the tile, which are discussed in Section 5.3. The different convergence behavior on the surface of the tile and at the center plane leads us to believe that the through thickness wave interactions in the tile are under resolved. Since the 0.125 mm/cell simulations took almost 60 h on 512 processors of the DoD Copper machine, and produced an output data file that used 2.0 TB of disk space, we did not pursue a higher resolution calculation to see if the 0.125 mm/cell calculation was converged. Simulations of strong shocks (Tonge, 2014) using the artificial viscosity parameters used in this work showed that the shock front was spread over up to 9 computational cells. Since the plate is relativity thin compared to its other two dimensions, we believe that the artificial viscosity could have a greater effect on wave propagation through the thickness of the plate than in the plane of the plate. At the 0.5 mm/cell resolution the wave front would be spread from the center of the plate to the surface due to artificial viscosity alone. We recognize this as a limitation and are careful in the discussion in Section 5.3 to present results that we believe are consistent with expected wave mechanics solutions and not overly sensitive to the through thickness resolution. To balance computational effort and our ability to look at multiple sets of material parameters, we use the 0.25 mm/cell resolution for all subsequent results. Strassburger et al. (2005) computed the distance from the impact site to the boundary between the light and dark regions in the experimental images as a function of time. In Fig. 10 we plot their results (in black) along with the computed damage extent measured at the center plane and on the surface of the target. The experimentally observed damage front velocity was 8381 m/s (Strassburger, 2004), which is 89 percent of the 9367 m/s experimentally measured wave speed using crossed polarizers. In the simulations the longitudinal (uniaxial strain) wave speed (based on the material properties listed in Table 1) is 10,256 m/s and the computed damage velocity in the center of the tile is 9900 m/s, which is 97 percent of the longitudinal wave speed. These high damage velocities suggest that the damage nucleation and growth is driven by the longitudinal wave. The simulated damage front velocity seems to agree better with the crack nucleation front reported in experiments than the reported damage front velocity. The nucleation or activation of defects in the experiments could correlate with damage growth because the damage model tracks the activation and growth of pre-existing flaws. However a definitive conclusion separating the effects of crack growth, nucleation, and granular flow on the measured shadowgraph would require a simulation methodology that both tracks the free surfaces created by cracks and simulates the interaction of light with these free surfaces. Such a simulation capability is not currently available, but in the near future such comparisons should be possible. The computed surface damage lags behind the damage at the center of the plate suggesting that the arrival of damage at the surface depends on the behavior of the damaged material and on how the damage in the center of the tile interacts with the free surface boundary condition. The general agreement with experiments provides confidence that our material model and simulation approach capture the important physical processes within this loading regime. Similar comparisons have been made in Leavy et al. (2013). Using the Kayenta material model the authors observed some localization patterns that resembled the fingering that was experimentally observed. Using our model we observe heterogeneous damage patterns but we have not performed sufficient analysis on the nature of these heterogeneous regions to conclusively equate them to the fingers that were observed in the experiments. 5.2. The consequences of variability and the damage kinetics Strassburger et al. (2005) also observed the development and growth of a number of finger-like localized damage zones during the experiments. This transition from a relatively homogeneous damage zone to the propagation of discrete localized features is common in experimental observations of the failure of brittle materials. Our simulations also show heterogeneity,

Fig. 10. Comparison between simulated and experimental observations of the damage front propagation history. The damage front velocity, based on a linear best fit, is reported in parenthesis in the legend. The experimentally observed damage front velocity is between the simulated damage growth velocity measured at the top surface and in the center of the plate.

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

135

Fig. 11. Damage pattern at the center of the plate 6 μs after impact demonstrating the effect of variability. Macroscopic variability promotes the development of a heterogeneous damage pattern, which is observed in the experimental results. The development of heterogeneous damage patterns is a prerequisite for developing physically reasonable fragment distributions.

but not the localization into fingers. Riedel et al. (2010) demonstrated that the ability of a model to capture the localization of damage into discrete fingers depends on both the material model and the treatment of localization by the computational framework. The development of a heterogeneous damage field is necessary prior to developing the finger pattern, but the current implementation of CPDI within the Uintah framework may not provide sufficient support for weak discontinuities to support the formation of the observed fingers. The Kayenta simulations referenced in the previous section were run in a Lagrangian mode using the advanced finite element code Alegra. Traditional 8 node hexahedral elements naturally admit a weak discontinuity in displacement over element boundaries. In our material model we illustrate the importance of local fluctuations in the material strength in Fig. 11 by disabling the sampling procedure discussed in Section 2. As shown in Fig. 11, the heterogeneous damage pattern disappears when the variability is removed. Without variability, the damage pattern is purely the result of the stress interactions, and since this system is highly symmetric, that symmetry is carried into the damage pattern and it does not have any localized features. While the variability has a strong effect on the heterogeneous nature of the damage, it seems to have little effect on the general kinetics and extent of the damage zone. From this observation we conclude that the damage growth rate and general damage shape are controlled primarily by the macroscopic loading, boundary conditions, and average material behavior. This is a reasonable result, because one expects to see generally the same pattern of damage for two similar tiles; however one does not expect the detailed local damage pattern to be identical for two different tests. This comparison provides another example of when the symmetry breaking effect of local variability is important for obtaining physically realistic results from simulations of high rate brittle failure. In addition to investigating the effect of flaw sampling, we are able to look at the effect of changing parameters that define the material behavior. During the discussion of the material model in Section 3 we identified sources for many of the model parameters. However there is still uncertainty about the effective friction coefficient (A) in the granular flow relationship as well as the maximum crack growth speed (vm) for the microcracks in the micromechanics based damage model. We investigated the effect of changing the pressure sensitivity of the granular flow and changing the rate sensitivity of the material by reducing the maximum crack velocity (vm). The pressure sensitivity of the granular material (A) has little effect on the damage front location in the center plane of the tile. This is illustrated by comparing Fig. 11a and the center simulation in Fig. 12. Both of these simulations use a maximum crack velocity (vm) of 20 percent of the Rayleigh wave speed, but use different effective coefficients of friction in the granular flow rule. Similar results were observed with friction coefficients of 0.2 and 1.2. This reinforces the conclusion that the damage front location in the center of the plate is dominated by interactions resulting from the longitudinal wave, which are discussed in more detail in Section 5.3. The

Fig. 12. Damage distribution at the center of the plate after 6 μs for different damage growth rates (vm). As the maximum crack velocity increases from 1 percent of the Rayleigh wave speed (Cr) to the Rayleigh wave speed the damage pattern behind the damage front becomes more heterogeneous. The granular flow slope (A) is 0.8 in these images, which is greater than the value of 0.6 used in the previous figures. The location of the damage front is the same for 100 percent and 20 percent of the Rayleigh wave speed.

136

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

friction coefficient of 0.8 does produce some black regions that have a finger-like shape near the impact point behind the main damage front. The black color indicates that the damage level within the material has continued to grow after the activation of granular flow. In this case the continued damage growth occurs because the coefficient of friction for granular flow is greater than the crack face coefficient of friction. Therefore, under moderate pressure there is a range of deviatoric stresses that are sufficiently large to cause crack growth, but are insufficient to cause granular flow. Since damage growth is a softening process (especially at large damage levels), damage above 0.125 tends to localize quickly if it is not inhibited by the granular flow mechanism. This produces the black finger-like patterns in Fig. 12. Since this pattern is bordered by fully damaged material we do not believe that it represents the same physical process that is observed in the Strassburger experiments. In the color scheme used in Fig. 12 this would be indicated by red fingers surrounded by blue undamaged regions. The damage kinetics have a moderate effect on the location of the damage front in the center plane of the tile, but have a large effect on the damage pattern behind the damage front. In Fig. 12 we demonstrate this effect by showing the damage pattern on the center plane 6 μs after impact for three simulations where the maximum allowable crack velocity increases from 1 percent of the Rayleigh wave speed, on the left, to 20 percent in the center, and 100 percent on the right. We alter the damage kinetics by changing vm because this quantity directly scales the rate of damage growth for a given flaw size distribution, wing crack length, and applied load. Changing the crack density would change the damage growth rate, but also changes the effective variability in the problem. By changing only the maximum crack velocity, we are able to isolate the effect of the damage kinetics. All of these simulations use a coefficient of friction of 0.8 instead of the 0.6 used in the baseline simulations. From these images, we see that as the damage kinetics become faster, we see larger damage gradients behind the damage front. Additionally, the shape of the damage front changes from a smooth, almost circular front with vm = 0.01Cr to an angular almost trapezoidal front with vm ¼Cr. The rate sensitivity in the damage model has a regularizing effect and resists the formation of sharp damage gradients (a longer time to failure provides more time for neighboring material points to develop similar damage). Decreasing the maximum allowable crack velocity increases the rate sensitivity of the damage resulting in a more uniform damage pattern. Conversely, increasing the limiting crack growth velocity decreases the rate sensitivity of damage and promotes localization and the formation of sharp damage gradients. While the damage extent at the center plane is relatively insensitive to the granular flow slope, the damage extent on the surface is sensitive to the granular flow slope. The location of the damage extent on the surface of the tile changes when the pressure dependence of granular flow is altered. Comparing the computed arrival time of the damage front (Fig. 13) to that seen in the experiments, we find that A should be smaller than 0.8 and use A ¼0.6 in all subsequent calculations. 5.3. Center plane damage driven by the longitudinal wave The experimental images in Strassburger et al. (2005) present a 2 dimensional view of the damage propagation within the target plate; however, the development of damage is a three dimensional process. As discussed in the previous section, the location of the damage front at the center plane of the tile is in front of the damage front observed on the top surface of the tile. In this section we investigate the stress wave interactions that favor damage growth on the interior of the tile. Both the experimental measures of the damage velocity and the simulations indicate that the damage front moves at around 90 percent of the longitudinal (uniaxial strain) wave speed (Fig. 10). To explain the reasons for this behavior, we investigate the stress and damage pattern on a cross section through the target plate 6 μs after the impact event in Fig. 14. The schematic at the top of Fig. 14 shows the orientation of the cross section with respect to the tile and the impactor.

Fig. 13. The damage front location on the surface of the plate suggests that a granular flow coefficient (A) of 1.2 is too high. Although the granular flow coefficients of 0.2 and 0.8 produce different damage front location histories, there is not sufficient experimental data to favor one of these two over the other.

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

137

Fig. 14. Cross sections from the simulation after 6 μs illustrate nonuniform behavior in the through thickness direction due to wave interactions and inertial confinement. Release waves from the longitudinal loading wave interact to initiate damage below the surface of the tile. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

Looking at the pressure (top cross section image) and the equivalent stress (second cross section image) reveals that the longitudinal wave (shown in the dashed green line) reflects from the free surfaces and causes a region of tension when the reflected longitudinal waves from the two boundaries interact. The domains of the reflected longitudinal and shear waves are shown in the equivalent stress plot using black and red lines respectively. Comparing the damage pattern (third image down) with the pressure and equivalent stress indicates that the damage begins to grow first in the high shear region just under the surface of the tile. As the damage develops and granular flow is activated, the local stresses are relieved. The cross

138

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

section views of the deviatoric (fourth image down) and volumetric (bottom image) components of granular flow indicate that, for a fixed distance from the impact site, there is more granular flow away from the center of the tile than at the center. This occurs because of the effective inertial confinement at the center of the plate. The dilation associated with granular flow requires that the material on the surface of the plate moves away from the center line to accommodate the additional porosity. Accelerating the material away from the center plane of the plate provides a confining stress that resists granular flow.

6. Relationship to other similar ceramic failure models There are several material models available for modeling the response of ceramic (or geologic) materials subjected to dynamic loading. Some of the most widely used models are the ones developed by Johnson and Holmquist (1993), Johnson et al. (2003) and Holmquist and Johnson (2006). These models are designed for engineering applications with goals of robustness, speed, and the ability to fit existing empirical data. The generalized plasticity framework (Brannon et al., 2009) was also designed for engineering applications with input from microphysical theories, but with the explicit goal of fitting triaxial compression data for geological materials. Since these models lack a connection to microstructural features, it can be difficult to evaluate material design compromises using them. The multiplane microcracking model (Gailly and Espinosa, 2002) and the model developed by Dienes et al. (2006) are micromechanics based models that use energy based arguments to determine the microscale crack growth behavior. These models can track multiple crack orientations (something our proposed model does not do). The model developed in this work and the model by Deshpande et al. (2011) use the wing cracking approach developed in Ashby and Hallam (1986). The model in this work and the model developed by Deshpande and coworkers are very similar. The primary difference between these two models is that the new model provides for material variability and should capture Aleatory uncertainty (Brannon et al., 2007), while the Deshpande model does not. There are minor differences regarding the transition from damage growth to granular flow and the model for crack interactions. The new model is developed and qualitatively validated using a three dimensional implementation while the Deshpande model is limited to axisymmetric problems. As a result, the Deshpande model cannot capture the radial cracking that occurs during an impact event. The models developed by Guy et al. (2012) and Hild et al. (2003) are based on an “obscuration zone” hypothesis and a Weibull distribution of flaw activation stresses in uniaxial tension. These models are very attractive for tension dominated problems, especially the fractal fragmentation arguments (for severely over-driven problems) and the local to non-local transition that the authors developed the models to solve. Micromechanics based models use many material parameters that are directly measured using standard experimental techniques, but the remaining inputs such as the flaw distribution used in the model can be difficult to estimate without detailed microstructural analysis. The rapid development of advanced material characterization techniques such as 3D X-ray computed tomography should enable quantitate defect mapping at the required resolution in a non-destructive environment in the near future. Even with such an analysis, it is currently difficult to map real microstructural features, like a weak grain boundary, to idealized slit flaws. The material model discussed in this work was not extensively calibrated using all available experimental data, but instead was fit to dynamic uniaxial compression testing. The Kayenta model is more straightforward to calibrate from experimental measurements of the peak stress at various confining pressures. For the work on ALON (Leavy et al., 2013) the Kayenta model was calibrated using meso-scale simulations. The calibrated model was then used to simulate the EOI experiment using the advanced finite element program Alegra. Within Kayenta, the time that it takes a material volume (material point or integration point) to transition from intact to fully damaged is a model input based on the shear wave speed of the material. A similar timescale in the proposed model depends on the simulated number of flaws within the sub volume and the magnitude by which the stress exceeds the flaw activation stress. This is one difference in the material response that could contribute to the different observed fingering behavior. The rate effects in the proposed model should not be sufficient to prevent localization of damage into fingers (part II of this series shows sphere on cylinder simulations where radial localizations that look like cracks are observed). As discussed by Riedel et al. (2010) in the context of Smoothed Particle Hydrodynamics (SPH), the particular choice of how a numerical scheme represents localization can have a large effect on the localization pattern. The Uintah implementation of CPDI does not support the formation of a weak displacement discontinuity (Sadeghirad et al., 2013) while linear displacement based Lagrangian finite elements (which were used in the Alegra simulations) will support a weak displacement discontinuity. This is likely a contributing factor to the differences between the localizations observed using the proposed model and those using Kayenta. The first author is working to implement the proposed model in the same computational frameworks where Kayenta and the models by Johnson and coworkers have been implemented. When the models are implemented in the same host code they can be more effectively compared. Micromechanics based damage models introduce additional physics into failure modeling of materials, but also introduce an additional computational cost. In an attempt to quantify the computational cost of the proposed model, we have performed a timing study that compares the execution time using the Kayenta model and the proposed model with different numbers of flaw bins. For these tests Uintah was run using a single core of an 8-core Intel Sandy Bridge processor. Uintah was compiled with the Intel compiler suite with the optimization flags set to -g -O3 -axvex. The timing results were captured using Intel's vTune utility for a problem containing 174,176 particle domains. The results are summarized in Table 3. The cost of the proposed model depends heavily on the number of bins used to represent the local flaw population.

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

139

Table 3 Timing results for an optimized build of the Uintah framework using the proposed model and Kayenta to illustrate the computational cost of the proposed model. The reported numbers are the average time in seconds spent in each phase of the calculation per time-step for a simulation containing 174,176 particle domains. Program phase

Kayenta

1 bin

25 bins

100 bins

Total time Constitutive model time Particle loop time Damage calculation Self consistent calculation Crack growth

3.867 0.9813 N/A N/A N/A N/A

3.506 0.5141 0.5078 0.2688 0.2281 0.0086

3.888 0.9094 0.8734 0.4391 0.2109 0.1969

5.550 2.005 1.945 0.9656 0.2203 0.7203

Using 25 bins results in a material model that is moderately expensive, but can still be used for large scale problems.

7. Summary and future work We present a micromechanics based material model that incorporates flaw sampling statistics, damage growth from a distribution of interacting microcracks, a Mie–Grüneisen equation of state, and granular flow of the fully damaged material. A major feature of the model is the explicit incorporation of flaw sampling statistics by considering a Poisson process acting at each material point sub volume. This sampling process results in a distribution of effective material point strengths, which depends on both size and loading rate. We simulate an Edge on Impact experiment (Strassburger et al., 2005) to partially and qualitatively validate the model. We then use the model to understand why damage grows preferentially in the interior of the plate rather than on the surface in that experiment, and discuss the coupling of propagation of the damage zone with crack kinetics and the parameters for the granular flow model. The combination of this material model with a numerical technique such as CPDI allows us to perform detailed simulations of impact events where the failure process is linked to microstructural variables. The microstructural link is important for material design. In future work, it will be possible to use experimentally obtained microstructures to provide a flaw distribution, leading to the possibility of predictive material models. Similarly, since the microstructure is accounted for, we could perform a study using different flaw distributions to suggest promising investigation paths towards improved material performance.

Acknowledgment The effort of A.L. Tonge was partially supported by the National Science Foundation through IGERT 0801471. Part of this 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 authors also acknowledge many useful discussions with B. Leavy of the Army Research Laboratory regarding Kayenta and modeling of the Edge on Impact Experiments. The authors would also like to thank the Uintah development team for providing the Uintah framework as an open source project primarily developed by the University of Utah. The authors received particular help from J. Guilkey, J. Schmidt, and A. Humphrey regarding the use, installation, and organizational principals governing the Uintah framework.

Appendix A. Mie–Grüneisen EOS details We start by assuming that the internal energy (e) associated with volume changes (Je) can be separated into a thermal contribution eθ and a “cold” contribution ec, which is only a function of the volumetric deformation Je:

e (J e , θ ) = ec (J e ) + eθ (J e , θ ).

(A.1)

Following the arguments in Section 4.2 of Drumheller (1998) we arrive at the Grüneisen equation which relates the total pressure p to the “cold” reference pressure pc and the temperature θ through the use of the Gruneisen parameter Γ and the specific heat at constant entropy cη :

p = pc + ρ0

Γ cη θ Je

(A.2)

140

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

This equation could also be written in terms of the internal energy and the cold internal energy:

p = pc + ρ0

Γ ( e − ec ) Je

(A.3)

A suitable reference curve is the Principal Hugoniot, which is the locus of points that are achievable through a single shock process from ambient conditions (room temperature and pressure). The balance of mass, momentum, and energy across the shock front are referred to as the Rankine–Hugoniot shock jump conditions (Drumheller, 1998, Section 3.4)

ρ0 Us = ρ ( Us − Up )

(A.4)

p − p0 = ρ0 Us Up

(A.5)

1

pUp = 2 Us Up2 + ρ0 Us (e − e0 )

(A.6)

Eqs. (A.4) and (A.5) can be combined with the volumetric deformation (Je) to write the pressure on the Hugoniot as a function of Je and the empirical relationship between the shock speed Us and the volumetric deformation: 2

pH (J e ) = ρ0 ( Us (J e ) ) (1 − J e )

(A.7)

These equations (A.4)–(A.6) can be rearranged to write the internal energy as a function of the pressure (pH (J e )) and the volumetric deformation:

eH (J e ) =

pH (J e ) ( 1 − J e ) + e0 2ρ 0

(A.8)

Here e0 is the internal energy at ambiant conditions. We now apply Eq. (A.3) to compute the equilibrium pressure (p+ ) at the Hugoniot state:

p+ = pc (J e ) + ρ0

Γ (eH − e0 − ec (J e ) + e0 ) = pH Je

(A.9)

Solving for the cold pressure (pc) and using Eq. (A.8) results in:

⎡ ⎤ Γ pc (J e ) = pH (J e ) ⎢ 1 − e (1 − J e ) ⎥ + ρΓ (ec (J e ) − e0 ). ⎣ ⎦ 2J

(A.10)

Using Eq. (A.10) to evaluate Eq. (A.3) completes the transition from a reference curve as the cold curve to a reference to the principal Hugoniot:

⎛ ⎞⎞ Γ (J e ) ⎛ Γ (J e ) p (J e , θ ) = pH (J e ) ⎜ 1 − ⎜ 1 − J e ⎟ ⎟ + ρ0 e (e (J e , θ ) − e0 ). e 2J ⎝ J ⎝ ⎠⎠

(A.11)

Using Eq. (A.1) and the specific heat at constant entropy cη we rewrite the internal energy using the temperature, the cold energy, and the specific heat:

e (J e , θ ) = ec (J e ) + c η (η) θ .

(A.12)

We define the cold energy such that (ec (1) = 0). Therefore the reference energy e0 is given by

e0 = c η θ 0

(A.13)

Using Eqs. (A.12) and (A.13) we rewrite the pressure p (J e , θ ) using the temperature and cold energy:

⎛ ⎞⎞ Γ (J e ) ⎛ Γ (J e ) p (J e , θ ) = pH (J e ) ⎜ 1 − ⎜ 1 − J e ⎟ ⎟ + ρ0 e (ec (J e ) + c η (θ − θ 0 )) 2J e ⎝ J ⎝ ⎠⎠

(A.14)

The cold energy, by definition, is related to the cold pressure:

− ρ0

dec = pc . dJ e

(A.15)

Using this equation and Eq. (A.10) we have a first order linear ordinary differential equation for the cold energy:

⎞⎞ −pH ⎛ dec Γ Γ ⎛ + e (ec − e0 ) = ⎜ 1 − e ⎜ 1 − Je ⎟⎟ e dJ J 2J ⎝ ρ0 ⎝ ⎠⎠ This equation can be integrated once the reference pressure pH (J e ) and the Grüneisen parameter Γ (J e ) are specified.

(A.16)

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

141

We note that the specific heat is only a function of the entropy (cη (η)) and the Grüneisen parameter is only a function of the density (Γ (ρ)). The specific heat is a measure of the coupling between energy and temperature. The Grüneisen parameter quantifies the coupling between thermal energy (atomic vibrations) and pressure. The cold energy is the potential energy stored by the atomic bonds in the absence of atomic vibrations. This is the same measure as the strain energy except strain energy is referenced per unit volume not unit mass. In order to finish specifying the equation of state, we select functional dependences for Γ (J e ) and cη (η):

Γ (J e ) = J e Γ0

(A.17)

c η (η) = c η = cv

(A.18)

Here Γ0 and cv are empirical material constants. The choice of Γ (J e ) = J e Γ0 implies that Γρ = Γ0 ρ0 is a constant. This choice of Γ allow us to simplify Eq. (A.14) to

⎛ ⎞⎞ Γ ⎛ p = pH (J e ) ⎜ 1 − 0 ⎜ 1 − J e ⎟ ⎟ + ρ0 Γ0 ( ec + c η (θ − θ 0 ) ). ⎝ ⎠⎠ 2⎝

(A.19)

Many materials are accurately represented with a linear relationship between the shock speed Us and the particle velocity Up. If a linear relationship is assumed, the shock velocity can be written as

Us = C0 + SUp Us (J e ) =

(A.20)

C0 1 − S (1 − J e )

(A.21)

Shocks do not form under hydrostatic tension ( J e > 1.0) in most materials. Therefore, we assume a constant wave speed (Us = C0 ) in this regime. This specification of Us (J e ) results in a reference pressure curve given by

ρ0 C02 (1 − J e ) 2

pH (J e ) = { ( 1 − S (1 − J e ) ) ρ0 C02 (1

−

, J e < 1.0

Je )

.

otherwise

(A.22)

Using Eq. (A.17), the differential equation for the cold energy (Eq. (A.16)) becomes

−pH (J e ) ⎛ ⎞⎞ Γ0 ⎛⎜ dec ⎜1 − + Γ0 (ec − e0 ) = 1 − Je ⎟⎟ e ⎠⎠ ρ0 ⎝ dJ 2⎝

(A.23)

This equation can be solved by the method of integrating factors:

ec (J e ) = exp (Γ0 (1 − J e )) ( cv ( θ H (J e ) − θ 0 )) + cv θ 0.

(A.24)

Rearranging this equation into a more convenient form:

ec (J e ) = exp (Γ0 (1 − J e )) cv θ H (J e ) + cv θ 0 ( 1 − exp (Γ0 (1 − J e )) ) We have introduced the term θH

⎛ cv θ H (J e ) = ⎜ − ⎝

∫1

Je

(J e ) ,

(A.25)

which is an effective temperature given by:

⎞ ⎞ pH (J ) ⎛ Γ ⎜ 1 − 0 (1 − J ) ⎟ exp ( − Γ0 (1 − J )) dJ ⎟. ρ0 ⎝ 2 ⎠ ⎠

(A.26)

By introducing the compressive volume strain ϵc = 1 − J e with dJ′ = − d ϵ′c we write this equation in the more compact form:

cv θ H (J e ) =

∫0

ϵc

⎡p ⎛ ⎞⎤ 1 exp ( − Γ0 ϵ′c ) ⎢ H ⎜ 1 − Γ0 ϵ′c ⎟ ⎥ dϵ′c 2 ⎣ ρ0 ⎝ ⎠⎦

(A.27)

Under compression this integral must be evaluated numerically, but under tension (ϵc < 0) the cold energy (ec) is

ec (ϵc ≤ 0) = C02 ϵ2c + cv θ 0 (1 − exp (Γ0 ϵc ))

(A.28)

Once the cold energy is known the pressure at a given temperature and volumetric compression is given by Eqs. (A.19) and (A.12):

⎡ ⎞⎤ Γ ⎛ p (J e , θ ) = pH (J e ) ⎢ 1 − 0 ⎜ 1 − J e ⎟ ⎥ + ρ0 Γ0 ⎡⎣ ec (J e ) + c η (θ − θ 0 ) ⎤⎦ ⎣ ⎠⎦ 2⎝

(A.29)

142

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

Appendix B. Calculation of the isotropic strain energy density function In this section we derive the isotropic strain energy density function associated with the general strain energy density defined in Eq. (19). The isotropic strain energy density is defined by ρ (n ) = 1. Making this substitution and splitting the integral into pieces results in:

f=

D D 1 τ: 0 : τ + Z r τ · τ : n ⊗ n dω − 2 8π ω 8π D D dω + Z c τ : I ⊗ n ⊗ n : τ dω − 8π ω 8π

∬

∬

∬ω (Zr − Zn ) τ: n ⊗ n ⊗ n ⊗ n: τ dω + 8Dπ ∬ω Zc τ: n ⊗ n ⊗ I: τ

∬ω 2Zc τ: n ⊗ n ⊗ n ⊗ n: τ dω

(B.1)

In these integrals the only dependence on orientation is in n . By rearranging the integrals to move all of the orientation independent terms outside of the integral we have:

f=

⎛ ⎞ ⎛ ⎞ 1 D D D τ: 0 : τ + (Z − Z n ) τ: ⎜ Zr τ ·τ: ⎜ n ⊗ n dω⎟ − n ⊗ n ⊗ n ⊗ n dω⎟: τ + Zc τ ⎝ ω ⎠ 8π r ⎝ ω ⎠ 2 8π 8π ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ D D :⎜ 2Z c τ : ⎜ n ⊗ n dω⎟ ⊗ I: τ + Zc τ : I ⊗ ⎜ n ⊗ n dω⎟: τ − n ⊗ n ⊗ n ⊗ n dω⎟: τ ⎝ ω ⎠ ⎝ ω ⎠ ⎝ ω ⎠ 8π 8π

∬

∬

∬

∬

∬

(B.2)

There are two integrals that we need to evaluate. The approach for evaluating these intergals is due to conversations with R. Brannon. We start by evaluating

B1 =

1 4π

∬ω n ⊗ n dω

(B.3)

We can interpret n as a position vector that extends from the origin to the boundary of the unit sphere:

B1 =

1 4π

∬∂Ω x ⊗ n ds.

(B.4)

Here x is a position vector in the unit sphere (Ω) and n is a normal vector to the sphere surface. From the Gauss divergence theorem we convert the surface integral into a volume integral:

B1 =

1 4π

∭Ω ∇x x dV .

(B.5)

The gradient of a position vector is I . The integral reduces to computing the volume of the sphere

B1 =

1 IV (r ). 4π

(B.6)

In this integral the volume is the unit sphere and therefore the total integral is given by:

B1 =

1 3

I.

(B.7)

A similar procedure can be used to evaluate the other integral tensor:

2 =

1 4π

∬ω n ⊗ n ⊗ n ⊗ n dω.

(B.8)

Writing this tensor in terms of its components:

( B2 )ijkl = 41π ∬ω ni nj nk nl dω. We interpret the first three normal vectors as position vectors in the unit sphere

( B2 )ijkl = 41π ∬∂Ω xi xj xk nl dS.

(B.9)

Ω: (B.10)

Applying the generalized divergence theorem:

( B2 )ijkl = 41π ∭Ω

∂(xi xj xk ) dV . ∂xl

(B.11)

Expanding using the derivative:

( B2 )ijkl = 41π ∭Ω xi, l xj xk + xj, l xi xk + xk, l xi xj dV

(B.12)

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

143

We represent one of the position vectors in each of the terms in the sum as an orientation n multiplied by a radius r and write the volume integral as a surface integral and a radial integral:

( B2 )ijkl = 41π ∫0

1⎛

⎜ ⎝

∬∂Ω

r

⎞ δil xj nk + δjl xi nk + δkl xi nj dS⎟ r dr . ⎠

(B.13)

The evaluation of B1 indicates that the surface integral evaluates to:

∬∂Ω

δil xj nk dS = r

1 4π r 3 δil δjk 3 3

(B.14)

Comparing this result to Eq. (B.13) results in

( B2 )ijkl = 13 ( δil δjk + δjl δik + δkl δij ) ∫0

1

r 4 dr .

(B.15)

These three terms are the three isotropic fourth order tensors ( , ¯ , I ⊗ I ) and the integral evaluates to

( B2 )ijkl = 151 (  + ¯ + I ⊗ I ) = 151 I ⊗ I + 152 

1 : 5

(B.16)

Appendix C. Summary of Paliwal–Ramesh micromechanics damage model The self consistent problem can be solved analytically, as discussed in Hu et al. (2015, Appendix B) and provided for reference in Appendix D, for the stress inside the ellipse σ e . This effective stress provides the local loading environment for potential crack activation and growth. The local stress is resolved into tractions normal and tangent to crack faces. Since the frictional contact forces can only resist a finite shear traction across the crack faces, the excess traction is converted into a wedging force (Fw). This wedging force depends on the initial crack length (s), the angle (ϕ) between the crack face normal and the compression direction e^1, and the crack face coefficient of friction (μ). As discussed by Paliwal and Ramesh (2008), the wedging force is given by

⎡ e e e Fw = 2s ⎢⎣ μ σ11 cos2 (ϕ) + σ22 sin2 (ϕ) + σ12 sin (2ϕ) −

(

) ( (σ 1 2

e 11

⎤ e e − σ22 sin (2ϕ) − σ12 cos (2ϕ) ⎥⎦

The combination of the wedging force and the direct contribution from the wing crack tip given by

KI =

Fw π (l + 0.27s )

)

)

e s22

(C.1)

results in an effective stress intensity factor at

e + σ22 π (l + sin (ϕ) s ) .

(C.2)

When the local stress intensity factor (KI) exceeds the microscale fracture toughness (KIC), crack growth begins. Since we are particularly interested in capturing the competition between a low density of large flaws and a high density of small flaws, we require a crack growth law that captures the dynamic effects related to a moving crack tip. A suitable crack growth law was developed by Freund et al. (1972). In this model, the crack tip velocity (l )̇ asymptotically approaches the maximum crack growth velocity (vm), which is an experimentally determined fraction of the Rayleigh wave speed ( Cr ), as the stress intensity αc

factor increases. An additional empirical fitting parameter approaches its maximum value. The crack growth law is

l̇=

γ Cr ⎛ KI − KIC ⎞ c ⎟ . ⎜ αc ⎝ KI − 0.5KIC ⎠

γc is added to control the rate that the crack growth velocity

(C.3)

The dimensionless constants αc and γc can be determined experimentally by measuring the crack velocity as a function of the applied stress intensity factor. By adopting this crack growth law we are applying a macroscopic crack growth relationship in a micromechanics model. This application assumes that the microscale environment around the cracks can be treated as a continuum and that the Rayleigh wave speed and fracture toughness are well defined at this scale. These are reasonable assumptions; however, additional experimental measurements of the microscale fracture toughness of materials could help improve this mode. We assume that the microscale fracture toughness is similar to the macroscale fracture toughness and discard microstructural details such as the variation in fracture toughness between cleavage planes and grain boundaries. Under tensile loading conditions the crack faces will separate and the wedging force should be set to zero. Instead of explicitly checking the normal traction across the crack face, we set the wedging force to zero if Eq. (C.3) evaluates to a e negative number. When Fw is set to zero we recover KI = σ22 π (l + sin (ϕ) s ) , which is consistent with the stress intensity factor for a crack in a large plate under tension.

144

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

Appendix D. Corrected self consistent method solution (Hu et al., 2015, Appendix B) The compliance of the matrix material is given by equation (26); however, it is convenient to write out the non-zero components of the compliance tensor:

s1111 =

1 E0

(D.1)

s2222 =

8 (1 − ν02 ) 1 + ( D ( 4 − 2ν0 )) E0 3E0 (2 − ν0 )

(D.2)

s3333 =

1 E0

(D.3)

s1122 = s2233 =

2 (1 − ν02 ) −ν 0 − D E0 3E 0

(D.4)

s1212 =

8 (1 − ν02 ) 1 + ν0 D + 2E 0 3E0 (2 − ν0 )

(D.5)

s1133 =

−ν 0 E0

(D.6)

s1313 =

1 + ν0 2E 0

(D.7)

s2323 =

8 (1 − ν02 ) 1 + ν0 + D 2E 0 3E0 (2 − ν0 )

(D.8)

To use plane strain instead of plane stress, we need to compute the stiffness tensor (cijkl) by inverting sijkl:

⎡ c1111 c1122 c1133 ⎤ ⎡ s1111 s1122 s1133⎤−1 ⎥ ⎥ = ⎢s ⎢c c c s s ⎢ 1122 2222 2233 ⎥ ⎢ 1122 2222 2233⎥ ⎣ c1133 c2233 c3333⎦ ⎣ s1133 s2233 s3333⎦

(D.9)

From the stiffness tensor we can compute the reduced stiffness tensor. We then invert the reduced stiffness tensor to compute the reduced compliance tensor:

⎡ s 11 s 11⎤ ⎡ c ⎤−1 ⎢ 11 22 ⎥ = ⎢ 1111 c1122 ⎥ 11 22 ⎣ c c 1122 2222 ⎦ ⎢⎣ s22 s22 ⎥⎦

(D.10)

The inversion of the 2  2 matrix can be done analytically which gives the planar compliance components in terms of the components of the stiffness tensor: 11 s11 =

c2222 2 c1111c2222 − c1122

(D.11)

22 s22 =

c1111 2 c1111c2222 − c1122

(D.12)

11 s22 =

−c1122 2 c1111c2222 − c1122

(D.13)

12 s12 = s1212

(D.14)

For plane stress the reduction occurs for the compliance tensor not the stiffness tensor. This self consistent solution is an application of the work in “Theoretical Elasticity” by Green and Zerna. The following solution was developed by Hu et al. (2015) as a correction to Paliwal and Ramesh (2008), it is provided here for completeness. Starting from Section 9.1 of Green and Zerna we define an Airy stress function ϕ (z ) with the general form:

ϕ = Ω (z1) + Ω¯ (z¯1) + ω (z2 ) + ω¯ (z¯2 )

(D.15)

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

145

Here we have defined the complex variables zk and z¯k with k ¼1, 2 as (Eqs. (6.9.1) and (6.9.2)):

(

)

22 4 11 12 11 0 = s22 α − 2 s22 + 2s12 α 2 + s11

α12

α22

=

=

⎛ 11 ⎞ 12 ⎟+ 2 ⎜ s22 + 2s12 ⎝ ⎠

(

11 12 4 s22 + 2s12

(D.16) 2

)

11 22 s22 − 4s11

22 2s22

⎛ 11 ⎞ 12 ⎟− 2 ⎜ s22 + 2s12 ⎝ ⎠

(

11 12 4 s22 + 2s12

(D.17) 2

)

11 22 s22 − 4s11

22 2s22

(D.18)

αk =

αk2

(D.19)

γk =

αk − 1 αk + 1

(D.20)

zk = z + γk z¯

(D.21)

z¯k = z¯ + γ¯k z

(D.22) kl

Here the material parameters sij are the components of the planar compliance tensor (see Section 5.12 for the relation to the common stiffness matrix in Voigt notation). From Section 9.1 the displacements are given by (Eq. (9.1.4)):

D = u + iv = δ1Ω′(z1) + ρ1 Ω¯ ′(z¯1) + δ2 ω′(z2 ) + ρ2 ω¯ ′(z¯2 )

(D.23)

Here (Eqs. (6.9.4) and (6.9.5)): 11 22 2 βk = s22 − s22 αk

(D.24)

δ1 = (1 + γ1) β2 − (1 − γ1) β1

(D.25)

δ2 = (1 + γ2 ) β1 − (1 − γ2 ) β2

(D.26)

ρ¯1 = (1 + γ1) β2 + (1 − γ1) β1

(D.27)

ρ¯2 = (1 + γ2 ) β1 + (1 − γ2 ) β2

(D.28)

And the tractions are (9.1.8)

P = X + iY = 2i ( γ1Ω′(z1) + Ω¯ ′(z¯1) + γ2 ω′(z2 ) + ω¯ ′(z¯2 ) )

(D.29)

To account for the elliptical geometry we define the mapping:

z = z (ζ ) = cζ +

d ζ

(D.30)

c=

a+b 2

(D.31)

d=

a−b 2

(D.32)

zk = (c + γk d) ζk +

(d + γk c ) ζk

(D.33)

We now defined functions f (ζ1) and g (ζ2 ) such that

Ω′(z1) = f (ζ1)

(D.34)

ω′(z2 ) = g (ζ2 )

(D.35)

146

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

We try potentials of the form

f (ζ1) = H1ζ1 +

G1 ζ1

g (ζ2 ) = H2 ζ2 +

(D.36)

G2 ζ2

(D.37)

H1 = (B + iC )(c + γ1d)

(D.38)

H2 = (B′ + iC′)(c + γ2 d)

(D.39)

Evaluating the stress at large ζ, and requiring bounded displacements leads to Eqs. (9.1.16) and (9.1.17) in Green and Zerna, which can be used to solve for B , B′, C , and C′ (note N1 = σ 3 and N2 = σ1):

N1 + N2 = 4γ1 (B + iC ) + 4γ¯1 (B − iC ) + 4γ2 (B′ + iC′) + 4γ¯2 (B′ − iC′)

(D.40)

N1 − N2 = − 4γ12 (B + iC ) − 4 (B − iC ) − 4γ22 (B′ + iC′) − 4 (B′ − iC′)

(D.41)

Equation (9.1.18), which results from bounded displacements at infinity, is also required

0 = (1 − γ1γ¯1) ⎡⎣ (γ1 − γ2 )(1 − γ1γ¯2 )(B + iC ) − (γ¯1 − γ¯2 )(1 − γ¯1γ2 )(B − iC ) ⎤⎦ − (1 − γ2 γ¯2 ) ⎡⎣ (γ1 − γ2 )(1 − γ¯1γ2 )(B′ + iC′) − (γ¯1 − γ¯2 )(1 − γ1γ¯2 )(B′ − iC′) ⎤⎦

(D.42)

We can rearrange and group terms to set up the solution for B , B′, C , and C′:

N1 + N2 = 4 (γ1 + γ¯1) B + 4i (γ1 − γ¯1) C + 4 (γ2 + γ¯2 ) B′ + 4i (γ2 − γ¯2 ) C′

(D.43)

N1 − N2 = − 4 (γ12 − 1) B − 4i (γ12 − 1) C − 4 (γ22 − 1) B − 4i (γ22 − 1) C′

(D.44)

In the matrix at the boundary between the inclusion and the matrix we have ∥ ζ ∥ = 1 (∥ z ∥ = 1 = ∥ ζ ∥,

1 ζ

= ζ¯ and the

displacements are given by

D = δ1Ω′(z1) + ρ1 Ω¯ ′(z¯1) + δ2 ω′(z2 ) + ρ2 ω¯ ′(z¯2 )

(

(D.45)

)

(

D = δ1 ( H1ζ + G1ζ¯ ) + ρ1 H¯ 1ζ¯ + G¯ 1ζ + δ2 ( H2 ζ + G2 ζ¯ ) + ρ2 H¯ 2 ζ¯ + G¯ ζ

)

(D.46)

D = (δ1H1 + ρ1 G¯ 1 + δ2 H2 + ρ2 G¯ 2 ) ζ + (δ1G1 + ρ1 H¯ 1 + δ2 G2 + ρ2 H¯ 2 ) ζ¯

(D.47)

And the resulting force over an arc is given by

P = 2i ( γ1Ω′(z1) + Ω¯ ′(z¯1) + γ2 ω′(z2 ) + ω¯ ′(z¯2 ) )

( (

) (

)

(

(D.48)

) (

P = 2i γ1 H1ζ + G1ζ¯ + H¯ 1ζ¯ + G¯ 1ζ + γ2 H2 ζ + G2 ζ¯ + H¯ 2 ζ¯ + G¯ 2 ζ

(

P = 2i (γ1H1 + G¯ 1 + γ2 H2 + G¯ 2 ) ζ + (γ1G1 + H¯ 1 + γ2 G2 + H¯ 2 ) ζ¯

))

(D.49)

)

(D.50)

The inclusion is an isotropic material. If we assume that the complex potentials Ωi (z ) and ω¯ i (z¯ ) describe the material, the displacements and tractions are given by

μDi = kΩi (z ) − zΩ¯ i′(z¯ ) − ω¯ i (z¯ )

(D.51)

Pi = 2i (zΩ¯ i′(z¯ ) + Ωi (z ) + ω¯ i (z¯ ))

(D.52)

The material constant k is given by k = (3 − ν ) /(1 + ν ) for plane stress and k = 3 − 4ν for plane strain. These calculations are done for plane stress conditions because for plane strain the planar compliance tensor must be recalculated (compute the 3D stiffness tensor, perform the planar reduction, then invert that matrix to compute the planar compliance tensor). We try potentials of the form Ωi (z ) = A1 z , Ωi′(z ) = A1, and ω¯ i′(z¯ ) = A2 z¯ and apply the same conformational mapping (z = cζ + d/ζ ). Now we can evaluate the traction and displacement at the boundary of the ellipse (∥ z ∥ = 1 = ∥ ζ ∥, 1 = ζ¯ ): ζ

μDi =

( ( k − 1) A1c − A2 d) ζ + ( ( k − 1) A1d − A2 c ) ζ¯

Pi = 2i ( ( 2A1c + A2 d) ζ + ( 2A1d + A2 c ) ζ¯ )

(D.53) (D.54)

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

Displacements and tractions must be continuous across the ellipse boundary, and along the boundary dependent so we have 4 complex equations:

147

ζ and ζ¯ are in-

μ (δ1H1 + ρ1 G¯ 1 + δ2 H2 + ρ2 G¯ 2 ) =

( ( k − 1) A1c − A2 d)

(D.55)

μ (δ1G1 + ρ1 H¯ 1 + δ2 G2 + ρ2 H¯ 2 ) =

( ( k − 1) A1d − A2 c )

(D.56)

(γ1H1 + G¯ 1 + γ2 H2 + G¯ 2 ) = ( 2A1c + A2 d)

(D.57)

(γ1G1 + H¯ 1 + γ2 G2 + H¯ 2 ) = ( 2A1d + A2 c )

(D.58)

We can separate the known (H1 and H2) quantities from the unknowns (G1, G2, A1, and A2):

(k − 1) cA1 − dA2 − μρ¯1 G1 − μρ¯2 G2 = μδ¯1H¯ 1 + μδ¯2 H¯ 2

(D.59)

(k − 1) dA1 − cA2 − μδ1G1 − μδ2 G2 = μρ1 H¯ 1 + μρ2 H¯ 2

(D.60)

2cA1 + dA2 − G1 − G2 = γ¯1H¯ 1 + γ¯2 H¯ 2

(D.61)

2dA1 + cA2 − γ1G1 − γ2 G2 = H¯ 1 + H¯ 2

(D.62)

These 4 complex equations can be solved for the necessary complex constants, then the stresses in the ellipse are given by Eq. (8.1.18) in Green and Zerna (evaluated inside the ellipse on the real axis z = z¯ ):

σrr + σθθ = 4 ( Ωi′(z ) + Ω¯ i′(z¯ ) )

(D.63)

⎛ ⎞ z¯ σrr − σθθ + 2iσrθ = − 4 ⎜ z¯Ω¯ i″(z¯ ) + ω¯ i″(z¯ ) ⎟ ⎝ ⎠ z

(D.64)

Evaluating inside the ellipse these derivatives are

Ωi′ = A1

(D.65)

Ωi″ = 0

(D.66)

ωi′ = A2 z

(D.67)

ωi″ = A2

(D.68)

(

σrr + σθθ = 4 A1 + A¯1 )

)

(D.69)

σrr + σθθ = 8 Re (A1)

(D.70)

⎛ z¯ ⎞ σrr − σθθ + 2iσrθ = − 4 ⎜ z¯0 + A¯2 ⎟ ⎝ z ⎠

(D.71)

σrr − σθθ + 2iσrθ = − 4 Re (A2 ) + 4i Im (A2 )

(D.72)

With the coordinate system that we have set up σ11 = σrr , σ 22 = σθθ , and σ12 = σrθ :

σ11 = 4 Re (A1) − 2 Re (A2 )

(D.73)

σ 22 = 4 Re (A1) + 2 Re (A2 )

(D.74)

σ12 = − 2 Im (A2 ).

(D.75)

References Andrade, Jose E., Chen, Qiushi, Le, Phong H., Avila, Carlos F., Evans, T. Matthew, 2012. On the rheology of dilative granular media: bridging solid- and fluidlike behavior. J. Mech. Phys. Solids 60 (6), 1122–1136. Armero, F., Simo, J.C., 1993. A priori stability estimates and unconditionally stable product formula algorithms for nonlinear coupled thermoplasticity. Int. J. Plast. 9 (6), 749–782.

148

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

Ashby, M.F., Hallam, S.D., 1986. The failure of brittle solids containing small cracks under compressive stress states. Acta Metall. 34 (3), 497–510. Basista, M., Gross, D., 1998. The sliding crack model of brittle deformation: an internal variable approach. Int. J. Solids Struct. 35 (5–6), 487–509. Benz, W., Asphaug, E., 1994. Impact simulations with fracture. I. Method and tests. Icarus 107 (1), 98–116. Bolton, M.D., 1986. Strength and dilatancy of sands. Geotechnique 36 (1), 65–78. Brannon Rebecca, M., Gowen Travis, J., 2014. Aleatory quantile surfaces in damage mechanics. Journal of the European Ceramic Society 34.11, 2643–2653. Brannon, R.M., Wells, J.M., Strack, O.E., 2007. Validating theories for brittle damage. Metall. Mater. Trans. A: Phys. Metall. Mater. Sci. 38 A (12), 2861–2868. Brannon, R.M., Fossum, A.F., Strack, O.E., 2009. Kayenta: Theory and User's Guide. Technical Report SAND2009-2282. Sandia National Laboratories. Brannon, Rebecca M., 2007. Elements of Phenomenological Plasticity: Geometrical Insight, Computational Algorithms, and Topics in Shock Physics chapter Elements of Phenomenological Plasticity: Geometrical Insight, Computational Algorithms, and Topics in Shock Physics, ShockWave Science and Technology Reference Library, Springer, Berlin, Heidelberg, ISBN 978-3-540-22364-1, pp. 189–274. Budiansky, B., O'Connell, R.J., 1976. Elastic moduli of a cracked solid. Int. J. Solids Struct. 12 (2), 81–97. Carroll, M., Holt, A.C., 1972. Suggested modification of the p-α model for porous materials. J. Appl. Phys. 43 (2), 759–761. Chen, W., Luo, H., 2004. Dynamic compressive responses of intact and damaged ceramics from a single split Hopkinson pressure bar experiment. Exp. Mech. 44 (3), 295–299. Chocron, S., Anderson Jr., C.E., Dannemann, K.A., Nicholls, A.E., King, N.L., 2012. Intact and predamaged boron carbide strength under moderate confinement pressures. J. Am. Ceram. Soc. 95 (1), 350–357. Chocron, S., Anderson Jr., C.E., Dannemann, K.A., Nicholls, A.E., 2013. Pressure effects on the compressive response of confined intact and damaged sodalime glass. Exp. Mech. 53 (1), 77–89. Clayton, J.D., 2010. Deformation, fracture, and fragmentation in brittle geologic solids. Int. J. Fract. 163 (1–2), 151–172, http://dx.doi.org/10.1007/ s10704-009-9409-5. Curran, D.R., Seaman, L., Cooper, T., Shockey, D.A., 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), 53–83. Das, R., Cleary, P.W., 2010. Effect of rock shapes on brittle fracture using smoothed particle hydrodynamics. Theor. Appl. Fract. Mech. 53 (1), 47–60, http: //dx.doi.org/10.1016/j.tafmec.2009.12.004. Daux, C., Moës, N., Dolbow, J., Sukumar, N., Belytschko, T., 2000. Arbitrary branched and intersecting cracks with the extended finite element method. Int. J. Numer. Methods Eng. 48 (12), 1741–1760. Deng, H., Nemat-Nasser, S., 1992. Microcrack arrays in isotropic solids. Mech. Mater. 13 (1), 15–36. Deshpande, V.S., Evans, A.G., 2008. Inelastic deformation and energy dissipation in ceramics: a mechanism-based constitutive model. J. Mech. Phys. Solids 56 (10), 3077–3100. Deshpande, Vikram S., Nell Gamble, E.A., Compton, Brett G., McMeeking, Robert M., Evans, Anthony G., Zok, Frank W., 2011. A constitutive description of the inelastic response of ceramics. J. Am. Ceram. Soc. 94, s204–s214. Dienes, J.K., Zuo, Q.H., Kershner, J.D., 2006. Impact initiation of explosives and propellants via statistical crack mechanics. J. Mech. Phys. Solids 54 (6), 1237–1275. Drumheller, D.S., 1998. Introduction to Wave Propagation in Nonlinear Fluids and Solids. Cambridge University Press, Cambridge, UK. Falk, M.L., Needleman, A., Rice, J.R., 2001. A critical evaluation of cohesive zone models of dynamic fracture. J. Phys. IV France 11 (2001) Pr5-43-Pr5-50 http://dx.doi.org/10.1051/jp4:2001506. Foster, C.D., Regueiro, R.A., Fossum, A.F., Borja, R.I., 2005. Implicit numerical integration of a three-invariant, isotropic/kinematic hardening cap plasticity model for geomaterials. Comput. Methods Appl. Mech. Eng. 194 (50–52), 5109–5138. ISSN 0045-7825. http://dx.doi.org/10.1016/j.cma.2005.01.001. Freund, L.B., 1972. Crack propagation in an elastic solid subjected to general loading. II. Non-uniform rate of extension. J. Mech. Phys. Solids 20 (3), 141–152. Gailly, Benjamin A., Espinosa, Horacio D., 2002. Modelling of failure mode transition in ballistic penetration with a continuum model describing microcracking and flow of pulverized media. Int. J. Numer. Methods Eng. 54 (3), 365–398, http://dx.doi.org/10.1002/nme.427. ISSN 1097-0207. Graham-Brady, Lori, 2010. Statistical characterization of meso-scale uniaxial compressive strength in brittle materials with randomly occurring flaws. Int. J. Solids Struct. 47 (18–19), 2398–2413. Grechka, V., Kachanov, M., 2006. Effective elasticity of rocks with closely spaced and intersecting cracks. Geophysics 71 (3), D85–D91. Grechka, V., Kachanov, M., 2006. Effective elasticity of fractured rocks: a snapshot of the work in progress. Geophysics 71 (6). Guy, N., Seyedi, D.M., Hild, F., 2012. A probabilistic nonlocal model for crack initiation and propagation in heterogeneous brittle materials. Int. J. Numer. Methods Eng. 90 (8), 1053–1072. Hild, Franois, Denoual, Christophe, Forquin, Pascal, Brajer, Xavier, 2003. On the probabilistic deterministic transition involved in a fragmentation process of brittle materials. Comput. Struct. 81 (12), 1241–1253. ISSN 0045-7949. http://dx.doi.org/10.1016/S0045-7949(03)00039-7. Advanced Computational Models and Techniques in Dynamics. Holmquist, T.J., Johnson, G.R., 2006. Characterization and evaluation of boron carbide for plate-impact conditions. J. Appl. Phys. 100 (9). Holmquist, Timothy J., Johnson, Gordon R., 2011. A computational constitutive model for glass subjected to large strains, high strain rates and high pressures. J. Appl. Mech. 78 (5), 051003. Hu, Guangli, Liu, Junwei, Graham-Brady, Lori, Ramesh, K.T., 2015. A 3d mechanistic model for brittle materials containing evolving flaw distributions under dynamic multiaxial loading. J. Mech. Phys. Solids 78, 269–297, http://dx.doi.org/10.1016/j.jmps.2015.02.014. ISSN 0022-5096. Huang, Chengyi, Subhash, Ghatu, 2003. Influence of lateral confinement on dynamic damage evolution during uniaxial compressive response of brittle solids. J. Mech. Phys. Solids 51 (6), 1089–1105. Johnson, G.R., Holmquist, T.J., 1993. An improved computational constitutive model for brittle materials. High Press. Sci. Technol., 981–984. Johnson, G.R., Holmquist, T.J., Beissel, S.R., 2003. Response of aluminum nitride (including a phase change) to large strains, high strain rates, and high pressures. J. Appl. Phys. 94 (3), 1639–1646. Kimberley, J., Ramesh, K.T., Daphalapurkar, N.P., 2013. A scaling law for the dynamic strength of brittle solids. Acta Mater. 61 (9), 3509–3521. Knuth, Donald E., 1970. The art of computer programming. In: Semi-numerical Algorithms of Addison–Westley Series in Computer Science and Information Processing, vol. 2, second ed., Addison-Westley Pub. Co, Reading, Mass. Leavy, B., Strack, E., Brannon, R., Jensen, R., Houskamp, J., 2009. Simulation of experimental variability with spatially heterogeneous models. In: Society for Experimental Mechanics – SEM Annual Conference and Exposition on Experimental and Applied Mechanics, vol. 4, 2009, pp. 2290–2292. Leavy, R.B., Clayton, J.D., Strack, O.E., Brannon, R.M., Strassburger, E., 2013. Edge on impact simulations and experiments. Proced. Eng. 58, 445–452. Love, E., Sulsky, D.L., 2006. An unconditionally stable, energy–momentum consistent implementation of the material-point method. Comput. Methods Appl. Mech. Eng. 195, 3903–3925. Luo, Huiyang, Chen, Weinong W., Rajendran, A.M., 2006. Dynamic compressive response of damaged and interlocked SiC–N ceramics. J. Am. Ceram. Soc. 89 (1), 266–273. Macon, David James, Brannon, Rebecca Moss, Strack, Otto Erik, 2014. Plastic Cap Evolution Law Derived from Induced Transverse Isotropy in Dilatational Triaxial Compression. Technical Report SAND2014-1217. Sandia National Laboratories, Sandia National Laboratories (SNL), Albuquerque, NM, Livermore, CA, United States. Martin, B.E., Kabir, M.E., Chen, W., 2013. Undrained high-pressure and high strain-rate response of dry sand under triaxial loading. Int. J. Impact Eng. 54, 51–63. McCauley, J.W., Patel, P., Chen, M., Gilde, G., Strassburger, E., Paliwal, B., Ramesh, K.T., Dandekar, D.P., 2009. Alon: a brief history of its emergence and evolution. J. Eur. Ceram. Soc. 29 (2), 223–236. Meyer Jr., H.W., Brannon, R.M., 2012. A model for statistical variation of fracture properties in a continuum mechanics code. Int. J. Impact Eng. 42, 48–58. Molinari, J.F., Gazonas, G., Raghupathy, R., Rusinek, A., Zhou, F., 2007. The cohesive element approach to dynamic fragmentation: the question of energy convergence. Int. J. Numer. Methods Eng. 69 (3), 484–503.

A.L. Tonge, K.T. Ramesh / J. Mech. Phys. Solids 86 (2016) 117–149

149

Nemat-Nasser, S., Horii, H., 1982. Compression-induced nonplanar crack extension with application to splitting, exfoliation, and rockburst. J. Geophys. Res. 87 (B8), 6805–6821. Nemat-Nasser, S., Obata, M., 1988. Microcrack model of dilatancy in brittle materials. J. Appl. Mech., Trans. ASME 55 (1), 24–35. Nie, X., Chen, W.W., 2011. The influence of temperature and confinement pressure on the dynamic response of damaged borosilicate glass. In: Swab, Jeffery (Ed.), Advances in Ceramic Armor VII: Ceramic Engineering and Science Proceedings, vol. 32. Wiley, New York, pp. 3–10. Norton, Robert.L., 2000. Machine Design—An Integrated Approach, second ed. Prentice-Hall, Inc, Upper Saddle River, NJ. Paliwal, B., Ramesh, K.T., 2008. An interacting micro-crack damage model for failure of brittle materials under compression. J. Mech. Phys. Solids 56 (3), 896–923. Paliwal, B., Ramesh, K.T., McCauley, J.W., 2006. Direct observation of the dynamic compressive failure of a transparent polycrystalline ceramic (aion). J. Am. Ceram. Soc. 89 (7), 2128–2133. Paliwal, B., Ramesh, K.T., McCauley, J.W., Chen, M., 2008. Dynamic compressive failure of alon under controlled planar confinement. J. Am. Ceram. Soc. 91 (11), 3619–3629. Pandolfi, A., Ortiz, M., 2012. An eigenerosion approach to brittle fracture. Int. J. Numer. Methods Eng. 92 (8), 694–714. Parker, Steve, de St. Germain, J. Davison, Schmidt, John, Harman, Todd, Guilkey, James, et al., 2011. Uintah Website. Poncelet, E.F., 1946. Fracture and comminution of brittle solids. Trans. SME/AIME 169 (1946): 37-56. Radovitzky, R., Seagraves, A., Tupek, M., Noels, L., 2011. A scalable 3d fracture and fragmentation algorithm based on a hybrid, discontinuous Galerkin, cohesive element method. Comput. Methods Appl. Mech. Eng. 200 (1–4), 326–344. Riedel, W., Hiermaier, S., Thoma, K., 2010. Transient stress and failure analysis of impact experiments with ceramics. Mater. Sci. Eng. B 173 (1–3), 139–147. ISSN 0921-5107. http://dx.doi.org/10.1016/j.mseb.2009.10.038. 3rd International Conference on Science and Technology of Advanced Ceramics. Sadeghirad, A., Brannon, R.M., Burghardt, J., 2011. A convected particle domain interpolation technique to extend applicability of the material point method for problems involving massive deformations. Int. J. Numer. Methods Eng. 86 (12), 1435–1456. Sadeghirad, A., Brannon, R.M., Guilkey, J.E., 2013. Second-order convected particle domain interpolation (cpdi2) with enrichment for weak discontinuities at material interfaces. Int. J. Numer. Methods Eng. 95 (11), 928–952. Sevostianov, I., Kachanov, M., 2002. On elastic compliances of irregularly shaped cracks. Int. J. Fract. 114 (3), 245–257. Shockey, D.A., Marchand, A.H., Skaggs, S.R., Cort, G.E., Burkett, M.W., Parker, R., 1990. Failure phenomenology of confined ceramic targets and impacting rods. Int. J. Impact Eng. 9 (3), 263–275. Simo, J.C., 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), 221–245. Strassburger, E., 2004. Visualization of impact damage in ceramics using the edge-on impact technique. Int. J. Appl. Ceram. Technol. 1 (3), 235–242. Strassburger, Elmar, Patel, Parimal, McCauley, James, W., Templeton, Douglas, W., 2005. Visualization of wave propagation and impact damage in a polycrystalline transparent ceramic—ALON. In: 22nd International Symposium on Ballistics. Strassburger, E., Patel, P., McCauley, J.W., Templeton, D.W., 2006. High-speed photographic study of wave propagation and impact damage in fused silica and alon using the edge-on impact eoi method. AIP Conf. Proc. 845 (II), 892–895. Sulsky, D., Chen, Z., Schreyer, H.L., 1994. A particle method for history-dependent materials. Comput. Methods Appl. Mech. Eng. 118 (1–2), 179–196. Terzaghi, Karl, Peck, Ralph B., Mesri, Gholamreza, 1996. Soil Mechanics in Engineering Practice, third ed. John Wiley & Sons, Inc., New York. Tonge, Andrew L., 2014. A Unified Framework which Uses Multi-scale Microstructure Information for Modeling Dynamic Failure in Brittle Materials (Ph.D. thesis). The Johns Hopkins University. Wallstedt, P.C., Guilkey, J.E., 2008. An evaluation of explicit time integration schemes for use with the generalized interpolation material point method. J. Comput. Phys. 227 (22), 9628–9642. Warner, Charles T., Hartnett, Thomas M., Fisher, Donald., Sunne, Wayne., 2005. Characterization of alon optical ceramic. SPIE Proc., Dev. Spinel Alum. Oxynitride 5786, 95–111, http://dx.doi.org/10.1117/12.607596. Wilkins, MarkL., 1980. Use of artificial viscosity in multidimensional fluid dynamics calculations. J. Comput. Phys. 36, 281–303.