Aleatory Uncertainty and Scale Effects in Computational 

Damage Models for Failure and Fragmentation

by OE Strack, RB Leavy, and RM Brannon

ARL-RP-0497 September 2014

A reprint from Int. J. Numer. Meth. Engng., 3 June 2014.

Approved for public release; distribution is unlimited.
NOTICES

Disclaimers

The findings in this report are not to be construed as an official Department of the Army position unless 

so designated by other authorized documents.

Citation of manufacturer’s or trade names does not constitute an official endorsement or approval of the 

use thereof.

Destroy this report when it is no longer needed. Do not return it to the originator.
Army Research Laboratory

Aberdeen Proving Ground, MD 21005-5069

ARL-RP-0497 September 2014

Aleatory Uncertainty and Scale Effects in Computational 

Damage Models for Failure and Fragmentation

RB Leavy

Weapons and Materials Research Directorate, ARL

OE Strack

Sandia National Laboratories

RM Brannon

University of Utah

A reprint from Int. J. Numer. Meth. Engng., 3 June 2014.

Approved for public release; distribution is unlimited.
ii

REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering 

and maintaining the data needed, and completing and reviewing the collection information. Send comments regarding this burden estimate or any other aspect of this collection of information, 

including suggestions for reducing the burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson 

Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to 

comply with a collection of information if it does not display a currently valid OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD-MM-YYYY)

September 2014

2. REPORT TYPE

Reprint

3. DATES COVERED (From - To)

January 2011–January 2013

4. TITLE AND SUBTITLE

Aleatory Uncertainty and Scale Effects in Computational Damage Models for 

Failure and Fragmentation

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

OE Strack, RB Leavy, and RM Brannon

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

U.S. Army Research Laboratory

ATTN: RDRL-WMP-C

Aberdeen Proving Ground, MD 21005-5069

8. PERFORMING ORGANIZATION

 REPORT NUMBER

ARL-RP-0497

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR'S ACRONYM(S)

11. SPONSOR/MONITOR'S REPORT

 NUMBER(S)

12. DISTRIBUTION/AVAILABILITY STATEMENT

Approved for public release; distribution is unlimited.

13. SUPPLEMENTARY NOTES

A reprint from Int. J. Numer. Meth. Engng., 3 June 2014.

14. ABSTRACT

Stress concentrations near grain boundaries, precipitates, and similar micro-heterogeneities nucleate instabilities leading to macroscale 

fracture. As it is not practical to model each flaw explicitly, their ensemble effect is modeled statistically. Accounting for this aleatory 

uncertainty requires smaller specimens (e.g., small finite elements) to have generally higher and more variable strengths, which is 

necessary for the initial failure probability of a finite domain to be unaffected by its discretization into elements. Localization itself, 

which might be attributed to constitutive instability, requires realistic numerical perturbations to predict bifurcations such as radial 

cracking in axisymmetric problems. These perturbations, stemming from microscale heterogeneity, are incorporated in simulations by 

imposing statistical spatial variability in the parameters of an otherwise conventional (deterministic and scale-independent) damage 

model. This approach is attractive for its algorithmic simplicity and straightforward calibration from standard strength tests. In 

addition, it results in virtually no loss of efficiency or robustness relative to deterministic models and accommodates general threedimensional loading. Despite these advantages, some significant challenges remain and are discussed. However, it is demonstrated that 

including aleatory uncertainty with associated scale effects significantly improves predictiveness on large-scale computational 

domains, where it is impractical to resolve each crack or localization zone.

15. SUBJECT TERMS

constitutive equations, damage, fracture and fragmentation, fuzzy probabilistic methods, impact, rate dependence, scale/size effects, 

mesh dependence, verification, validation

16. SECURITY CLASSIFICATION OF: 

17. LIMITATION

OF ABSTRACT

UU

18. NUMBER

OF PAGES

34

19a. NAME OF RESPONSIBLE PERSON

R. B. Leavy

a. REPORT

Unclassified

b. ABSTRACT

Unclassified

c. THIS PAGE

Unclassified

19b. TELEPHONE NUMBER (Include area code)

410-278-7982

Standard Form 298 (Rev. 8/98)

Prescribed by ANSI Std. Z39.18
INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING

Int. J. Numer. Meth. Engng (2014)

Published online in Wiley Online Library (wileyonlinelibrary.com). DOI: 10.1002/nme.4699

Aleatory uncertainty and scale effects in computational damage

models for failure and fragmentation

O. E. Strack1, R. B. Leavy2 and R. M. Brannon3,*,†

1Sandia National Laboratories, Computational Shock and Multiphysics, Albuquerque, NM, USA 2Army Research Laboratory, Impact Physics, Aberdeen Proving Ground, MD, USA 3University of Utah, Mechanical Engineering, Salt Lake City, UT, USA

SUMMARY

Stress concentrations near grain boundaries, precipitates, and similar micro-heterogeneities nucleate

instabilities leading to macroscale fracture. As it is not practical to model each flaw explicitly, their ensemble effect is modeled statistically. Accounting for this aleatory uncertainty requires smaller specimens (e.g.,

small finite elements) to have generally higher and more variable strengths, which is necessary for the

initial failure probability of a finite domain to be unaffected by its discretization into elements. Localization itself, which might be attributed to constitutive instability, requires realistic numerical perturbations

to predict bifurcations such as radial cracking in axisymmetric problems. These perturbations, stemming

from microscale heterogeneity, are incorporated in simulations by imposing statistical spatial variability

in the parameters of an otherwise conventional (deterministic and scale-independent) damage model. This

approach is attractive for its algorithmic simplicity and straightforward calibration from standard strength

tests. In addition, it results in virtually no loss of efficiency or robustness relative to deterministic models

and accommodates general three-dimensional loading. Despite these advantages, some significant challenges remain and are discussed. However, it is demonstrated that including aleatory uncertainty with

associated scale effects significantly improves predictiveness on large-scale computational domains, where

it is impractical to resolve each crack or localization zone. Copyright © 2014 John Wiley & Sons, Ltd.

Received 31 January 2014; Revised 29 April 2014; Accepted 29 April 2014

KEY WORDS: constitutive equations; damage; fracture and fragmentation; fuzzy probabilistic methods;

impact; rate dependence; scale/size effects; mesh dependence; verification; validation

1. INTRODUCTION

This paper has been submitted as part of a special issue in honor of Professor Ted Belytschko,

whose contributions are pervasive in the finite-element community. As an example, the production

finite-element code that we use in this paper [1] is based on his seminal publications on hexahedral elements [2, 3]. Among his numerous contributions in the area of fracture mechanics, the

extended FEM [4] has enabled accurate modeling of discontinuities in finite-element meshes (often

introduced by cracks). As computational resources continue to evolve, methods such as extended

FEM might ultimately be used ubiquitously to model failure. Meanwhile, damage models, such as

the one summarized in this paper, remain useful for modeling ballistic impact and similar highenergy events. We have selected this topic for this special issue because of Professor Belytschko’s

continual support of our efforts to include strength statistics and scale effects in engineering damage

simulations.

This paper describes a practical method for extending a computational framework to accommodate aleatory uncertainty and associated size effects (e.g., spatial variability in strength and other

material properties) without prohibitively increasing computational overhead or model development

*Correspondence to: R. M. Brannon, University of Utah, Mechanical Engineering, Salt Lake City, UT, USA.

†E-mail: Rebecca.Brannon@utah.edu

Copyright © 2014 John Wiley & Sons, Ltd.
O. E. STRACK, R. B. LEAVY AND R. M. BRANNON

Figure 1. Tungsten-carbide (WC) sphere impacting a silicon carbide (SiC-N) ceramic at 375 m/s. All

simulations use the same constitutive model on the same three mesh resolutions, but those in the bottom row

are initialized using experimental data for aleatory uncertainty and associated size effects, giving relatively

mesh-insensitive radial cracking patterns similar to those in the (truly representative) experiment [7].

time relative to “conventional” (deterministic and scale-independent) damage models. As illustrated

in Figure 1, a conventional model’s predictions can be made more realistic and less mesh sensitive

when it is run using properties that include statistical variability and scale effects that are evident in

standard laboratory strength testing. Heterogeneity is applied to an existing conventional smeareddamage model at the host-code level so that few (if any) alterations to the underlying material

model’s source code are needed. Merits and drawbacks of the conventional damage model used

for this study [5, 6] receive only slight attention in this paper, whereas the physical foundations,

algorithms, and calibration procedures for uncertainty and size effect are covered in detail.

Aleatory uncertainty refers to natural spatial variation in flaw or inclusion morphology that causes

local regions (at a scale far smaller than a finite element‡) to be weaker or stronger than surrounding

material. Because it is impractical to precisely measure and model each flaw size and orientation,

certain material properties such as strength are distributed statistically in a manner consistent with

macroscale observations for a variety of specimen sizes. Slight spatial variability, which is relatively

inconsequential in stable deformations [9, 10], activates localized zones of damage (e.g., shear bands

or fractures) as the intensity of loading approaches an unstable state.

Conventional plasticity and damage models have long been known to admit bifurcation to a

localized state (when, e.g., the acoustic tensor has a zero eigenvalue; cf. [12, 13]). In the laboratory,

however, damage begins at stress levels well below this mathematical stability limit [14, 15]. The

situation is analogous to buckling: a column will buckle at a load well below the theoretical limit

load if its heterogeneities (holes and pits) are large enough to induce non-infinitesimal perturbations

in the stress field. If mesh morphology, material properties, and loading in a simulation are perfectly

homogeneous, then conventional deterministic damage models (including those with bifurcation

and decohesion criteria) are mathematically incapable of predicting heterogeneous fragmentation

of the type depicted in Figure 2. Without perturbations, an idealized simulation of homogeneous

loading (i.e., run at “infinite precision”) would reach a stability limit at the same moment for all

elements, producing unrealistic failure in every element in a computational domain, which is an

undesired material response that has thwarted numerous attempts to simulate dynamic failure of

armor ceramics [16, 17]. Because constitutive instabilities appear to be extremely sensitive to the

nature of perturbations, realistic perturbation fields must be used rather than relying on spurious

mesh-dependent and platform-dependent disturbances such as numerical round-off (cf. [18, 19]).

‡Simulations presented in this paper were run using a production-quality finite-element code [1], but the methods apply

equally well to other PDE solvers. The same techniques have, for example, been confirmed to give similar results when

finite elements are reinterpreted as particle domains in a material-point-method code [8].

Copyright © 2014 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2014)

DOI: 10.1002/nme
ALEATORY UNCERTAINTY AND SCALE EFFECTS IN COMPUTATIONAL DAMAGE MODELS

Figure 2. Heterogeneous dynamic response to nominally homogeneous (and initially quasistatic) unconfined

compression of a ceramic (SiC-N). The simulation uses aleatory uncertainty in strength (i.e., macroscale

statistical variability attributed to micro-heterogeneity) to capture the tendency of these experiments to fail

at only one end, typically leaving an intact fragment at the other end [11].

Figure 3. (a) Spurious mesh-texture-biased radial cracking in a conventional damage model compared with

(b) realistic randomly oriented cracking (with statistically reproducible spacing) with aleatory uncertainty.

For a SiC-N cylinder impacted by a WC sphere at 375 m/s, the first two simulations in sub-figure (b) use the

same mesh as sub-figure (a)—each with a different random seed, whereas the third uses a “paved” mesh.

A damage model is successful if it is able to reproduce (without mesh sensitivity) observed scaledependent and rate-dependent statistical bifurcation into damaged and intact zones with realistic

expended fragmentation energy.

In a systematic investigation of correlations between standard material properties (such as

hardness, grain size, porosity, etc.) and ballistic performance of a variety of SiC-N ceramics, Ray

et al. [20] found that “a wide variety of SiC-based materials can give good ballistic results, contradicting some of the theories about what is important to improve ceramics for armor.” Their

investigation revealed a possible correlation between the Weibull modulus and ballistic performance,

but they were quick to note that this result might have been coincidental because there was also

no correlation to the weak tail of the strength distribution. If these conclusions are independently

validated through further testing, then variability itself (not the value of the lowest strength) might

be the key factor in ballistic performance, suggesting that the salient numerical requirement is to

realistically seed an instability, perhaps of the nature discussed by Grinfeld [21].

For axisymmetric loading, such as indentation of a macroscopically homogeneous cylinder,

symmetry-breaking radial cracks originate from micromorphological instability [21]. Just as a

mathematical perturbation-free analysis of column buckling produces unrealistic simultaneous

failure at every point in the column, a perturbation-free mathematical analysis of axisymmetric

loading produces simultaneous failure of all points equidistant from the symmetry axis or excessive

cracking in a deterministic simulation such as Figure 3(a). Round-off noise represents a numerically infinitesimal perturbation many orders of magnitude smaller than the perturbations evident

in experimental data. As illustrated in Figure 3, measured statistics in strength should be used to

establish heterogeneous material property fields of statistically homogeneous media. When simulations rely exclusively on round-off error to stimulate instabilities, the results might break symmetry

without explicit control of perturbations, but such predictions have no physical basis, as they are typically caused by the following: (1) mesh-dependent or platform-dependent differences in the fields

or (2) order of operation errors in the code (where, for example the stress field at certain locations

is updated at the wrong time). In Figure 3(b), realistic random crack orientations exhibit reduced

Copyright © 2014 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2014)

DOI: 10.1002/nme
O. E. STRACK, R. B. LEAVY AND R. M. BRANNON

texture bias on the radial-trisection meshes. Timmel [22] notes that other mesh types might promote realistic failure patterns, but algorithms that minimize this mesh-texture bias are preferable.

The “paved” mesh in the rightmost image in Figure 3(b), for example, encourages slight curvature in the crack paths without any other pronounced differences in comparison with the radial

trisection results.

In addition to needing realistic statistical strength perturbations in simulations, an equally

important feature is the well-known concomitant dependence of strength on specimen size [23]. As

material failure tends to originate at the weakest point in a material, large samples will be weaker

(on average) because they are more likely to contain a critically large or critically oriented flaw.

Because the location of the weakest points in a material may not coincide with the points of highest stress, localized material failure does not generally originate at the highest stress concentrations,

nor does it necessarily precipitate cascading formation of macrocracks. With aleatory uncertainty

in simulations, non-cascading local damage (i.e., non-weakest-link behavior) is evident in the

“speckled” failure points in Figure 3(b) (which are very similar to those observed in computed

tomography data; cf. [24], which likewise do not all coalesce into macrocracks).

Using statistics in constitutive modeling dates back to the 19th century, where work was driven by

geomaterials because of their heterogeneity even at the macroscale. Statistics have been also used

extensively to predict concrete fracture [25] and munitions fragmentation or “hot-spot” formation

in energetic media [26–28]. Convergence of the statistical distribution of fragments generated by a

uniformly expanding ring was demonstrated by including Weibull variability of the yield strength

(without a size effect) [29]. Weibull statistics and size effects applied to cohesive zone models

improve their numerical convergence and predictive properties [30], without necessarily implying

that the overall strength is itself Weibull distributed.

In Section 2, we summarize a recent study [31] showing that strength distributions strongly

depend on the tensorial character of the loading mode, which leads to the mathematically rigorous

concept in Section 3 of “quantile” (isoprobability) surfaces that generalize the limit functions

of deterministic damage theories. Stress dependence of failure is characterized in Section 4 and

Section 5 from examination of experimental data, which suggest that both scale effects and aleatory

uncertainty decrease with increasing pressure. These observations motivate using a Weibull strength

distribution only for spherical (isotropic) tension, whereas non-Weibull distributions (as well as

decreased statistical variability with pressure) are naturally predicted for states having a nonzero

stress deviator. Data delocalization methods accounting for non-uniform stress fields in laboratory testing are discussed in Section 6, along with recent work on so-called data relocalization

needed in simulations to compensate for the inability of low-order shape functions to describe stress

variation on an element [32]. The constitutive damage model used in our case studies is briefly discussed in Section 7 (though our approach applies to any damage model). Finally, although failure

initiation receives the most attention in this paper, Section 8 emphasizes the need for including

scale-dependent damage rate effects during the failure progression phase of the simulation. Failure

initiation (e.g., the first growth of a microcrack leading to initiation of stiffness degradation) is

herein treated as a weakest link phenomenon, but there is no assumption that initiation of failure

will necessarily progress to a complete loss of stiffness or strength. Thus, overall structural failure

is not generally a weakest link event.

2. EXTENSION OF WEIBULL THEORY FOR GENERALIZED STRESS STATES

Weakest link theory (which we do not assume applies to all loading modes) is governed by Weibull

strength distributions, leading to Weibull scale effects in which large specimens (and hence large

finite elements) are weaker, on average, than small specimens (or small finite elements). For general

loading trajectories, our phenomenological extension of Weibull theory will be shown to exhibit

realistic non-Weibull statistical strength distributions and non-Weibull scale effects.

2.1. Failure theories and implications for specimen size dependence

If it were possible to create a sample containing exactly one small flaw, then onset of failure under a

given tensor stress state T (depicted in Figure 4 as a 3D Mohr’s circle) relates to the likelihood that

Copyright © 2014 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2014)

DOI: 10.1002/nme
ALEATORY UNCERTAINTY AND SCALE EFFECTS IN COMPUTATIONAL DAMAGE MODELS

Figure 4. Uncertainty in a single crack’s orientation (quantified by its unit normal, which is a point on the

unit sphere) allows the stress state (Mohr’s circle) to fall outside a failure line, depicted here as a straight

line corresponding to a criterion that the crack will fail if the resolved shear stress (minus friction) reaches a

critical value. Statistically variable crack sizes and failure modes perturb the failure line itself (cf. [31]).

the flaw is critically oriented or critically large. Figure 4, for example, shows that even a simplistic

Mohr–Coulomb failure criterion has uncertainty of failure if the crack orientation is random [31].

Similar observations have been made using more sophisticated single-crack failure theories [33],

and the same concepts have been applied to account for variable slip-plane orientations in plasticity

[34]. Some flaw orientations are safe from failure, whereas others (dark regions on the spheres in

Figure 4) are not. Brannon and Gowen [31] limit their scope to failure initiation in a statistically

isotropic medium (for which flaw orientation is initially uniformly random so that the initial failure

probability is the ratio of solid angle of the zone of critical orientations to the area of the sphere).

Using analytical methods, they derive failure statistics (probabilities that samples will fail under specific loading conditions) from uncertainties in flaw orientations, and they illustrate the corresponding

concept of strength quantile surfaces that generalize the notion of a limit failure surface used in

classical damage theories.

Dienes et al. [28] and Nemeth [35] accommodate intrinsic and deformation-induced anisotropy by

allowing flaw size distributions to be correlated to flaw orientation either from the outset or through

oriented crack growth. Macon et al. [36] further point out that induced anisotropy (leading to what

Nemeth [35] refers to as “equatorial” flaw distributions) may represent a way to model dilatation in

triaxial compression in a way that avoids anomalous predictions of hydrostatic strengthening under

dilatation seen in some isotropic hardening models (e.g., [37]). Rather than further expanding on the

physical motivations for these so-called unit-sphere representations of failure statistics, we focus on

actual implementations of these theories in engineering simulations.

2.2. Weibull theory (and its extension) for generalized stress states

Consider a stress state T D N, where is the magnitude of stress and N is a unit tensor called

the stress trajectory or loading mode. For example, uniaxial stress has N D DIAGŒ1; 0; 0, whereas

spherical tension has N D DIAGŒ1; 1; 1=p3 and biaxial tension has N D DIAGŒ1; 1; 0=p2. For

simplicity, our discussions might presume a straight path through stress space (where N is constant),

but the concept of strength quantile surfaces in six-dimensional (6D) stress space [31] admits more

tortuous tensorial paths through stress space in which the scalar (called strength) is merely the

value of the path parameter at failure. Just as deterministic plasticity yield surfaces generally evolve

with plastic flow, strength quantile surfaces are also expected to change in a history-dependent

manner so that probabilistic strength, , generally depends implicitly on evolving internal state variables such as crack density in addition to the loading mode. Let be the median strength observed

when testing a sample of volume V under a given loading mode N. To account for microscale heterogeneities (e.g., stress concentrations at flaws) and to accommodate associated scale effects, a

theory is required for predicting the probability, Pfail.T;V/, that a sample of a volume V (generally distinct from V ) will begin failing under an applied stress T. The Pfail.T;V/ function is

essentially a cumulative distribution function (CDF) in 6D stress space, so it is a Heaviside step

function in a deterministic theory (i.e., its gradient is singular at the limit surface). Otherwise, for

non-deterministic theories, the directed derivative along a stress trajectory gives the probability

Copyright © 2014 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2014)

DOI: 10.1002/nme
O. E. STRACK, R. B. LEAVY AND R. M. BRANNON

Figure 5. CDFs for normalized failure stress † for various normalized sample volumes,  ̨, and various

values of the Weibull modulus m.

density function (PDF) for trajectory-dependent strength as p.; V / @Pfail=@T



W N, which is the

key observation needed to construct the physically meaningful tensor-based function, Pfail.T;V/,

whose level sets are dubbed as “aleatory quantile surfaces”.

For a fixed stress trajectory N, define non-dimensionalized stress and size ratios:

† D 



and  ̨ D V

V : (1)

Weibull theory (which we soon argue does not apply to general loading modes but is nevertheless

useful for illustrating concepts in a familiar context) predicts that a sample of non-dimensionalized

volume  ̨ will begin failing under non-dimensionalized stress † with probability

Pfail D 1  Psafe D 1  2 ̨†m

; (2)

which is the initial CDF for strength, in which m is a constant called the Weibull modulus (or

shape parameter). The complementary CDF Psafe.T;V/ represents the probability that a sample of

volume V will not have begun failing under an applied stress T. Plots of Equation (2) are shown in

Figure 5 for various relative sample volumes and Weibull moduli, demonstrating that small values

of m correspond to large variability in strength. A purely deterministic model corresponds to a unit

step function in the limit m ! 1. Because the plots for small specimens (small  ̨) fall to the right

of those for larger specimens, small specimens have, on average, higher strengths. Equivalently,

because curves for small specimens fall below those of larger specimens, a small specimen is less

likely to fail at any given applied stress. Furthermore, the plots show that small specimens have

greater spread in their strengths.

2.3. Mode, median, and standard deviation of the Weibull distribution

For the Weibull distribution introduced in Equation (2), the median †median is the 50% quantile,

representing the stress level at which half of the specimens in a large population will have initiated

failure. Hence, †median is the value of † at which Pfail D 1=2, and referring to Equation (2), it

follows that the median strength depends on specimen size according to

†median D  ̨1=m: (3)

The probability density function (PDF), p.†;  ̨/ is the derivative of the CDF with respect to †:

p.†;  ̨/ @Pfail

@† D .ln 2/ 

2 ̨†m 

 ̨m†m1 (4)

The so-called Weibull scale parameter †scale, the mode †mode (at the peak of the PDF), and the

mean, †mean. ̨/ D R 1

0 † p.†;  ̨/ d† are found from

†scale

†median D .ln 2/ 1

m ; †mode

†median D

m  1

m ln 2

 1

m

; and

†mean

†median D  

1 C 1

m

	

Œln 21=m (5)

where  is the gamma function. Variability (spread) in the data may be quantified by the standard

deviation, †s.d., illustrated in Figure 6(a) and defined by

Copyright © 2014 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2014)

DOI: 10.1002/nme
ALEATORY UNCERTAINTY AND SCALE EFFECTS IN COMPUTATIONAL DAMAGE MODELS

(a) Standard deviation (b) Initialization

Figure 6. The graphs in (a) imply that, in addition to higher strength, greater variability must apply at smaller

elements within a computational domain, such as those near the centerline in (b), giving an unavoidable

impression of heterogeneity in what is nevertheless a statistically homogeneous strength distribution so that a

physical domain’s failure initiation probability under homogeneous loading is unaffected by the meshing of

that domain.

†s.d.

†median D

s






2 C m

m

 






1 C

1

m

2

: 



p6m for large m (6)

The plot of the standard deviation †s.d. in Figure 6(a) illustrates our earlier assertion that one

should not rely on round-off error to initiate mathematical instabilities. Even for the largest Weibull

modulus in Figure 6(a), the standard deviation in strength can be on the order of 15%. In fact, for the

left-hand side of Equation (6) to equal a typical round-off error of 1016, the Weibull modulus would

need to be more than 10 orders of magnitude larger than reported for any polycrystalline material

to our knowledge. Hence, perturbations that precipitate instabilities or bifurcations associated with

localization in a numerical code must be realistically seeded in order to obtain realistic results. A

deterministic model unrealistically has no variability (zero standard deviation) in strength.

Hild et al. [38] identify a characteristic sample (element) size marking a statistical-deterministic

transition such that larger elements are treated deterministically while variability is applied for

smaller elements. That approach admits nonphysical jump discontinuities in the standard deviation,

whereas Figure 6(a) shows that seeding all elements with a Weibull distribution results in a smooth

transition toward a deterministic response (zero standard deviation) as the element size increases.

Both higher strength and greater variability are evident in the time-zero initialization of a statistically homogeneous material in the computational domain of Figure 6(b), which was generated using

concepts developed in the next section to accommodate smooth, trajectory-dependent, and generally

non-Weibull seeding of strength for 3D stress states.

Failure probability functions are valuable in engineering if one seeks a probability of failure

initiation for the purpose of, say, component reliability and safety assessment [39, 40]. Theories

are also available for expected values of fragment sizes in impact and shock loading after failure

initiation [41, 42]. However, if one is interested in detailed analysis of the dynamic effects of the

failure process, then realizations of strength, consistent with the Psafe.T;V/ function, are required

to generate realistic crack spacings and fragment sizes [43]. Simulating realistic damage and fragmentation events is especially important for assessing degradation of strength that could lead to

particular vulnerabilities to subsequent impact events. The next section of this paper recognizes

and accounts for the existence of a Psafe.T;V/ function without requiring that it ever be explicitly

constructed in simulations. Accordingly, we will now focus on the goal of converting failure

probabilities into strength realizations on a computational domain.

Copyright © 2014 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2014)

DOI: 10.1002/nme
O. E. STRACK, R. B. LEAVY AND R. M. BRANNON

3. “FUZZY FAILURE” THEORY AND STRENGTH REALIZATIONS

Classical deterministic plasticity models presume existence of a scalar-valued yield function,

f .T;V/, such that a stress state T applied to a specimen of volume V is elastic if f .T;V/<0. The

yield surface is the set of T for which f .T;V/ D 0. Scale-insensitive models omit dependence on

V . The yield function also implicitly depends on internal variables, such as dislocation density and

crack texture parameters, which can evolve to allow the yield surface to expand (harden) or contract

(soften) in response to plastic loading. The yield surface bounds the domain of elastically attainable

states; engineering damage theories also use a so-called limit function F.T;V/ that defines the

boundary of stress states attainable by any quasistatic means [44]. A quasi-brittle model employs

distinct yield and limit functions .f ¤ F /, whereas a purely brittle model uses coincident functions .f D F /. Damage corresponds to collapse of the limit surface and concomitant degradation

of elastic stiffness.

These conventional sharp-threshold plasticity models have required considerable time and

resources to develop. Many (such as [45] and [46]) are also equipped with large material property

databases, compiled over many years. By presuming that existing discrete threshold functions are

actually fitted only to median experimental data, so-called “fuzzy failure theory” incorporates spatial variability and scale effects in these models. Moreover, because of the size effect, these fits apply

only to the specimen size used in the calibration experiments, where specimen size is defined using

“data-delocalization” methods as explained in Section 6.

3.1. Aleatory uncertainty theory

Conceptually, spatial variability is incorporated into these conventional plasticity models by

replacing a limit function F.T;V/ by a scale-dependent probability field Psafe.T;V/ that ranges

continuously from a value of 1 at stress states that are safe from failure to a value of 0 at stress

states that are certain to induce failure. The conventional limit function could, for example, be identified with the 50% quantile isosurface: F.T;V/ D Psafe.T;V/  1=2. Of course, similar revision

of the yield function applies to quasi-brittle models. Strength is statistically variable because of

aleatory uncertainty in the micromorphology (such as crack densities and textures). Recalling that

F.T;V/ depends implicitly on internal variables that directly or indirectly account for micromorphology, it follows that a natural way to indirectly incorporate the Psafe.T;V/ in a simulation is

to retain an existing code infrastructure for a scale-insensitive F.T/, but to initialize each finite

element with different F.T/ functions by statistically perturbing the parameters defining F.T/ in

a way that recovers observed strength probabilities captured in Psafe.T;V/, where V is taken as

the finite-element volume. As this initialization phase of the simulation is the only point at which

Psafe.T;V/ is needed, there is no algorithmic burden to otherwise incorporate this function into

run-time calculations, making the implementation of aleatory uncertainty quite manageable.

The Psafe.T;V/ function allows the effects of micromorphology to be treated in ensemble through

empirical strength CDFs. With this engineering approach, the median (which has been presumably

already fit to laboratory data) is retained but scaled appropriately with specimen (finite-element)

size. Statistically perturbed realizations of strength surfaces in stress space are centered about

the scale-adjusted median (i.e., 50% quantile) limit surface according to the statistics evident in

laboratory data (and consistent with limit surface admissibility requirements such as convexity).

“Fuzziness” in the limit surfaces of conventional plasticity models has been formalized by

Brannon and Gowen [31], who introduce the concept of quantile (isoprobability) surfaces in stress

space, replacing the single failure surface used in deterministic models by the family of surfaces

in Figure 7(a). The red in Figure 7(a) depicts the median (50% quantile) isoprobability surface

viewed down the hydrostat, with the 5% and 95% ‘quantile surfaces depicted as blue. The images

in Figure 7(b) are similar except that quantile surface colors are replaced with statistical dot densities centered about the median quantile surface. The “weak tail” of a strength distribution refers to

the smaller quantile surfaces, which (being closer to the origin in stress space) represent instances

of low-stress failure. For the idealized micromorphology used to generate Figure 7(a), the median

50% quantile surface (colored red) is a rounded triangle (modeled using the Willam–Warnke thirdCopyright © 2014 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2014)

DOI: 10.1002/nme
ALEATORY UNCERTAINTY AND SCALE EFFECTS IN COMPUTATIONAL DAMAGE MODELS

(a) Quantile surfaces

Mohr-Coulomb Willam-Warnke

(b) Fuzzy failure surfaces (octahedral view)

Figure 7. Accounting for variability in flaw size affects the location of the failure threshold itself, producing

quantile surfaces in stress space depicted here looking down the hydrostat and visualized according to initial

failure probability density by (a) colors showing theoretical shapes of each quantile surface [31] or (b) dot

densities that were generated using self-similar quantile shapes (more like the current implementation).

invariant function common in geomechanics; cf. [6]), whereas the weak 5% quantile is the blue

distorted hexagon that is easily modeled with Mohr–Coulomb third-invariant dependence. The

quantile functions in Figure 7(b) did not include this recent finding of Lode-angle sensitivity.

In our fuzzy failure theory, natural subscale heterogeneities are accommodated in a macroscale

simulation by using different realizations of strength in each finite element, such that the median

strength of a finite volume is preserved, with a key goal being that failure-initiation statistics ought

to be unaffected by mesh refinement. In numerical implementations, each finite element is regarded

as a finite-sized specimen requiring assignment of strength consistent with scale effects observed in

the laboratory. For illustration purposes, suppose that the stress at failure for a specimen size V

is known to be Weibull distributed with median for a given stress trajectory N (keep in mind that

non-Weibull strengths are still possible for other trajectories). To assign strength properties to a finite

element, one must generate a uniformly distributed random number R on the interval 0 6 R<1.

The realization of failure stress can then be found by replacing Psafe with R in Equation (2) and

solving for †. The result is

† D


 ln R

 ̨ ln.1=2/1=m

i.e., D 

"

V ln R

V ln.1=2/#1=m

: (7)

For a different sample volume V ¤ V



, this equation predicts a different strength. In our case, the

strength assigned to a finite element depends on the volume V of the element relative to the sample

volume used in laboratory testing. Equal-sized regions always have the same probability of failure

initiation under homogeneous loading regardless of whether or not those regions are subdivided

into many or few elements. As a special case, Equation (7) ensures that a collection of small elements having total volume V will have the same failure initiation statistics (under uniform loading)

as observed in the original laboratory sample of volume V (non-uniform loading is addressed in

Section 6). Equation (7) therefore contains an implicit length scale that ensures mesh-insensitive

predictions for the onset of failure under homogeneous loading. This is verified numerically for

a spherical stress trajectory, N D DIAGŒ1; 1; 1=p3, in Section 5 (Figure 14). As discussed in

Section 7, mesh insensitivity for the subsequent progression of failure is an active area of research in

which nonlocal effects and energy release considerations are important for mathematical, numerical,

and physical reasons.

Plots of the material strength associated with each realization parameter R in Equation (7), and

strength distribution curves associated with multiple realizations, are shown in Figure 8.

3.2. Algorithm for implementing aleatory uncertainty

A computational procedure to support aleatory uncertainty (fuzzy thresholds) in material properties

can be designed such that an existing code infrastructure is impacted minimally and so that existing

Copyright © 2014 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2014)

DOI: 10.1002/nme
O. E. STRACK, R. B. LEAVY AND R. M. BRANNON

Figure 8. Strength realization functions for various Weibull moduli, given by Equation (7), and probability

density functions for failure stress, given by Equation (4).

conventional deterministic material models require no internal revisions of their subroutines. Our

approach is outlined below

(1) Modify the user input parser so that material parameters that were formerly deterministic and

specified (for example) by a syntax

PROP D value

are treated statistically whenever the user employs a syntax (for example) of the form

PROP D Weibull(median, Vref, Wmod)

where median is the median measured in the laboratory using a sample of volume Vref and

Wmod is the Weibull modulus. Because median usually equals value from the deterministic

model, this approach to aleatory uncertainty requires two new user inputs (Vref and Wmod),

both of which will be available if the original strength measurements were accompanied by

specimen size and statistically significant standard deviations (as should be the case if testing

followed ASTM and similar standards).

(2) For each statistical user input, allocate storage for one field variable, giving it a name such as

PROPSR (where “PROP” is the property name and “SR” stands for statistical realization).

(3) Within an initialization loop over each finite element, compute a uniformly distributed random

number R on Œ0; 1/, and use the finite-element volume V in Equation (7) (or any other distribution of interest) to assign a realization value “” for the statistical property. Save the result

in the PROPSR field variable array.

(4) During each time step, before calling the constitutive model, construct a copy of the material

property array in which PROP is replaced by PROPSR, and send this revised property array to

the material model.

Virtually all of the CPU overhead resides in step 3, which is a once-only initialization setup task

typically requiring only a tiny fraction of overall CPU for the calculation. Because step 4 requires

nothing more than a small amount of additional memory storage and an array copy, this approach

to aleatory uncertainty entails almost no change in overall CPU costs relative to the deterministic

model. Because the constitutive model is treated as a “black box”, the necessary code architecture

may be written in a way that supports all constitutive models in the code and not just the damage

Copyright © 2014 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2014)

DOI: 10.1002/nme
ALEATORY UNCERTAINTY AND SCALE EFFECTS IN COMPUTATIONAL DAMAGE MODELS

model of interest. Users must be cautioned, however, that carelessly or arbitrarily perturbing parameters might violate the model’s admissibility conditions (such as the convexity requirements). As

discussed in the next section, achieving physically reasonable distributions for material parameters

might be facilitated by introducing a change of variables.

4. DIRECTIONALLY-WEIBULL THREE-INVARIANT THEORY

The aleatory quantile function, Psafe.T;V/, is the probability that a sample of size V is safe from

failure at an applied stress T. This function is used indirectly in implementations of aleatory uncertainty. Specifically, the material parameters implicitly contained in the limit function F.T;V/ are

statistically perturbed, thus giving each finite element its own perturbed limit surface distinct from

other elements. If F.T;V/ depends on only a single material parameter (as for von Mises or Tresca

theory), then model calibration requires measured strength statistics along only a single loading trajectory N (such as simple shear). In contrast to calibration, validation requires additional data along

other trajectories (such as uniaxial compression, biaxial tension, etc.). Such efforts might imply the

need for a three-parameter limit function, for which calibration requires strength data along three

non-crossing trajectories, as later discussed in Section 5.2.

Measured strengths are interpreted in the context of aleatory quantile (i.e., constant failure initiation probability) surfaces embedded in full 6D tensorial stress space. These multi-dimensional

surfaces are typically visualized as 2D cross-sections of 6D stress space. The “octahedral diagram”

in Figure 7, for example, is a visualization as seen looking down the hydrostat in principal stress

space, where the median (50%) quantile surface is surrounded by other quantile surfaces to define a

“cloud” that replaces the discrete limit surface in deterministic models. The “meridional diagram”

in Figure 9 depicts a different 2D slice of 6D stress space to show that equivalent shear strength

 pJ2 increases with pressure p I1=3, where I1 and J2 are the first and second invariants

of the stress tensor T, defined I1 D trT and J2 D 1

2 tr.devT/2. Factors of p2 and p3 in the axes

labels make the plot isomorphic to stress space (ensuring that the 6D unit tensor N maps to a 2D unit

vector that also forms the correct angle with the hydrostat). In upcoming plots, isomorphic scaling

factors are omitted whenever geometric accuracy is not needed. Quantile surfaces in Figure 9 are

depicted as a gray-scale “cloud”. The red line is the median (50% quantile), taken from the limit

surface in a deterministic damage model. The 5% and 95% quantiles are lighter-gray contours.

The loading increment direction is typically fixed in model calibration experiments, meaning that

the stress tensor is of the form T D M C C, where is a time varying load intensity parameter,

M is a unit tensor in the loading increment direction, and C is a constant tensor that is usually an

elastic prestress reference state. If, for example, the experiment involves a hydrostatic compression

to a pressure p followed by axisymmetric compression with lateral stress held constant at p, then

Figure 9. Loading trajectory through a fuzzy failure “cloud”.

Copyright © 2014 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2014)

DOI: 10.1002/nme
O. E. STRACK, R. B. LEAVY AND R. M. BRANNON

C D DIAGŒp; p; p and M D DIAGŒ1; 0; 0. If C D 0; M is the same as N defined in Section 2.2

and sketched in Figure 9. For these experiments, strength distributions are actually distributions

for the scalar loading parameter , and the quantiles for are then simply mapped to full 6D

stress space to generate points for fitting the corresponding quantile surfaces. If is the median,

for example, then T D M C C must be a point on the median quantile surface. Points on other

quantile surfaces in 6D stress space are similarly identified. The process of finding stress tensors on

quantile surfaces must be repeated for at least as many distinct loading paths as there are distinct

material parameters in the limit function F.T;V/, thus providing a means of setting all material

parameters corresponding to any quantile realization of strength. The strength in Equation (7) does

not necessarily represent a single component of stress or even any particular stress invariant. It is

merely a scalar parameter in the range from 0 to 1, whose distribution p.; V / indirectly quantifies

variation of the tensor-based CDF Pfail.T;V/ along a fixed loading trajectory M. Specifically, the

generally meaningless scalar is mapped to the physically significant stress tensor, T D M C C.

The scalar-based PDF p.; V / then provides information about the tensor-based CDF Pfail.T;V/

through the directed derivative relationship, p.; V / D .@Pfail=@T/ W M.

It should be clear from examination of Figure 9 that the absolute spread in failure strength data

(quantified by variance of ) will be affected by the angle of attack in stress space. More broadly,

each different loading direction M (i.e., each test, such as uniaxial tension, compression, shear, etc.)

is associated with a different distribution of strength. They are connected into a unified theory by fitting stress tensors (not ) on each quantile surface to the limit function F.T;V/ of a “conventional”

deterministic damage model.

For our own baseline “conventional” deterministic damage constitutive theory, we created a new

model, called Kayenta [6], which revised an existing geomaterial model [47] to incorporate loss

of strength in a manner similar to the Johnson–Holmquist class of ceramic models [45, 48–50].

Kayenta is a generalized isotropic damage model (as opposed to an anisotropic model, e.g., [28],

or a combined plasticity/linear-fracture-mechanics model, e.g., [51]). Similar to the strain-based

Johnson–Cook fracture [46] and Johnson–Holmquist models, Kayenta allows equivalent shear

strength  pJ2 in axisymmetric compression to depend on pressure p I1=3 according to

 D a1  a3ea2p (8)

where 1a1; a2; a3o are empirical fitting parameters (an additional term, Ca4p, is also available for

linear Drucker-Prager (and similar) functions but is omitted here for simplicity).

The Johnson–Holmquist and Johnson–Cook models fall into the category of nonlinear Drucker–

Prager generalized plasticity (i.e., their octahedral profile is always a circle). Kayenta incorporates

the third invariant (the Lode angle ) by multiplying  by a function of  and p that equals unity in

axisymmetric compression but is larger than unity at other Lode angles (giving non-circular octahederal profiles such as those in Figure 7). In Kayenta, the Lode angle multiplier function also

includes pressure dependence that allows the octahedral profile to vary with pressure, smoothly (and

automatically) transitioning from a maximum principal stress (or strain) criterion at low pressures

to a von Mises criterion at high pressures. A cap for porous media is also available as an option but

is not used in the sample calculations herein.

The limit surface defined by Equation (8) could be perturbed by changing one or more of the

1a1; a2; a3o parameters. However, these parameters ought not to be varied independently. The user

must ensure that PDFs used to perturb material properties will not violate any fundamental physical

constraints. Moreover, the user must recognize that achieving an observed Weibull (or other) distribution in strength might require complicated non-Weibull joint distribution functions for the model

parameters. For example, the constraint a1 > a3 must be satisfied to ensure positive shear strength

at zero pressure. Other considerations, such as convexity of the limit surface, place additional coupled constraints on the 1a1; a2; a3o parameters. None of these strength-defining parameters are

candidates for independent Weibull (or similar) perturbation. To address this inconvenience, we

re-parameterized Equation (8) into the form

 D Y

h

1  es.pCh/=Y i

; (9)

Copyright © 2014 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2014)

DOI: 10.1002/nme
ALEATORY UNCERTAINTY AND SCALE EFFECTS IN COMPUTATIONAL DAMAGE MODELS

Figure 10. Meridional profiles (including probability clouds depicting variability) corresponding to (a)

shear dominant variability obtained by perturbing Y , which failed to match observed damage patterns in

ceramic fragmentation simulations, and (b) pressure-dominant variability obtained by perturbing h, which

is consistent with both spall and Hugoniot elastic limit data and gives more realistic damage patterns

in simulations.

where Y is the shear strength limit, h is the spherical tensile strength, and s is the initial slope. The

meanings of these alternative parameters are labeled in Figure 10, which also shows probability

density clouds resulting from statistically perturbing Y or h individually.

The Weibull modulus for SiC was reported to be 17 by the manufacturer Coorstek-Vista [52],

but spall data [53] for SiC-N seem to have a Weibull modulus of only 4 (i.e., much more variable).

Brazilian tensile strengths (presented later in Figure 17) appear to have a Weibull modulus of about

8. Hugoniot data for B4C [54] has a Weibull modulus of 60, with similar numbers for SiC. Moreover,

it is not even clear that some of these data sets are very well fitted to a Weibull distribution. To match

the general character of observed strength data (i.e., lessening of variability with pressure and nonWeibull strength distributions for most loading trajectories), it is proposed that the limit surface be

perturbed as sketched in Figure 10(b).

Unlike the original 1a1; a2; a3o fitting parameters, there is no coupling of the constraints on

Kayenta’s alternative limit-surface parameters 1Y; h; so, which only must be positive. Hence, any of

them may be Weibull perturbed independently according to Equation (7) (where in that equation

is replaced with Y;h, or s). Consistent with a similar statistical uncertainty assumption used to analyze the nearly spherical tension states that exist during spall [55], and also consistent with the

tendency (evident in Figure 13) of strength measurements to be more variable at low pressure, we

take the h parameter to be Weibull distributed, whereas the other limit surface parameters (Y and

s) are deterministic. This choice leads to non-Weibull statistics and non-Weibull scale effects for

all non-spherical loading paths. This approach, which truly accounts for the 6D tensorial nature of

stress, differs significantly from related efforts [56] that apply a 1D Weibull distribution to a single

stress invariant such as equivalent shear stress, maximum principal stress, and so on. We conjecture

that variability decreases with pressure because the additional confinement provides friction that

can lock crack faces and hence suppress brittle failure modes in shear. Of course, high-pressure data

often comes from high-rate experiments, but we have tentatively ruled out the possibility that variability decreases with increasing loading rate; this choice is motivated by the fact that quasistatic

indentation experiments also seem to exhibit pressure-based variability because they tend not to

bifurcate into symmetry-breaking radial cracks in regions of high confinement such as immediately

beneath an indenter [57]. Further insight into the influence of pressure and loading rate is provided

in a recent investigation of rate effects in brittle media [58], which applies micromorphological

dimensional analysis to obtain a universal curve which closely fits rate-dependent and statistically

variable strength for several different cracked and porous brittle materials (including various rocks,

ceramics, and concrete).

In Figure 11(a), the value of the spherical tension limit h is exaggerated to illustrate the character

of perturbation effects. The actual value of h for brittle ceramics is typically so small in comparison

with the high pressure strength Y that the actual median limit curve might appear to pass through

Copyright © 2014 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2014)

DOI: 10.1002/nme
O. E. STRACK, R. B. LEAVY AND R. M. BRANNON

(a) Sample size effect (b) Tension plot

Figure 11. The strength size effect in “directed-Weibull” theory varies with stress space trajectory. Here,

strength can increase without bound for progressively smaller samples in spherical tension (horizontal

green loading line in (a) and straight green line in (b)), but it naturally saturates to a classical plasticity

limit in uniaxial tension (black angled loading line in (a) and black curved line in (b)), with clear

deviation from Equation (3) for non-spherical loading. This general framework naturally accommodates

any other experimentally observed distributions of quantile surfaces, including non-monotonic scaling

of strength.

the origin for lab scale .cm3/ samples, as illustrated in Figure 13(a) for a large .> 1 cm3/ sample

of SiC-N. Using Figure 7 to generate strength envelopes means that realizations of the 1Y; h; so

parameters are scale dependent. Smaller samples have, on average, higher strength. To illustrate,

Figure 13(b) shows the median limit surface for a grain-scale .m3/ sample.

When only the h parameter is Weibull-perturbed, the median spherical tensile strength approaches

infinity as the element size approaches zero. However (as seen in Figure 11), the median uniaxial

tension and shear strengths have finite limits, and the strength realizations in these non-spherical

directions deviate realistically from a Weibull distribution. This “directionally-Weibull” approach

gives a scaling effect similar to that reported by Bažant [23] on the basis of experimental data

for concrete.

Bažant has amassed direct experimental data for concrete supporting deviation from the Weibull

size dependence similar to what is depicted in Figure 11(b). No similar direct evidence of the

deviation from the Weibull size effect has (to our knowledge) been collected for sufficiently small

specimens of structural ceramics (whose “aggregate” size is orders of magnitude smaller than that

of concrete). We arrived at this effect indirectly by recognizing that macroscale tensile strengths are

far more variable than compressive strengths (Figure 13), suggesting that limit surface realizations

require perturbations of the type in Figure 10(b).

The nature of the limit surface “variability cloud” and its median master surface strongly influences the appearance of cone or radial cracks in dynamic indentation. In indentation experiments

(discussed in Section 8), cone cracks tend to emerge around and under the area of contact where

the pressure is extremely high. They are also symmetric [59], indicating a low degree of spatial

variability in the high-pressure realm. Radial cracks, on the other hand, tend to appear in relatively low-pressure regions, and they significantly disrupt symmetry, indicating a high degree of

spatial variability in strength akin to Figure 10(b). As seen in Figure 12, similar observations have

been made for brittle coatings [57, 60], glass [61], and boron carbide plates [62]. The order of the

appearance of radial cracks and cone cracks can be altered by changes in materials and indenter

geometries.

Mapping out aleatory quantile surfaces is clearly labor intensive. Given that it is often impossible

or impractical to obtain statistically significant samples of laboratory-measured strengths, judicious

use of virtual data (e.g., [64, 65]) might help to complete the aleatory model parameterization and/or

to identify appropriate functional forms of the limit function F.T;V/ about which perturbations

must be applied. Virtual data (if validated against available data and verified against classical idealized exact solutions) provides insight into the microphysical sources of aleatory uncertainty, which

likewise aids in the formulation of path-dependent macroscale models for evolution of internal state

variables that control the loss of stiffness and collapse of the limit surface from failure.

Copyright © 2014 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2014)

DOI: 10.1002/nme
ALEATORY UNCERTAINTY AND SCALE EFFECTS IN COMPUTATIONAL DAMAGE MODELS

Figure 12. (a) Simulated dynamic indentation: transition from conical damage to radial cracking beneath

the cone, similar in character to observations in glass indentation [57, 60]. (b) Quasi-static indentation:

symmetry-breaking radial cracks beneath symmetry-preserving ring cracks [61, 63].

5. MODEL CALIBRATION OF VARIABILITY AND SCALE EFFECTS

Methods to account for strength data from differently sized experimental samples and methods to fit

repeated strength measurements to a Weibull (or similar size-dependent) distribution are discussed

in this section. These techniques provide median strength values and distribution properties (such

as the mean and standard deviation) for a fixed experimental configuration.

5.1. Normalizing experimental data to a standard reference volume

Ideally, the same size sample would be used in a variety of strength tests for various loading modes

such as triaxial compression, simple shear, buttonhead tension, and so on. If median strengths are

measured in these different loading modes using the same sample size, then the results directly map

to a strength limit surface in stress space. However, if different loading modes are measured using

samples of different sizes, then all of the data must be scaled to a common reference volume before

mapping the median strength surface.

Weibull scaling in Equation (2), including the reformulated versions in Equation (7), and

Equation (11), requires the volume to be expressed relative to a reference volume V . Similarly, the

median stress corresponds to measurements using a sample of size V . Recognizing that different

experiments might be performed using samples of different volumes, this section briefly explains

how all data from homogeneous loading experiments can be scaled to a common reference volume.

The far more subtle problem of scaling data from heterogeneous loading experiments (such as

Brazilian or bending tests) remains unsolved, but the relevant issues are introduced in Section 6.

Consider a set of experiments that measure a median strength  and Weibull modulus m for

a sample volume V . According to Equation (2), the same set of experiments conducted using a

sample with the reference volume V would produce the same Weibull modulus .m D m/ and a

volume-scaled median strength given by

D 

"

V

V 

#1=m

: (10)

By using this formula, experiments that use different specimen sizes may all be standardized to a

single reference volume V . The choice for the reference volume is arbitrary, but a reference volume is essential to generate failure threshold realizations by comparing the finite-element size to

the reference volume size used in laboratory calibration experiments. In geomechanics, it is common to employ two different sets of material properties, one for the so-called laboratory scale and

another for “field scale” (cf. [66]). The outlined implementation of Weibull theory in Equation (10)

automates scale effects so that only one material data set spans a large range of size scales.

The directionally-Weibull theory developed in Section 4 produces a different apparent (and

approximate) Weibull modulus m for different loading trajectories, which should be accounted for

when scaling sample sizes. The Weibull modulus for a selected model and loading direction can be

determined by numerically (or analytically) generating a series of realizations of the limit surface,

Copyright © 2014 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2014)

DOI: 10.1002/nme
O. E. STRACK, R. B. LEAVY AND R. M. BRANNON

collecting the corresponding CDF (failure probability) values along the fixed loading trajectory, and

then applying the analysis outlined below.

5.2. Measuring strength parameters

The following procedure for measuring the strength parameters is performed under two idealized

conditions. One is that the loading is statistically homogeneous with no significant macroscale stress

gradients prior to failure (e.g., quasi-static triaxial compression and extension and buttonhead tension tests). The second condition presumes that the same reference specimen size V is used in all

experiments. Even though these conditions can be difficult to meet, understanding the method used

to calibrate data under these idealizations is a prerequisite to subsequent discussions in Section 6

that release these assumptions by accounting for stress gradients within the specimen.

Although expanded testing and analysis are likely to be ongoing for many years, preliminary

scale-normalized strength data for SiC-N are shown in Figure 13 for a variety of pressures, with

median scaled data fitted to Kayenta’s limit surface function. Without scaling all data to a common

specimen size, no fitting exercise would be meaningful. The raw data measured using a sample of

volume V were scaled to a common reference volume V via Equation (10). For experiments having

a nonuniform stress field, the specimen volume V is understood to be the effective volume from a

data delocalization scheme (Section 6).

A suite of experiments using different loading directions must be run to identify the median

failure envelope, as it is loading path-dependent. This conventional part of the calibration process

identifies the dependence of on stress triaxiality (or, more generally, on N). Once the loading path

for a particular experiment is selected, the experiment must be conducted multiple times to obtain

a statistically significant collection of strength realizations 11; 2;:::;no. Deciding if the data fit

a Weibull distribution requires assigning values to the reference volume V , the median strength 

at the reference volume V , and the Weibull modulus m. Specifically, the data must be fit to the

“linearized” version of the Weibull equation

y D mx C b where y D ln ln 1

Psafe

 ln ln.2/; x D ln



; and b D ln

V

V ; (11)

which is obtained by taking the natural logarithm of Equation (2) twice and rearranging. Any data,

whether Weibull distributed or not, may be depicted in a plot of y versus x. Such data are proved to

be Weibull distributed if the data fall on a straight line with respect to these Weibull axes, though we

caution that even exactly Weibull-distributed data will exhibit stair-stepped deviations from the line

as a result of finite-sampling errors (Figures 14 and 16). As seen in Figure 16, finite sampling errors

in exactly Weibull-distributed data tend to produce a “dip” away from the lower-left portion of the

theoretical Weibull line. On the other hand, simple microphysically derived strength CDFs (cf. [31])

also have a dip that is not the result of finite sampling error. Resolving whether a dip is a genuine

(a) Data clouds and experimental data scaled to a large

sample reference volume V. 

(b) Data clouds and notional experimental data scaled to

a small sample reference volume V.

Figure 13. Examples of strength scaling and data clouds for comparison with experiments of silicon carbide.

a) Data clouds and experimental data scaled to a large sample reference volume V ; (b) data clouds and

notional experimental data scaled to a small sample reference volume V .

Copyright © 2014 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2014)

DOI: 10.1002/nme
ALEATORY UNCERTAINTY AND SCALE EFFECTS IN COMPUTATIONAL DAMAGE MODELS

(a) 600 points (b) 6000 points

Figure 14. Verification of mesh independence for the onset of failure. The figures show simulated .xk; yk/

data from Equation (13) for three mesh resolutions to demonstrate that failure statistics for the onset of

failure are mesh insensitive and approach the Weibull distribution as the number of realizations is increased.

(a) 600 points and (b) 6000 points.

attribute of the failure statistics or merely a byproduct of finite sampling is a strong motivation to

increase the number of strength tests in laboratory testing.

As discussed in Sections 2.2 and 4, “strength” is regarded as the value of the loading parameter

along a loading trajectory at failure. Here, “failure” is the initiation of stiffness and strength

reduction upon attainment of a peak value of , but the definition could be altered to mean phase

transformation, or any other “event” marking a distinct irreversible change in material behavior.

Repeating the experiment n times will give a set of strength realizations that can be put into canonical order, 1 6 2 6 :::k  6 n. Associated cumulative probabilities for sample number k are

given by

Pfail

k D k  1=2

n

and Psafe

k D 1  Pfailk : (12)

This estimator is advocated by Sullivan and Lauzon [67], who also recommend using no fewer than

20 measurements and preferably more than 50. An even larger number of measurements is required

to determine if the Weibull distribution is, in fact, the actual distribution. For example, more than

1000 tests would be required to distinguish between the Weibull distribution and the log-normal

distribution [68]. Other possibilities such as the Gumbel distribution [69] undoubtedly require an

equally large sample to validate. Consequently, no part of this paper presupposes that a Weibull

distribution applies to measured strength distributions. A Weibull distribution is used primarily as

a familiar example to illustrate considerations in data fitting. As explained in Section 4, the empirical framework naturally predicts strengths fitting no standard distribution, Weibull or otherwise.

Minimally, distributions must be used that will generate physically admissible realizations.

For fitting data to a Weibull distribution, the volumes V and V are both set equal to the specimen

volume used in the experiment, making b D 0 in Equation (11). From the collection of k; Psafe

k



data pairs, is identified to be the median (i.e., the value corresponding to Psafe D 1=2). Thus, only

the Weibull modulus m in Equation (11) remains yet to be inferred from the measured data. To obtain

m, the .k; Psafe

k / data pairs are converted to .xk; yk/ data pairs by the following transformations:

xk D ln

k



and yk D ln ln

1

Psafe

k

 ln ln 2: (13)

The Weibull modulus m is then determined via true linear (not affine) regression fitting of these

data pairs to y D mx. Namely,

m D

Pn

kD1

ynxn

Pn

kD1

x2

n

(14)

Copyright © 2014 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2014)

DOI: 10.1002/nme
O. E. STRACK, R. B. LEAVY AND R. M. BRANNON

Naturally, the regression error is a direct quantitative indicator of the extent to which a Weibull

distribution is appropriate. Sullivan and Lauzon [67] offer an alternative to simple linear regression.

However, by analyzing synthetically generated exactly Weibull-distributed data, we find their

method to be superior to linear regression only when the total number of strength measurements far

exceeds their recommended minimum of 20 data points. For data sets of fewer than 14 points, we

similarly find that Equation (14) more accurately predicts the Weibull modulus than a scheme that

sets the two (shape and scale) Weibull parameters by matching the sample mean and standard deviation to Equations 5 and 6. Even more than 14 points are needed for this alternative fitting scheme

to be superior whenever data beyond 2 or 3 standard deviations from the mean are misidentified

as “outliers”. For a typical engineering sample, which rarely exceeds 20 data points and is often

far below what is needed for statistical significance, these finite-sampling errors and inadvertent

removal of outliers seem to be best managed using the standard CDF estimator in Equation (12) and

the simple linear regression in Equation (14). Overall, the median is the most convenient measure

of central tendency because our approach aims to identify generally non-Weibull quantile surfaces

(for which the median is the 50% quantile).

Figure 14 verifies mesh independence (and recovery of the user-prescribed Weibull strength distribution) for the onset of failure under idealized spherical tension. The effect of finite sampling on a

Weibull plot (i.e., a plot of y vs. x) is apparent in the numerical data; finite sampling tends to cause

dips away from the Weibull line near the tail of the distribution (see also Figure 16).

In Figure 15(a), a ceramic flyer (the darker color on right side of both images) impacts a ceramic

target twice the thickness of the flyer. The damaged (spalled) material has been hidden from the

visualization of the simulation on the left, and the velocities are shown on the right (blue to red

color scale indicates low to high velocities, respectively). Spatial variability induces non-uniform

velocity on the surface of the target that would not occur in a deterministic simulation. Spall is

a dynamic tensile failure produced by the interaction of two rarefaction release waves, making it

a truly internal failure mode that is unaffected by surface smoothness at the specimen boundary.

Because of the strongly non-uniform stress fields occurring in spall tests, the data delocalization

methods discussed in Section 6 must be applied to iterate on the specimen volume associated with

measured spall strengths.

In Figure 16, Weibull plots and plots of the complementary CDF Psafe are shown for two different

experiments: the Hugoniot elastic limit strengths for B4C ceramic [54] and previously unpublished

tensile strengths for SiC-N for the Brazilian diametrical compression test. Both appear to be Weibulldistributed, but the B4C shock compression data clearly exhibit much lower variability than the

Brazilian data. These data deviate from the best Weibull fit but only by amounts well within the

range of deviations of Weibull data pulled from the exact analytical Weibull distribution for the same

modulus m (also shown in the figure). As seen, finite-sampling deviations in the exactly Weibulldistributed data are similar in magnitude and character to deviations of the measured data from the

Weibull fit. Thus, measured data may be closer to a Weibull distribution than might seem evident

from a Weibull plot.

(a) Spall test

SIMULATION

SIMULATION

(b) Brazillian test

Figure 15. Simulations with variability demonstrating bifurcation from nominal symmetries. (a) Spall test

and (b) Brazilian test.

Copyright © 2014 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2014)

DOI: 10.1002/nme
ALEATORY UNCERTAINTY AND SCALE EFFECTS IN COMPUTATIONAL DAMAGE MODELS

y

y

Figure 16. Measured strength data (thick black line) and best Weibull fit (thick straight gray line), along

with the same number of exactly Weibull distributed simulated data points (thin colored lines).

(a) Brazilian experiments (b) Brazilian simulations

Figure 17. (a) Laboratory strength data show effects of three specimen sizes (five strengths measured for

the largest specimen are near the left-most linear regression fit, twenty strengths for the medium-sized

specimen are in the center, and six for the smallest specimen are on the right). (b) Simulations for the

same sample sizes, with initial spherical strength perturbed using a Weibull modulus of 6, update those

in [70].

Finite-sampling effects (the “dipping” deviation from the Weibull line) are especially apparent

in the Brazilian strength data plotted in Figure 17, because those experiments acquired fewer data

points than the experiments in the left side Hugoniot elastic limit data of Figure 16. Examples of

Brazilian experiments and comparable simulations are in Figure 15(b).

The strength size scaling of the Brazilian material samples is noted in both the experiments and

simulations in Figure 17. These plots omit the typical normalization of the abscissa by the median

stress at failure so that the relative strengths and slopes can more readily be seen for the different

sample sizes. An interesting aspect of Figure 17 is that the experimental data appear to show not only

an increase in strength with decreasing specimen size but also a decreasing Weibull modulus (i.e.,

decreasing slope and hence greater relative variability in strength) with decreasing specimen size.

However, because of the small number of experimental samples and the impact of finite-sampling

errors, this trend is inconclusive. The numerical simulations in Figure 17 do not show the same

Copyright © 2014 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2014)

DOI: 10.1002/nme
O. E. STRACK, R. B. LEAVY AND R. M. BRANNON

trend, and further work is warranted to understand the scaling of the apparent Weibull modulus as a

function of sample size for the Brazilian test (both experimentally and numerically).

Although it would be possible to tune the parameters of the model to better match the data, doing

so would be grossly misleading, because the simulations failed to exhibit the same mitigation of

mesh sensitivity that was seen in the dynamic indentation simulations of Figure 1 (i.e., changing

the mesh would require inadmissible re-tuning of material properties). In the Brazilian simulations,

the methods of Section 6 have been confirmed to provide a high rate of convergence for the onset

of failure in the Brazilian test. However, further research breakthroughs are required to likewise

identify regularizations for the subsequent progression of failure (both at the element and structural

level) while also predicting realistic strength distributions and scale effects.

6. NONLOCAL CORRECTIONS OF DISCRETIZATION ERRORS

It has long been asserted that finite-element analysis of distributed damage requires some type of

nonlocal theory [71, 72]. Nonlocality is attributed to many physical sources such as crack interactions [73], which might be important on fine discretizations where aleatory uncertainty in crack

locations implies the need for correlated strength distributions. However, recent work [32] suggests

purely numerical motivations for nonlocality.

When attempting to calibrate scale-dependent strength from laboratory experiments, so-called

data delocalization techniques are recommended to assign the corresponding specimen volume. If

the stress field is uniform, then the specimen volume is the actual volume. However, for bending, Brazilian indirect tension, spall, and other experiments having a nonuniform stress field,

data-delocalization theory [38, 74–76] sets the specimen volume using a stress-weighted average. For a uniformly loaded and statistically homogeneous specimen, we adopt the axiom that

the probability of a local failure must scale with volume, not area (i.e., with V , not V 2=3). However, our own experimental data for the Brazilian indirect tension tests in Figure 17 scale with

area, not volume. This apparent area-based scaling is a natural result when the data are analyzed

using data-delocalization theory. Specifically, whenever stress is nonuniform (so that only a portion of the actual volume V is loaded to an appreciable stress level), the “effective volume” is

defined by

Veff D

Z

V

.=peak/

mdV (15)

where is a scalar measure of stress (such as the maximum principal stress) and peak is the

maximum of in the domain V . In the Brazilian test, for example, higher stresses are located over a

narrow region at the center of the domain, making this region have higher weight in Equation (15),

thus predicting that Veff scales with area as observed in the experiments.

Of course, evaluating Veff requires a description of the non-uniform stress field prior to failure.

Iterations are also required to find Veff because it depends on one of the properties being measured

(the Weibull modulus, m). The explicit presence of the Weibull modulus corresponds to an implicit

assumption that the strength is Weibull distributed. Furthermore, because Equation (15) depends on

a single scalar stress invariant, it fails to account for the tensorial nature of stress, and it implicitly

introduces a failure criterion into a situation where the failure criterion is part of the information

being sought in the experiment itself. To alleviate some of these shortcomings, we expect that a

practical generalization of Equation (15) might be

Veff D

Z

V

X.T/

X.T/peak

dV (16)

where the form and function of X.T/ should be regarded as the primary information sought in laboratory testing. As explained in an idealized context, namely Equation (11) of [31], the generalized

function “X” is expected to be a function of the loading mode and a functional of underlying

micromorphology (e.g., evolving flaw size and orientation statistics).

Copyright © 2014 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2014)

DOI: 10.1002/nme
ALEATORY UNCERTAINTY AND SCALE EFFECTS IN COMPUTATIONAL DAMAGE MODELS

Data-delocalization is required to interpret data from experiments having a nonuniform stress

field to assign a volume-dependent strength under uniform stress. Kamojjala and Brannon [32] have

developed a “data-relocalization” technique to reduce errors due to under-resolution of stress fields.

They explain that, to determine the failure initiation probability, the effective uniform stress acting

on a finite element should be set as an m-norm of the actual pre-failure stress field:

Psafe D 2V

V . eff

/

m

where eff 

” m dV

V

1=m

(17)

Here, m is the Weibull modulus, implying that this theory would require revision for nonWeibull distributions, perhaps replacing m with a micro-physically based X.T/ function, as in

the delocalization proposal of Equation (16). Also, is the pre-failure analytical principal stress

invariant field over the finite element, which is known for the Brazilian test and the other simple

cases studied by Kamojjala and Brannon. To adapt this relocalization formula to cases for which the

analytical stress field is unknown, a nonlocality algorithm is required to estimate the stress concentrations within an element on the basis of stresses in surrounding elements. This recommendation is

a purely numerical motivation for nonlocality, and it might be accommodated in a “semi-local” way

(i.e., without needing to query fields in neighboring elements) by applying an m-norm to the nonconstant stress field available within higher-order elements (Kamojjala and Brannon’s exploration

of higher-order elements essentially missed this opportunity by treating the stress field as constant

in the neighborhood of the Gauss points).

The previous sections have focused on statistics for the onset of failure, which controls the

statistically variable strength assigned at initialization and data-relocalization methods to adjust

the strength of an element in response to stress gradients. In a rigorous analysis of quasistatic

mode–I opening of a sharp crack, Bažant [23] addresses additional considerations associated with

failure progression that may or may not lead to cascading structural failure after failure initiation.

Specifically, Bažant pointed out that stress singularities near a crack tip increase more rapidly with

mesh refinement than element strengthening associated with the size effect, thus encouraging crack

growth. The related topic of damage progression is the subject of the next section.

7. STRENGTH AND STIFFNESS DEGRADATION

The preceding sections focused almost exclusively on the need to assign different values of material

strength to each element in a computational domain. The strength assigned to different elements

is statistically variable and size-scale dependent, consistent with laboratory observations. This

variability in strength provides perturbations in the stress field necessary to stimulate inherent

instabilities such as bifurcation of axisymmetric loading into radial cracking. However, once cascading failure has begun, additional scale effects are well known to be necessary to mitigate mesh

dependence during the progression of failure [77].

Seminal cohesive zone formulations [78, 79] and more recent models (cf. [30, 80, 81]) introduce

a scale effect by requiring a certain amount of energy be required to advance a crack. Until recently,

cohesive zone formulations have been developed on the basis of an implicit assumption that a finite

domain contains exactly one localized zone. However, high-rate loading can produce a network of

actively growing cracks within a domain, therefore calling for an increase in required fracture energy

[58, 77]. As explained in this section, this overall effect is accommodated in our model by requiring

a finite time to elapse before an element can lose strength.

To simulate crack growth, conventional isotropic damage models usually employ a scalar damage

parameter, D, that varies from 0 for undamaged material to 1 for fully failed material [44]. These

models lead to a drop in both stiffness and strength with increasing strain (i.e., softening) and

accordingly induce localization of material response. Without such localization, the governing

equations can change type, potentially requiring solution of a diffusion equation [82]. On the other

hand, traditional viscous rate dependence in the governing equations can have a regularizing effect

[83]. Similarly, the extreme localization of brittle and quasi-brittle materials (i.e., formation of

macrocracks) allows the governing equations to remain hyperbolic as long as boundary conditions

Copyright © 2014 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2014)

DOI: 10.1002/nme
O. E. STRACK, R. B. LEAVY AND R. M. BRANNON

at the new surfaces are handled properly. In a smeared damage model, these boundary conditions

are approximated by a boundary layer of failed material that offers resistance to compression but

progressively reduced resistance to tension and ultimately only frictional resistance to shear. Further discourse about the character of conventional damage models is beyond the scope of this paper,

where our focus here is on seeding these models realistically with uncertainty and scale effects. The

conventional damage model used for our case study (Kayenta) is well documented [6], so we merely

call attention to some aspects that differ from similar models available in commercial codes.

The Kayenta damage model employs a distinctive approach to degrading strength. Taking a cue

from experimental observations [84] that cracks in high-stress environments grow at approximately

constant speed (the effect of which is seen in Figure 18), the model smoothly decreases strength

and elastic stiffness to failed values over a scale-dependent time interval (conceptually similar to

the method applied in [85]). Consistent with the notion that element failure corresponds to traversal

of a macroscale “train” of coalescing microcracks across the element, a small element fails sooner

than a large element, and hence, the model keys the “time-to-failure” to the element size. For simplicity, elements with large aspect ratios are avoided so that the time-to-failure has a length scale

of approximately V 1=3, where V is the volume of the finite element that is softening. Not only will

larger elements be initialized with lower median strength, this straightforward (and computationally efficient) method for scaling the time-to-failure ensures that larger elements might still require

greater fracture energy, depending on how the material parameters are set. Implementation details,

including formulas for degrading strength and stiffness with damage, are provided in [6].

The time-to-failure scaling is insensitive to the loading rate. Accordingly, a specimen loaded

slowly will seem to lose strength suddenly because the time-to-failure will be very small relative to

the total loading time. If the same specimen is loaded very rapidly to the same strain level, it will

lose strength over a larger fraction of the total loading interval (because a fixed amount of time must

pass before the stiffness and strength degradations fully take effect). Thus, when loaded rapidly,

our model will appear to have a finite plastic strain-to-failure (which is similar in character to the

Johnson–Holmquist class of ceramic models [45, 48–50] whenever those models are compared with

ours at a given loading rate, but a key distinction is that our model’s use of time-to-failure provides

different apparent strain-to-failure at different loading rates).

Figure 18. Laboratory observations of the effect of impact speed on the number of radial cracks [70]. Shown

are B4C samples impacted at speeds ranging from 100 to 400 m/s from left to right.

Figure 19. Realistic trends in damage patterns when including spatial variability, scale effects, and time-tofailure. Changing the WC impactor speed (from 200 to 500 m/s from left to right, using the same mesh and

property seeding) produces a corresponding increased number of radial cracks in simulations of SiC targets.

Copyright © 2014 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2014)

DOI: 10.1002/nme
ALEATORY UNCERTAINTY AND SCALE EFFECTS IN COMPUTATIONAL DAMAGE MODELS

Because our model requires a fixed amount of time to pass before loss of strength occurs, the

higher strain rate caused by higher impact speed allows more material to reach the failure threshold

before release waves (or Mott waves—see also the related “event horizon” [76]) can arrive from the

radial cracks that formed slightly sooner at weaker points. The well-known experimental observation

that the number of radial cracks in axisymmetric loading increases with loading rate [86, 87] (see,

for example, observations for B4C in Figure 18) is therefore easily predicted with this time-to-failure

softening algorithm because it essentially allows a larger number of flaws to become activated at

high rates (see similar trends in Figure 19 for simulations of SiC). Consistent with [58], this higher

probability of flaw activation at higher loading rates was also mentioned by Brannon and Gowen

[31] as being a simple case of a larger portion of the 3D stress state (as depicted in our Figure 4)

falling transiently farther above each crack’s size-dependent failure line, hence naturally activating

larger populations of both flaw sizes and orientations.

8. NECESSITY OF BOTH STRENGTH VARIABILITY AND TIME-TO-FAILURE SCALING

As seen in Figure 20, imposing both strength variability and time-to-failure scaling significantly

reduces mesh dependencies of damaged zones, crack patterns, and rate of crack propagation for the

dynamic indentation of ceramics. The images in Figure 20 depict simulations of a 6.35 mm diameter

tungsten-carbide sphere launched at 375 m/s into a 25.4 mm diameter SiC-N ceramic cylinder,

described in [7]. Each set of images in Figure 20(a)–(d) contains a low resolution simulation next to

successively higher resolutions. All simulation parameters were kept the same, except for the mesh

resolution and the feature being explored in the model (which was changed for each subset).

Figure 20(a) contains results for a conventional model (Kayenta) calibrated to the median of the

experimental data and ignoring size effects. Here, the high-pressure limit surface parameter Y was

calibrated to the data used in [49], the tensile strength parameter h was calibrated to buttonhead tension results in [88], and the shear strength parameter s was calibrated to triaxial compression results

in [11]. To create an objective comparison between different features of the model, the time-tofailure was computed for the average element size of the mesh in the medium resolution simulation

for the fully enhanced model in Figure 20(d) (the middle image) and applied as a constant time-tofailure for the conventional model. In Figure 20(a), the number of radial cracks and especially the

depth of penetration depend on the element size in the mesh. For the simulations in Figure 20(b), the

time-to-failure was scaled according to element size as discussed in Section 7, while other parameters were the same as in the conventional model in Figure 20(a). Note the reduction in the number

(a) Conventional model (no variability or scaling) (b) Time-to-failure scaling without variability

(c) Variability without time-to-failure scaling (d) Variability and time-to-failure scaling

Figure 20. Influence of strength variability and time-to-failure scaling in dynamic indentation. Subfigures

each show three results at the same moment in time on progressively finer mesh resolutions.

Copyright © 2014 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2014)

DOI: 10.1002/nme
O. E. STRACK, R. B. LEAVY AND R. M. BRANNON

of radial cracks in the high-resolution simulation on the right due to the lower time-to-failure for

the smaller elements in that mesh. For the simulations in Figure 20(c), a Weibull variability of

strength, including the volume scaling term in Equation (7), was applied to the h parameter in

the Kayenta model, as discussed in Section 4. The other parameters (including constant time-tofailure) were the same as in the conventional model in Figure 20(a). Note the significant effect of

including spatial variability. The size scaling effects have resulted in a significantly higher median

element strength (and hence more realistic ballistic performance) while the variability has, in principle (modulo issues with the progression of damage discussed in Section 7 that can be encountered

in simulations of calibration experiments), preserved the match of the model to the suite of calibration experiments discussed in Section 5. Also, using a size-independent time-to-failure causes the

horizontal lateral cracking (on the perimeter toward the bottom of the target) to be located at successively higher elevations as the mesh resolution increases (contrast this to the location of these

cracks in Figure 20(d)).

The model in Figure 20(d) includes both the spatial variability of strength discussed in

Figure 20(c) and the time-to-failure scaling discussed for Figure 20(b). The two effects work

together, resulting in a model with substantially reduced mesh dependence, both in the amount of

damage induced, and in the rate of growth of damage. The approximate number of radial cracks and

location of lateral cracking are similar between the meshes. Mesh dependence has been reduced, but

not eliminated. Significant work is still required for truly predictive ballistic simulations of ceramics, but spatial variability, size scaling, and time-to-failure scaling represent technical advances in

comparison with predecessor theories. Moreover, these features can be added without significant

increases in computational (or implementation/development) overhead.

9. CONCLUSIONS

Enhancing an existing conventional deterministic damage model to support spatial variability

requires (in principle) no revisions of the existing constitutive model’s source code. Code changes

are applied at the host-code level, where fields are allocated (during initialization) for spatially

varying material parameters. For example, each element in a simulation might have its own randomly assigned scale-dependent strength. These perturbed parameters are assigned values by

applying whatever statistical distribution is evident in laboratory testing of finite-sized specimens

but with care given to anticipate parameter correlations required to produce admissible material

realizations (e.g., to avoid inadvertently generating non-convex strength limit surfaces).

If a material parameter 	 happens to be observed in the laboratory to be Weibull distributed and if

it also has corresponding Weibull scale effects, then the spatially variable and size-dependent value

assigned to a finite element of volume V is found by computing a random number R and then saving

	 D 	

"

V ln R

V ln.1=2/#1=m

(18)

into a field array, which then provides each finite element with its own perturbed material property

	. Here, 	 is the median value of 	 observed in the laboratory using a specimen of volume V , and

m is the Weibull modulus quantifying observed variability in the property. Appropriate revisions

naturally apply if the parameter is not Weibull distributed.

Given the overarching goal to introduce variability into an already available deterministic damage model, each material property 	 that is already available in existing property databases is merely

re-interpreted as the median about which uncertainty and scale effects are imposed. Thus, the only

new inputs are V and (if Weibull-distributed) m, both of which can be determined by reinspecting the laboratory data originally used to determine 	. “Data delocalization” schemes [38, 74–76]

are needed to set the effective specimen volume V for any strength experiments involving nonuniform stress fields, whereas “relocalization” [32] essentially re-introduces stress variation that is

not resolved on low-order finite elements.

If the value of a property 	 was determined from only a single experiment, then the experiment

would need to be repeated enough times to quantify the scale-dependent central tendency (e.g.,

Copyright © 2014 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2014)

DOI: 10.1002/nme
ALEATORY UNCERTAINTY AND SCALE EFFECTS IN COMPUTATIONAL DAMAGE MODELS

median) and the aleatory uncertainty (e.g., standard deviation). At least 20 (and preferably more

than 50) such tests are recommended for each loading mode (shear, tension, compression, etc.), but

indirect methods, such as “line VISAR” in spall testing seem promising to quantify variability in a

single experiment [53, 89].

Overall, parameterizing material models from observed aleatory uncertainty and associated scale

effects shows promise for enhancing the predictive capability of numerical methods of fracture and

fragmentation while still keeping computational and development costs manageable. An ordinary

deterministic limit surface is effectively replaced with a family of level sets (called aleatory quantile surfaces) associated with a notional failure probability function. Various scalar-based strengths

(e.g., buttonhead tension, triaxial compression, spall, simple shear, Brazilian indirect tension, etc.)

each have different distributions, yet they are all simultaneously modeled by reinterpreting them

as directed derivatives of the single, unified tensor-based quantile function. The tensorial nature of

stress is therefore directly and conveniently accounted for, giving non-Weibull, pressure-dependent,

and generally loading-mode-dependent strength distributions based on laboratory evidence. A scaledependent time-to-failure algorithm for degradation of stiffness (and strength) is shown to not only

alleviate mesh sensitivity in dynamic indentation but it also naturally ensures activation of larger

populations of flaws during higher-rate loading, automatically giving more fragments with increased

loading rate.

ACKNOWLEDGEMENTS

In addition to administrative, technical, and financial support by Sandia laboratory managers (especially

Tom Pfeifle and Randy Summers) and Army Research Laboratory managers (especially Bill Bruchey, Scott

Schoenfeld, and Todd Bjerke), funding and support for this work provided by TACOM Project Manager

John Rowe is deeply appreciated. We are additionally grateful to the following individuals: Rich Becker (for

his 2002 unpublished study of the effect of strength perturbations in plane stress), Tracy Vogler and Lalit

Chhabildas (for providing data and constructive feedback), Richard Jensen and Josh Houskamp (for exercising the model in Brazilian simulations and beta-testing), Moo Lee and David Bronowski (for Brazilian

and Triaxial testing of SiC-N ceramics), Mike Veilleux (for Monte Carlo verification testing of the Weibull

perturbations), and Tim Fuller, Mike Wong, and Greg Sharp (for extensive testing and improvement of

the Kayenta model and its implementation in several host codes). The majority of this research was performed at Sandia National Laboratories and the Army Research Laboratory, with some additional support

from Schlumberger corporation. Sandia National Laboratories is a multi-program laboratory managed and

operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the US

Department of Energy’s National Nuclear Security Administration under contract DE-AC04-94AL85000.

SAND2014-0865J.

REFERENCES

1. Robinson AC, Brunner TA, Carroll S, Drake R, Garasi CJ, Gardiner T, Haill T, Hanshaw H, Hensinger D, Labreche

D, Lemke R, Love E, Luchini C, Mosso S, Niederhaus J, Ober CC, Petney S, Rider WJ, Scovazzi G, Strack OE,

Summers R, Trucano T, Weirs VG, Wong M. ALEGRA: an arbitrary Lagrangian–Eulerian multimaterial, multiphysics code. Proceedings of the 46th AIAA Aerospace Sciences Meeting, Vol. AIAA-2008-1235, Reno, NV, 2008;

1–39.

2. Flanagan DP, Belytschko T. A uniform strain hexahedron and quadrilateral with orthogonal hourglass control.

International Journal for Numerical Methods in Engineering 1981; 17:679–706.

3. Belytschko T, Tsay C. A stabilization procedure for the quadrilateral plate element with one point quadrature.

International Journal for Numerical Methods in Engineering 1983; 19:405–419.

4. Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. International Journal

for Numerical Methods in Engineering 1999; 46:131–150.

5. Fossum AF, Brannon RM. On a viscoplastic model for rocks with mechanism-dependent characteristic times. Acta

Geotechnica 2006; 1:89–106.

6. Fossum AF, Brannon RM, Strack OE. Kayenta: theory and user’s guide. Technical Report SAND2009-2282, Sandia

National Laboratories: Albuquerque, NM, 2009.

7. LaSalvia JC, Normandia MJ, Miller HT, MacKenzie DE. Sphere impact induced damage in ceramics: I. armor-grade

SiC and TiB2. Ceramic Engineering and Science Proceedings 2005; 26(7):171–181.

8. de St. Germain J, McCorquodale J, Parker S, Johnson C. Uintah: a massively parallel problem solving environment.

In Ninth IEEE International Symposium on High Performance and Distributed Computing. IEEE: Piscataway, NJ,

2000; 33–41.

Copyright © 2014 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2014)

DOI: 10.1002/nme
O. E. STRACK, R. B. LEAVY AND R. M. BRANNON

9. Mauge C, Kachanov M. Mechanics of anisotropic materials with multiple cracks. Key Engineering Materials 1996;

120-121:3–46.

10. Abramov AA, Ul’yanova VI, Yukhno LF. Certain methods for linear problems with inexact initial data. Journal of

Computational and Applied Mathematics 2006; 192:2–10.

11. Lee MY, Brannon RM, Bronowski DR. Uniaxial & triaxial compression tests of silicon carbide ceramics under

quasi-static loading condition. Technical Report SAND2004-6005, Sandia National Laboratories: Albuquerque, NM,

2004.

12. Mandel J. Propagation des surfaces de discontinuite dans un milieu elastoplastique. In Stress Waves in Anelastic

Solids, Kolsky H, Prager W (eds). Springer: Berlin, 1964; 331–340.

13. Rice JR. The localization of plastic deformation. In Theoretical and Applied Mechanics, Vol. 1, Koiter WT (ed.).

North-Holland: Delft, 1976; 207–220.

14. Holcomb DJ, Costin LS. Detecting damage surfaces in brittle materials using acoustic emissions. Transactions of the

ASME Journal of Applied Mechanics 1986; 53(3):536–544.

15. Dannemann KA, Chocron S, Nicholls AE, Anderson CE Jr. Compressive damage development in confined

borosilicate glass. Materials Science and Engineering A 2008; 478:340–350.

16. Brannon RM, Wells JM, Strack OE. Validating theories for brittle damage. Metallurgical and Materials Transactions

A 2007; 38A:2861–2868.

17. Leavy RB, Rickter B, Normandia MJ. Modeling dynamically impacted ceramic material experiments. Ceramic

Engineering and Science Proceedings 2005; 26(7):11–18.

18. Becker R. The effect of porosity distribution on ductile fracture. Journal of the Mechanics and Physics of Solids

1987; 35(5):577–599.

19. Becker R. An analysis of shear localization during bending of a polycrystalline sheet. Journal of Applied MechanicsTranasactions of the ASME 1992; 59(3):491–496.

20. Ray D, Flinders RM, Anderson A, Cutler RA, Campbell J, Adams JW. Effect of microstructure and mechanical properties on the ballistic performance of SiC-based ceramics. In Advances in Ceramic Armor II: Ceramic Engineering

and Science Proceedings. chap. 7. John Wiley and Sons, Inc., 2008; 85–96.

21. Grinfeld MA, Schoenfeld SE, Wright TW. Morphological instability of failure fronts. Applied Physics Letters 2006;

88(10):104102.1–104102.3.

22. Timmel M, Kolling S, Osterrieder P, Bois PD. A finite element model for impact simulation with laminated glass.

International Journal of Impact Engineering 2007; 34:1465–1478.

23. Bažant ZP. Fracture and Size Effect in Concrete and Other Quasibrittle Materials. CRC Press: Boca Raton, Boston,

London, New York, Washington D.C., 1998.

24. Bourne NK, Millett JCF, Chen M, McCauley JW, Dandekar D. On the Hugoniot elastic limit in polycrystalline

alumina. Journal of Applied Physics 2007; 102:073514–1–9.

25. Bažant ZP, Xi Y, Reid S. Statistical size effect in quasi-brittle structures: I. Is weibull theory applicable. Journal of

Engineering Mechanics 1991; 117(11):2609–2622.

26. Grady DE, Kipp ME. Geometric statistics and dynamic fragmentation. Journal of Applied Physics 1985; 58(3):

1210–1222.

27. Baer MR. Mesoscale modeling of shocks in heterogeneous reactive materials. In Shockwave Science and Technology

Reference Library, Horie Y (ed.). Springer: Berlin Heidelberg, 2007; 321–356.

28. Dienes JK, Zuo QH, Kershner JD. Impact initiation of explosives and propellants via statistical crack mechanics.

Journal of the Mechanics and Physics of Solids 2006; 54:1237–1275.

29. Bishop JE, Strack OE. A statistical method for verifying mesh convergence in monte carlo simulations with

application to fragmentation. International Journal for Numerical Methods in Engineering 2011; 88(3):279–306.

30. Zhou F, Molinari JF. Dynamic crack propagation with cohesive elements: a methodology to address mesh

dependency. International Journal for Numerical Methods in Engineering 2004; 59(1):1–24.

31. Brannon RM, Gowen TJ. Aleatory quantile surfaces in damage mechanics. Journal of the European Ceramic Society

2014; 34(11):2643–2653.

32. Kamojjala K, Brannon RM. Data relocalization to mitigate slow convergence caused by under-resolved stress fields

in computational damage mechanics. Journal of the European Ceramic Society 2014; 34(11):2723–2730.

33. Batdorf SB, Heinisch HL. Weakest link theory reformulated for arbitrary fracture criterion. Journal of the American

Ceramic Society 1978; 61(7-8):355–358.

34. Cupschalk SG. Probabilistic technique related to amplitude-dependent damping under combined stresses. ASTM

Special Technical Publication 1997; 1304:349–364.

35. Nemeth NN. Probability density distribution of the orientation of strength-controlling flaws from multiaxial loading

using the unit-sphere stochastic strength model for anisotropy. International Journal of Fracture 2014; 185(1-2):

97–114.

36. Macon DJ, Brannon RM, Strack OE. Plastic cap evolution law derived from induced transverse isotropy in dilatational triaxial compression. Technical Report SAND2014-1217, Sandia National Laboratories: Albuquerque, NM,

2014.

37. Fossum AF, Brannon RM. Unified compaction/dilatation, strain-rate sensitive constitutive model for rock mechanics

structural analysis applications. ARMA/NARMS 04-546, 6th North American Rock Mechanics Symposium, Houston,

Texas, June 5-9, 2004; 1–12.

38. Hild F, Denoual C, Forquin P, Brajer X. On the probabilisitc-deterministic transition involved in a fragmentation

process of brittle materials. Computers and Structures 2003; 81:1241–1253.

Copyright © 2014 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2014)

DOI: 10.1002/nme
ALEATORY UNCERTAINTY AND SCALE EFFECTS IN COMPUTATIONAL DAMAGE MODELS

39. Wellman GW. Failprob—a computer program to compute the probability of failure of a brittle component. Technical

Report SAND2002-0409, Sandia National Laboratories: Albuquerque, NM, 2002.

40. Nemeth NN, Powers LM, Janosik LA, Gyekenyesi JP. Cares/life ceramics analysis and reliability evaluation of

structures life prediction program. Technical Report NASA/TM 2003-106316, NASA: 7121 Standard Drive; Hanover,

MD 21076, 2003.

41. Grady DE. Local inertial effects in dynamic fragmentaion. Journal of Applied Physics 1982; 53:322–325.

42. Drugan WJ. Dynamic fragmentation of brittle materials: analytical mechanics-based models. Journal of the

Mechanics and Physics of Solids 2001; 49:1181–1208.

43. Meyer, HW Jr, Brannon RM. A model for statistical variation of fracture properties in a continuum mechanics code.

International Journal of Impact Engineering 2012; 42:48–58.

44. Brannon RM, Leelavanichkul S. Survey of four damage models for concrete. Technical Report SAND2009-5544,

Sandia National Laboratories: Albuquerque, NM, 2009.

45. Johnson GR, Holmquist TJ. A computational constitutive model for brittle materials subjected to large strains, high

strain rates and high pressures. In Proceedingedings of the 1990 International Conference on Shock-Wave and HighStrain-Rate Phenomena in Materials (EXPLOMET), Meyers MA, Murr LE, Staudhammer KP (eds). Marcel Dekker

Inc.: San Diego, 1992; 1075–1081.

46. Johnson GR, Cook WH. Fracture characteristics of three metals subjected to various strain rates, temperatures and

pressures. Engineering Fracture Mechanics 1985; 21(1):31–48.

47. Fossum AF, Brannon RM. The Sandia GeoModel theory and user’s guide. Technical Report SAND2004-3226, Sandia

National Laboratories: Albuquerque, NM, 2004.

48. Johnson GR, Holmquist TJ. An improved computational constitutive model for brittle materials. In High-Pressure

Science and Technology 1993, Schmidt SC, Schaner JW, Samar GA (eds). American Institute of Physics: New York,

1994; 981–984.

49. Holmquist TJ, Johnson GR. Response of silicon carbide to high velocity impact. Journal of Applied Physics 2002;

91(9):5858–5866.

50. Johnson GR, Holmquist TJ, Beissel SR. Response of aluminum nitride (including a phase change) to large strains,

high strain rates, and high pressures. Journal of Applied Physics 2003; 94(5):1639–1646.

51. Deshpande VS, Evans AG. Inelastic deformation and energy dissipation in ceramics: a mechanism-based constitutive

model. Journal of the Mechanics and Physics of Solids 2008; 56(10):3077–3100.

52. Cho K, Katz R, Bar-On I. Uniaxial strength and mode I fracture toughness of hot pressed silicon carbide (SiC)

materials. Technical Report ARL-TR-1406, Army Research Laboratory: Aberdeen Proving Ground, MD, 1997.

53. Vogler TJ, Trott WM, Reinhart WD, Alexander MD, Furnish MD, Knudson LC. Using the line-visar to study

multi-dimensional and mesoscale impact phenomena. International Journal of Impact Engineering 2008; 35(12):

1844–1852.

54. Vogler TJ, Reinhart WD, Chhabildas LC. Dynamic behavior of boron carbide. Journal of Applied Physics 2004;

95(8):4173–4183.

55. Wright TW, Ramesh KT. Dynamic void nucleation and growth in solids: a self-consistent statistical theory. Journal

of the Mechanics and Physics of Solids 2008; 56:336–359.

56. Oh KP. A Monte Carlo study of the strength of unidirectional fiber-reinforced composites. Journal of Composite

Materials 1979; 13:311.

57. Ford C, Bush MB, Hu X-Z, Zhao H. Numerical interpretation of cone crack initiation trends in a brittle coating on a

compliant substrate. Materials Science and Engineering A 2004; 380:137–142.

58. Kimberley J, Ramesh K, Daphalapurkar N. A scaling law for the dynamic strength of brittle solids. Acta Materialia

2013; 61:3509–3521.

59. Frank FC, Lawn BR. On the theory of Hertzian fracture. Proceedings of the Royal Society of London. Series A,

Mathematical and Physical Sciences 1967; 299(1458):291–306.

60. Rhee Y-W, Kim H-W, Deng Y, Lawn BR. Contact-induced damage in ceramic coatings on compliant substrates:

fracture mechanics and design. Journal of the American Ceramic Society 2001; 84(5):1066–1072.

61. Wereszczak A, Johanns K, Kirkland TP, Anderson CE, Jr., Behner T, Patel P, Templeton W. Strength and contact

damage responses in a soda-lime-silicate and a borosilicate glass. Proceedings of the 25th Army Science Conference,

Orlando, FL, 2006; 1–8. Paper FP-05.

62. Wereszczak AA, Kirkland TP, Strong KT Jr., Jadaan OM, Thompson GA. Size-scaling of tensile failure stress in

boron carbide. Advances in Applied Ceramics 2010; 109(8):487–492.

63. Chai H, Lawn BR, Wuttiphan S. Fracture modes in brittle coatings with large interlayer modulus mismatch. Journal

of Materials Research 1999; 14(9):3805–3817.

64. Graham-Brady L. Statistical characterization of meso-scale uniaxial compressive strength in brittle materials with

rarandom occurring flaws. International Journal of Solids and Structures 2010; 47:2398–2413.

65. Daphalapurkar N. Predicting variability in the dynamic strength of brittle materials considering pre-existing flaws.

Journal of the Mechanics and Physics of Solids 2011; 59:297–319.

66. Jaeger JC, Cook NGW. Fundamentals of Rock Mechanics (Second edn). John Wiley & Sons, Inc.: New York, 1976.

67. Sullivan JD, Lauzon PH. Experimental probability estimators for Weibull plots. Journal of Material Science Letters

1986; 5:1245–1247.

68. Weibull W. A statistical distribution function of wide applicability. Transactions of the ASME Journal of Applied

Mechanics 1951; 18:293–297.

69. Gumbel EJ. Statistics of Extremes. Columbia University Press: New York, 1957.

Copyright © 2014 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2014)

DOI: 10.1002/nme
O. E. STRACK, R. B. LEAVY AND R. M. BRANNON

70. Leavy RB, Brannon RM, Strack OE. The use of sphere indentation experiments to characterize ceramic damage

models. International Journal of Applied Ceramic Technology 2010; 7(5):606–615.

71. Bažant ZP. Imbricate continuum and its variational derivation. Journal of Engineering Mechanics, ASCE 1984;

110(12):1693–1712.

72. Bažant ZP. Mechanics of distributed cracking. Applied Mechanics Reviews 1986; 39:675–705.

73. Bažant ZP. Nonlocal damage theory based on micromechanics of crack interactions. Journal of Engineering

Mechanics 1994; 120(3):593–617.

74. Davies DGS. The statistical approach to engineering design in ceramics. Proceedings of the British Ceramic Society

1973; 22:429–452.

75. Bažant ZP, Xiang Y, Adley MD, Prat PC, Akers SA. Microplane model for concrete: II: Data delocalization and

verification. Journal of Engineering Mechanics 1996; 122(3):255–262.

76. Forquin P, Denoual C, Cottenot CE, Hild F. Experiments and modelling of the compressive behaviour of two SiC

ceramics. Mechanics of Materials 2003; 35(10):987–1002.

77. Molinari JF, Gazonas G, Raghupathy R, Rusinek A, Zhou F. The cohesive element approach to dynamic fragmentation: the question of energy convergence. International Journal for Numerical Methods in Engineering 2007;

69(3):484–503.

78. Dugdale DS. Yielding of steel sheets containing slits. Journal of the Mechanics and Physics of Solids 1960; 8:

100–104.

79. Barenblatt GI. The mathematical theory of equilibrium of cracks in brittle fracture. Advances in Applied Mechanics

1962; 7:55–129.

80. Elices M, Guinea GV, Gómez J, Planas J. The cohesive zone model: advantages, limitations, and challenges.

Engineering Fracture Mechanics 2002; 69:137–163.

81. de Borst R, Gutiérrez MA, Wells GN, Remmers JJC, Askes H. Cohesive-zone models, higher-order continuum

theories and reliability methods for computational failure analysis. International Journal for Numerical Methods in

Engineering 2004; 60:289–315.

82. Chen Z. Simulating the evolution of localization based on the diffusion of damage. International Journal of Solids

and Structures 2000; 37:7465–7479.

83. Simo J, Hughes T. Compuational Inelasticity (corrected edition). Springer, 2000.

84. Fineberg J, Marder M. Instability in dynamic fracture. Physics Reports 1999; 313:1–108. Table 3.

85. Becker R. Nonlocal control of fracture propagation in numerical simulations. Technical Report ARL-MR-842, Army

Research Laboratory: Aberdeen Proving Ground, MD, 2013.

86. Mott NF. Fragmentation of shell cases. Proceedings of the Royal Society of London. Series A, Mathematical and

Physical Sciences 1947; 189(1018):300–308.

87. Denoual C, Hild F. Dynamic fragmentation of brittle solids: a multi-scale model. European Journal of Mechanics -

A/Solids 2002; 21:105–120.

88. Wereszczak AA, Kirkland TP, Strong KT, Campbell J, LaSalvia JC, Miller HT. Size-scaling of tensile failure stress

in a hot-pressed silicon carbide. International Journal of Applied Ceramic Technology 2010; 7(5):635–642.

89. Furnish M, Gray I GT, Bingert J. Line visar and post-shot metallography comparisons for spall analysis. In Dynamic

Behavior of Materials, vol. 1, Proulx T (ed.)., Conference Proceedings of the Society for Experimental Mechanics

Series. Springer: New York, 2011; 351–352.

Copyright © 2014 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2014)

DOI: 10.1002/nme
29

1 DEFENSE TECHNICAL

(PDF) INFORMATION CTR

DTIC OCA

2 DIRECTOR

(PDF) US ARMY RESEARCH LAB

RDRL CIO LL

IMAL HRA MAIL & RECORDS MGMT

1 GOVT PRINTG OFC

(PDF) A MALHOTRA

1 PM ABRAMS 

(PDF) SFAE GCS HBCT

J ROWE

39 RDRL CIH C

(PDF) D GROVE

RDRL WM

P BAKER

J MCCAULEY

RDRL WMM B

B AYDELLOTE

B LOVE

RDRL WML H

G GAZONAS

B SCHUSTER

RDRL WMM E

J LASALVIA

J SWAB

RDRL WMP

D LYON

S SCHOENFELD

RDRL WMP B

C HOPPEL

S SATAPATHY

A SOKOLOW

T WEERASOORIYA

RDRL WMP C

R BECKER

S BILYK

T BJERKE

D CASEM

J CLAYTON 

D DANDEKAR

M GREENFIELD

R LEAVY

J LLOYD

M RAFTENBERG

C WILLIAMS

RDRL WMP D

R DONEY

B SCOTT

M ZELLNER

RDRL WMP E

S BARTUS

B CHAMISH

D HACKBARTH

J HOUSKAMP

E HORWATH

C KRAUTHAUSER

P SWOBODA

RDRL WMP F

D FOX

N GNIAZDOWSKI

RDRL WMP G

J STEWART
30

INTENTIONALLY LEFT BLANK.