NASA/TM—2013-217810

September 2013

An Erratum was added to this report March 2014.

Noel N. Nemeth

Glenn Research Center, Cleveland, Ohio

Unit-Sphere Anisotropic Multiaxial StochasticStrength Model Probability Density Distribution

for the Orientation of Critical Flaws
NASA STI Program . . . in Profi le

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientifi c and Technical Information (STI) 

program plays a key part in helping NASA maintain 

this important role.

The NASA STI Program operates under the auspices 

of the Agency Chief Information Offi cer. It collects, 

organizes, provides for archiving, and disseminates 

NASA’s STI. The NASA STI program provides access 

to the NASA Aeronautics and Space Database and 

its public interface, the NASA Technical Reports 

Server, thus providing one of the largest collections 

of aeronautical and space science STI in the world. 

Results are published in both non-NASA channels 

and by NASA in the NASA STI Report Series, which 

includes the following report types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major signifi cant phase 

of research that present the results of NASA 

programs and include extensive data or theoretical 

analysis. Includes compilations of signifi cant 

scientifi c and technical data and information 

deemed to be of continuing reference value. 

NASA counterpart of peer-reviewed formal 

professional papers but has less stringent 

limitations on manuscript length and extent of 

graphic presentations.

• TECHNICAL MEMORANDUM. Scientifi c 

and technical fi ndings that are preliminary or 

of specialized interest, e.g., quick release 

reports, working papers, and bibliographies that 

contain minimal annotation. Does not contain 

extensive analysis.

• CONTRACTOR REPORT. Scientifi c and 

technical fi ndings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientifi c and technical 

conferences, symposia, seminars, or other 

meetings sponsored or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientifi c, 

technical, or historical information from 

NASA programs, projects, and missions, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientifi c and 

technical material pertinent to NASA’s mission.

Specialized services also include creating custom 

thesauri, building customized databases, organizing 

and publishing research results.

For more information about the NASA STI 

program, see the following:

• Access the NASA STI program home page at 

http://www.sti.nasa.gov

• E-mail your question to help@sti.nasa.gov

• Fax your question to the NASA STI 

Information Desk at 443–757–5803

• Phone the NASA STI Information Desk at

 443–757–5802

• Write to:

 STI Information Desk

 NASA Center for AeroSpace Information

 7115 Standard Drive

 Hanover, MD 21076–1320
NASA/TM—2013-217810

September 2013

An Erratum was added to this report March 2014.

Noel N. Nemeth

Glenn Research Center, Cleveland, Ohio

Unit-Sphere Anisotropic Multiaxial StochasticStrength Model Probability Density Distribution

for the Orientation of Critical Flaws

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135
Acknowledgments

Available from

NASA Center for Aerospace Information

7115 Standard Drive

Hanover, MD 21076–1320

National Technical Information Service

5301 Shawnee Road

Alexandria, VA 22312

Available electronically at http://www.sti.nasa.gov

Level of Review: This material has been technically reviewed by technical management. 

This work was funded by the NASA Fundamental Aeronautics Program, Supersonics Project.

Trade names and trademarks are used in this report for identifi cation 

only. Their usage does not constitute an offi cial endorsement, 

either expressed or implied, by the National Aeronautics and 

Space Administration.

This work was sponsored by the Fundamental Aeronautics Program 

at the NASA Glenn Research Center.

Erratum

Issued March 2014 for

NASA/TM—2013-217810

Unit-Sphere Anisotropic Multiaxial Stochastic-Strength Model Probability Density Distribution for the Orientation of Critical Flaws

Noel N. Nemeth

September 2013

Page 11: Equation (27) should be

 

 

0 0

2

sin for

22 2 2

T

T

T T

T





       





                                      

 (27) 
NASA/TM—2013-217810 1 

Unit-Sphere Anisotropic Multiaxial Stochastic-Strength 

Model Probability Density Distribution for the 

Orientation of Critical Flaws 

Noel N. Nemeth 

National Aeronautics and Space Administration 

Glenn Research Center 

Cleveland, Ohio 44135 

Summary 

Models that predict the failure probability of monolithic glass and ceramic components under 

multiaxial loading have been developed by authors such as Batdorf, Evans, and Matsuo. These “unitsphere” failure models assume that the strength-controlling flaws are randomly oriented, noninteracting 

planar microcracks of specified geometry but of variable size. This report develops a formulation to 

describe the probability density distribution of the orientation of critical-strength-controlling flaws that 

results from an applied load. This distribution is a function of the multiaxial stress state, the shear 

sensitivity of the flaws, the Weibull modulus, and the strength anisotropy. Examples are provided 

showing the predicted response on the unit sphere for various stress states for isotropic and transversely 

isotropic (anisotropic) materials—including the most probable orientation of critical flaws for offset 

uniaxial loads with strength anisotropy. The author anticipates that this information could be used to 

determine anisotropic stiffness degradation or anisotropic damage evolution for individual brittle (or 

quasi-brittle) composite-material constituents within finite-element or micromechanics-based software. 

1.0 Introduction 

The term “unit sphere” as used herein refers to the models that were developed by Batdorf and Crose 

(1974), Batdorf and Heinisch (1978), Evans (1978), and Matsuo (1981) to predict the probability of 

failure of brittle materials under multiaxial loading. These models use a unit radius sphere representing 

the random orientation of flaws to calculate the effect of multiaxial stresses on material reliability. This 

approach assumes that the strength-controlling flaws are randomly oriented, noninteracting planar microcracks of specified geometry but of variable size. Fracture mechanics relationships for mixed-mode 

(modes I, II, and III) crack growth—combined with the weakest-link theory and integration over the 

surface area of a unit radius sphere representing all possible orientations of microcracks—are used to 

calculate the material probability of failure. 

This unit-sphere methodology was originally introduced within a statistical theory of brittle material 

strength by Weibull (1939), though without consideration for the mechanics of crack growth. Nemeth 

(2013) extended this methodology to predict anisotropic (transversely isotropic) brittle material strength 

response and demonstrated how this could be used as a failure criterion inside finite-element (or alternatively, micromechanics-based) software codes for individual brittle (or quasi-brittle) compositematerial constituents to predict the overall strength response. 

The unique feature of the unit-sphere methodology is the assumption that the presence of flaws in a 

brittle material drives the failure response and that these flaws are microcracks with an assumed orientation and geometry. In this manner, the stochastic strength response of the material can be predicted for 

an arbitrary multiaxial stress state. An interesting consequence of this methodology is that the orientation 

of the critical flaws (the failure-initiating flaws) relative to the applied multiaxial stress state can be 

predicted as a probabilistic distribution. 

This report develops the generalized formulation to describe the probability density distribution 

(PDD) of the orientation of critical-strength-controlling flaws for the unit-sphere model. The derivation 

starts from the most basic level with the description of the weakest link mechanism and the Weibull 
NASA/TM—2013-217810 2 

distribution. Following that is the description of the Batdorf unit-sphere methodology for the isotropic 

strength response. The extensions for the unit-sphere methodology with the anisotropic (transversely 

isotropic) strength response are provided to complete the description of all relevant parameters. These 

extensions are for (1) flaw orientation anisotropy, where a preexisting microcrack has a higher likelihood 

of being oriented in one direction than another direction, and (2) critical strength or fracture toughness 

anisotropy, where the level of critical strength Ic or fracture toughness KIc for mode I crack propagation 

changes with the orientation to the microstructure (for a tensile mode of failure only). This leads to the 

generalized formulation for the PDD of the critical flaw orientation over the unit sphere. Examples are 

provided to demonstrate this methodology for isotropic and anisotropic materials for various stress states 

and off-axis loading. 

If one knows the orientation distribution of critical flaws (and the subsequent direction of crack propagation), the anisotropic stiffness degradation (the anisotropic elastic constants associated with the damaged 

material) can be determined. The author anticipates that this will be helpful for a follow-on phase of this 

effort, not described here, of enabling the unit-sphere failure criterion methodology to work with NASA’s 

micromechanics analysis code/generalized method of cells (MAC/GMC) (Bednarcyk and Arnold (2002)), 

which is based on the GMC family of micromechanics theories, including doubly and triply periodic versions 

of the GMC (Aboudi (1995)) and the High-Fidelity Generalized Method of Cells (HFGMC) (Aboudi et al. 

(2003)). This incorporation will allow the full exercise of the unit-sphere methodology, including incremental 

time/load steps and fatigue analysis (as described in Nemeth et al. (2005)), to predict the durability (strength 

and lifetime) of composite laminates and woven composite structures. 

The unit-sphere methodology provides an improved mechanistic basis to the problem of predicting 

the strength response of an anisotropic material under multiaxial loading in comparison to polynomial 

interaction equation formulations such as Tsai-Wu, Tsai-Hill, and Hashin, among others. A unique 

physical attribute of this model is the ability to predict the critical-flaw-angle distribution for an arbitrary 

stress state. This distribution is derived and demonstrated in this report. 

2.0 Model Description 

2.1 Weibull Distribution 

Consider a stressed component containing many flaws (microcracks), and assume that failure is 

due to any number of independent and mutually exclusive mechanisms (links). Each link involves an 

infinitesimal probability of failure. Discretize the component into n incremental links. The probability of survival (PsV)i of the i

th link is related to the probability of failure (PfV)i of the ith link by 

(PsV)i = [1 – (PfV)i], and the resultant probability of survival of the whole structure is the product of 

the individual probabilities of survival: 

            













           

n

i

i fV i fV

n

i i fV

n

i

sV i

n

i

PsV P P P P

1 1 1 1

Π Π 1 Π exp exp (1) 

where the subscript V denotes volume-dependent terms. Equation (1) arises from the approximation to the 

exponential series 

... 2! 3! 4!1

2 3 4

     

x x x

e x x (for all real values of x) (2) 

where (PfV)i is assumed to be a small number such that the higher order terms in Equation (2) can be 

neglected. 

Assume the existence of a function V(), referred to as the crack-density function, representing the 

number of flaws per unit volume that have a strength equal to or less than . Under a local tensile stress 

i, the probability of failure of the i

th link, representing the incremental volume Vi, is 

 (PfV )i = [V(i) Vi] (3) 
NASA/TM—2013-217810 3 

where the incremental volume Vi is arbitrarily small such that the value of the expression within the 

brackets is much less than one. If one applies a uniform tensile stress , such that  = i for all 

incremental volumes Vi, then from Equation (1) the resultant probability of survival for material 

volume V, where V is the sum of all Vi, is 

PsV = exp[–V()V] (4) 

Equation (4) can also be derived from the Poisson probability density function. This distribution is 

described by (see Hoel et al. (1971), for example) 









   

  

0 elsewhere

0,1,2, ... !

exp( )

( ) ( ) x

x P X x f x

x

 (5) 

where  is a positive number. The real-valued function (x) is the discrete density function of random 

variable X where P(X = x) is the probability that a discrete real-valued random variable X equals a 

possible value x. The Poisson distribution approximates the binomial distribution for large values of n, 

where n is the number of Bernoulli trials with success probability p= /n at each trial. Equation (4) is 

obtained when P(X = 0) is computed for n = V and p = V(); hence, 

P P X   V sV V exp( ) exp η (σ) 0!

exp( ) ( 0)

0

         (6) 

Equation (6) calculates the probability that no flaws of strength  or less are present in the material 

volume V and, therefore, represents the survival probability of the material under applied load . 

The probability of failure for the uniformly stressed volume V is 

PfV V     1 exp   V    (7) 

where V is the total volume. 

Weibull introduced a three-parameter power function for the crack-density function V(): 

  mV

oV

uV

V 

















  

   (8) 

where uV is the threshold stress parameter. The scale parameter oV then corresponds to the stress level 

where 63.21 percent of tensile specimens with unit volumes would fracture. The scale parameter has 

dimensions of stress  (volume) 1/ , mV where mV is the shape parameter (Weibull modulus)—a dimensionless parameter that measures the degree of strength variability. As mV increases, the dispersion decreases. 

The threshold stress parameter uV is usually taken as zero for brittle materials (ceramics and glasses). 

This parameter is the value of the applied stress below which the failure probability is zero. When this 

parameter is zero, the two-parameter Weibull model is obtained in Equation (7). The two-parameter 

crack-density function is expressed as 



   















  V

V

m wV

oV

m

Vη ( ) k (9) 

where kwV = (oV)

–mV is the uniaxial Weibull crack-density coefficient (for an applied uniaxial load). 
NASA/TM—2013-217810 4 

2.2 Batdorf Unit-Sphere Model for Isotropic Strength Response 

Batdorf and Crose (1974) propose a statistical theory in which attention is focused on cracks and their 

failure under stress. Flaws (assumed to be microcracks in shape and size) are taken to be uniformly distributed and randomly oriented in the material bulk. Fracture is assumed to depend only on the tensile stress 

acting normal to the crack plane; hence, shear insensitivity is inherent to this first model. Subsequently, 

Batdorf and Heinisch (1978) included the detrimental effects of shear traction on a flaw plane. Their 

method applies fracture mechanics concepts by combining the crack geometry and a mixed-mode fracture 

criterion to describe the condition for crack growth. 

In the Batdorf unit-sphere theory, the incremental failure probability PfV under an applied multiaxial 

state of stress  at a given location in the component can be described as the product of two probabilities: 

  PfV Ieqc V P1V P2V  , ,   (10) 

where P1V is the probability of the existence in incremental volume V of a crack having an equivalent, 

or effective, critical stress between Ieqc and Ieqc + Ieqc. Critical stress Ic is defined as the remote, 

uniaxial fracture strength of a material containing a given crack in mode I loading. The term Ieqc denotes 

an effective (or equivalent) critical mode I stress from applied multiaxial stresses. The second probability, 

P2V, denotes the probability that a crack of critical stressIeqc will be oriented in a direction such that an 

effective stress Ieq (which is a function of fracture criterion, stress state, and crack configuration) satisfies 

the condition Ieq  Ieqc. The effective stress Ieq is defined as the equivalent mode I stress that a flaw 

would experience when subjected to a multiaxial stress state that results in mode I, II, and III crack surface displacements, and Ieqc is the threshold value of Ieq where unstable catastrophic crack growth 

ensues. An incremental volume V is used in Equation (10) because an infinitesimal volume dV cannot 

enclose a crack of critical stressIeqc and associated critical crack length ac. 

The effective stress Ieq represents an equivalent normal stress on the crack face from the combined 

action of the normal stress n and the shear stress . The microcrack orientation is defined by the angular 

coordinates  and , where the direction normal to the plane of the microcrack is specified by the radial 

line defined by  and  in space (see Fig. 1(a)). For the sake of brevity, the development of the effective 

stress equations is not shown (for details, see Nemeth et al. (2003, 2005)). For a tensile mode of failure 

assuming a penny-shaped crack with the Shetty mixed-mode fracture criterion (Shetty (1987)), the 

effective stress becomes 

 

2

2

I

1 4 2 2 eq n n

C

                      

 (11) 

where  is Poisson’s ratio and C is the Shetty shear-sensitivity coefficient, with values typically in the 

range 0.80 ≤ C ≤ 2.0. As C increases, the response becomes progressively more shear insensitive. Shear 

increases the equivalent stress as shown in Equation (11), and this has a deleterious effect on the predicted 

material strength. For a penny-shaped crack with a material having a Poisson’s ratio  ~ 0.22 and C = 0.80, 

0.85, 1.05, and 1.10, Equation (11) approximates, respectively, the following criteria: Ichikawa’s maximum 

energy-release-rate approximation (Ichikawa (1991)), the maximum tangential stress (Erdogan and Sih 

(1963)), the maximum strain-energy-release-rate formulation of Hellen and Blackburn (1975), and colinear 

crack extension. The value of C also can be fit empirically to experimental data—either on introduced 

cracks (as is done in Shetty (1987)) or on specimens being tested multiaxially. 
NASA/TM—2013-217810 5 

Figure 1.—Projection of equivalent stress onto a unit radius sphere. 

(a) Cauchy stress components on an infinitesimal tetrahedron 

resolving a normal stress, n, and a resultant shear stress,  on a 

plane normal to the direction defined by angular coordinates  and . 

(b) Projection of equivalent stress onto a unit radius sphere in the 

global coordinate system. The unit radius sphere represents all 

possible flaw orientations, where Ieq is an equivalent (or effective) 

stress; x, y, and z are normal orthogonal stress components; and 

xy, yz, and zx are shear stress components. An infinitesimal area, dA, 

on the unit sphere represents a particular flaw orientation (a direction 

normal to the flaw plane), and Ieq is a function of n and  for an 

assumed crack shape and multiaxial fracture criterion. 

 
NASA/TM—2013-217810 6 

A material is assumed to have n individual failure modes. Material reliability (probability of survival) 

is assumed to be the product of the survival probability of all of the failure modes. These failure modes 

can be from different flaw populations, from the flaws associated with different material constituents 

(i.e., composites as described in Nemeth (2013)), or from different physical mechanisms that govern the 

material failure such as tensile loading versus compressive loading. Here the compressive mode of failure 

is assumed to be of a different nature than the tensile failure mode—possibly involving the interaction of 

arrested cracks (flaws with initial crack growth that subsequently stops because of the angle of the crack 

growth and the interaction with the local stress field). Regardless, the Weibull distribution here is 

assumed to describe the stochastic strength response phenomenologically. This is argued further in 

Nemeth and Bratton (2011). The compressive failure mode is assumed to be controlled by shear stress, 

and a simple Tresca-like effective stress Ieq relation can be prescribed as 

Ieq = 2 (12) 

The multiplier of 2 in Equation (12) was chosen so that the maximum effective stress on the unit sphere in 

pure uniaxial compression is equal to the applied compressive stress. When the normal stress component n

on the crack face is tensile, then the value of Ieq in Equation (12) is set to zero in the unit sphere evaluation. A more rigorous treatment to account for compression would include the frictional effects of the 

opposed crack surfaces in contact in the effective stress relation. However, Equation (12) is sufficient for 

the demonstration purposes of this report. 

2.2.1 Relationship of Crack Probability of OccurrenceP1V to the Crack-Density Function V 

The strength of a component containing a flaw population is related to the critical flaw size, which is 

used implicitly in statistical fracture theories. Batdorf and Crose (1974) describe the crack probability of 

occurrence P1V as 

 

eqc

eqc

V eqc PV V I

I

I

1 d

d

d





 

   (13) 

This is better understood using Figure 2. In Figure 2(a), for example, a point on the curve V(Ieqc)V for 

the value Ieqc,1 represents the number of flaws in V having a strength equal to or less than Ieqc,1. Then, 

from Figure 2(a), 

V(Ieqc,2)V – V(Ieqc,1)V 

 = Probability that a flaw of strength Ieqc,2 or less is in V 

 – Probability that a flaw of strength Ieqc,1 or less is in V (14)

 = Probability that a flaw of strength between Ieqc,1 and Ieqc,2 is in V

Equivalently in Figure 2(b), 

 [V(Ieqc) + dV]V – V(Ieqc)V 

 = Probability that a flaw of strength (Ieqc + dIeqc) or less is in V

 – Probability that a flaw of strength Ieqc or less is in V (15) 

 = Probability that a flaw of strength between Ieqc and (Ieqc + dIeqc) is in V

 = V dV

The slope s of the relation V(Ieqc)V is s = (dV /dIeqc)V as shown in Figure 2(b); therefore, 

multiplying the slope s by dIeqc yields the result of Equation (13). 
NASA/TM—2013-217810 7 

Figure 2.—Crack-density function, V, versus the critical equivalent 

stress, Ieqc, where V is incremental volume, shown to explain 

Equation (13). (a) Incremental differences in values. (b) The same 

differences described in terms of infinitesimal values. 

2.2.2 Random Flaw Orientation as Described by a Unit Sphere 

The term P2V shown in Equation (10) is expressed as 

 



    4, I

2

eqc P V (16) 

where (, Ieqc) is the area of the solid angle projected onto the unit radius sphere in stress space (see 

Fig. 1(b)) containing all the crack orientations for which Ieq Ieqc for the applied far-field multiaxial 

stress state . The infinitesimal area, dA, on the unit sphere represents a particular flaw orientation 

(a direction normal to the flaw plane), and Ieq is an equivalent stress, which is a function of an assumed 

crack shape and multiaxial fracture criterion. The constant 4 is the surface area of a unit radius sphere 

and corresponds to a solid angle containing all possible flaw orientations. Equation (16) assumes that the 

flaws are randomly oriented with a uniform distribution (equal likelihood of being oriented at any angles 

 and ). If it is assumed that the other properties like crack geometry and KIc are also independent of 

 and  the strength response becomes isotropic (independent of the applied direction of load). 

2.2.3 Unit-Sphere Failure Probability Formulation for Isotropy 

Integrating Equation (10) with respect to Ieqc and using Equation (4) yields the probability of 

survival of a (uniformly stressed) volume element Vi: 
NASA/TM—2013-217810 8 

     

i

eqc

eqc

eqc V eqc

sV i

eq P V





































 



      

I ,max

0 I

I

I I d

d

d

4

, exp (17) 

where Ieq,max is the maximum effective stress that a randomly oriented flaw could experience from the 

given stress state. Hence, integrating over the entire volume of the component results in the failure 

probability: 

   





































 



        P V

V eqc

eqc

eqc V eqc

fV eq d d

d

d

4

, 1 exp

I ,max

0 I

I

I I (18) 

where Ieq,max is the maximum effective stress that a randomly oriented flaw could experience in the 

component and  is allowed to vary within the component. 

Because of the effect of multiaxial stress states, the crack-density function V(Ieqc) is altered in the 

Batdorf formulation. Batdorf and Heinisch, (1978) have typically used a power-law formulation for 

V(Ieqc), which leads to the Batdorf crack-density function of the form 

  V

V

V

m

oV

m

m BV eqc

V eqc BV eqc

k

k



  

       (x, y,z, , ) (x, y,z, , ) I

I I (19) 

where x, y, and z correspond to the location within the component and  and  are the unit-sphere 

orientation angles. The term kBV is the Batdorf crack-density coefficient, and kBV is the normalized 

Batdorf crack-density coefficient whereby kBV  kBV kwV as described in Nemeth et al. (2003, 2005). 

The material Batdorf crack-density coefficient kBV and the Weibull modulus mV are evaluated from 

experimental inert strength fracture data. By convention, kBV is calibrated to a uniaxial stress state. 

Batdorf and Crose (1974) initially proposed a Taylor series expansion for V(Ieqc), but this method has 

computational difficulties. Note that V(Ieqc) has units of inverse volume. 

Although the Weibull (Eq. (9)) and Batdorf (Eq. (19)) crack-density functions are similar in form, 

they are not the same. The Weibull function simply depends on the applied uniaxial stress distribution 

and is the only term other than mV and the volume necessary to calculate PfV. It does not extrapolate to 

other multiaxial stress states. The Batdorf function depends on the equivalent mode I strength of the 

crack Ieqc, which is probabilistic and must be integrated over a range of values for a given stress state. 

Furthermore, to obtain PfV, a crack orientation function, P2V, must be considered in addition to the density 

function and the volume. Finally, the Batdorf coefficient kBV cannot be calculated from inert strength data 

until a fracture criterion and crack shape are chosen—in contrast to the Weibull coefficient kwV, which 

depends only on normalizing or calibrating to experimental rupture data. 

To determine a component probability of failure, one must evaluate P2V (Eq. (16)) for each elemental 

volume Vi, within which a uniform multiaxial stress state  is assumed. The solid angle (, Ieqc) 

depends on the selected fracture criterion, the crack configuration, and the applied stress state. For 

multiaxial stress states, with few exceptions, (, Ieqc) must be determined numerically. For a sphere of 

unit radius (see Fig. 1(b)), an elemental surface area of the sphere is dA = sin  d d. If we project onto 

the spherical surface the equivalent (effective) stress Ieq(, , ), the solid angle (, Ieqc) will be the 

area of the sphere containing all the projected equivalent stresses satisfying Ieq Ieqc. Note that the 

symmetry of Ieq means that only one-eighth of the unit sphere needs to be integrated for an isotropic 

strength response. However, for transverse isotropy, one-half of the sphere should be integrated. When 

one-half of the unit sphere is considered, 
NASA/TM—2013-217810 9 

              

 





  

   2

0

2

0 I I

I , sin d d

2

1 4

,

eq eqc

eqc H (20) 

where 

H(Ieq, Ieqc) = 1 Ieq  Ieqc

H(Ieq, Ieqc) = 0 Ieq < Ieqc

The Heaviside function H(Ieq, Ieqc ) is a function of fracture criterion, stress state, crack configuration, and crack orientation. Substituting into Equation (18) and integrating first with respect to Ieqc for 

each individual angle or  on the unit sphere changes the component failure probability to (Batdorf 

(1978a,b)) 

  



 





    

       

V

PfV V eq sin d d dV

2

1 1 exp

2

0

2

0 I (21) 

where 

      x, y,z,, I I

mV

V eq BV eq k

For a given incremental volume, Ieq(x, y, z, , ) is the projected equivalent stress over the unit radius 

sphere in coordinate stress space as shown in Figure 1. 

Equation (21) circumvents the involved numerical integration of (, Ieqc) as was originally 

developed in Nemeth et al. (1990). Equations (18) and (21) are equivalent formulations; however, 

Equation (21) is more convenient for computational purposes with the few exceptions where specialcase closed-form solutions exist for particular stress states, crack geometries, and fracture criteria. 

2.3 Unit-Sphere Model for Transversely Isotropic Strength Response 

Two different physical mechanisms were considered in order to extend the unit-sphere model to 

account for anisotropic strength response: (1) flaw orientation anisotropy, whereby a preexisting microcrack has a higher likelihood of being oriented in one direction than another direction and (2) critical 

strength or fracture toughness anisotropy, whereby the level of critical strength Ic or fracture toughness 

KIc for mode I crack propagation changes with regard to the orientation of the microstructure (for a tensile 

mode of failure only). Flaw orientation anisotropy was previously considered by Buch et al. (1977), and 

critical strength anisotropy was previously considered by Batdorf (1973). Both models were developed to 

simulate the graphite anisotropic strength response. Extensions to these models, described and demonstrated in Nemeth (2013) for transversely isotropic strength response, are reproduced herein because 

terms from these formulations are used in the flaw orientation probability density function (PDF) 

formulation. The extensions include shear sensitivity for flaws and an improved functional form for the 

anisotropy equations. 

2.3.1 Flaw Orientation Anisotropy 

Flaw orientation anisotropy refers to the situation where a flaw has a higher likelihood of being 

oriented in one direction than another for a given critical strength. This means that a material will be 

stronger on average in one direction than another. An isotropic brittle material is equally strong in any 

direction, and thus its flaws are uniformly randomly oriented. However, in components made by 

processes such as extrusion or hot pressing, which induce texture, a bias will exist in the orientation 

distribution of processing flaws. Also, components finished by surface grinding will contain machining 
NASA/TM—2013-217810 10 

damage in the form of surface cracks that are oriented parallel and transverse to the grinding direction. 

In composite materials, the interface or interfacial layer between the fiber and the matrix can act as a flaw 

with an orientation bias that induces an anisotropic strength response. 

In the following unit-sphere formulation, consideration is given to a flaw that has a higher likelihood 

of being oriented in one direction than another. If an argument analogous to Equation (10) is used, the 

incremental probability of failure under the applied multiaxial state of stress  for a given  and  orientation of the (direction or vector) normal to the particular flaw plane can be written as the product of two 

probabilities: 

PfV(, , , Ieqc, V) = P1V (, ) d dH(Ieq, eqc) (22) 

where P1V is the probability that a crack having an equivalent critical stress between Ieqc and 

Ieqc + dIeqc exists in an incremental volume element V. The second probability, (, ) d d, 

denotes the probability that a crack of critical stressIeqc is oriented in the range between  and ( + d) 

and between  and ( + d), where the orientation of the microcrack is described by the vector that is 

normal to the plane of the microcrack. The Heaviside function H(Ieq, Ieqc ) is identical to that as defined in 

Equation (20), where H(Ieq, Ieqc) = 1 when Ieq(, ,) Ieqc and where H(Ieq, Ieqc ) = 0 when 

Ieq(, ,)  Ieqc. The effective stress Ieq(, ,) is defined as the equivalent mode I stress oriented at 

coordinates  and that a flaw would experience when subjected to a multiaxial stress state Σ that results 

in mode I, II, and III crack surface displacements, andIeqc is the threshold value of Ieq where unstable 

catastrophic crack growth ensues. 

Analogous to Equations (16) and (20) for the unit sphere, 

                2

0

2

0

P2V , H Ieq, Ic sin d d (23) 

where 

           

       2

0

2

0 , sin d d

( , ) ( , ) (24) 

and 

H(Ieq, Ic) = 1 Ieq  Ic

H(Ieq, Ic) = 0 Ieq < Ic

Equation (24) is explained further in Section 2.4. The function () describes the degree of anisotropy 

of the flaw orientation, where the normal direction to the flaw plane is given by angles  and . 

Equations (23) and (24) modify Equation (21) for probability of failure: 

    



 



               

V

PfV 1 exp , V eq sin d d dV

2

0

2

0 I (25) 

where 
NASA/TM—2013-217810 11 

      V V

V

m

oV

m

eq

BV m V eq BV eq k k 





















  

       x, y,z, , x, y,z, , I

I I

The normalized crack-density coefficient kBV now includes the effect of flaw orientation anisotropy and 

is evaluated numerically from Equation (25) for a uniaxial stress state applied in one of the material 

coordinate system axis directions. 

For a transversely isotropic strength response,  in Equation (24) is only a function of . Buch et al. 

(1977) introduced a cosine power function for () = [cos ()] where  is a constant. This relation was 

modified in Nemeth (2013) to enhance the functional flexibility and is expressed as 

 

 

cos 0

2

0

2

L

L

L

L

                      



       





 (26) 

Alternatively, 

 

 

0 0

2

sin for

22 2 2

T

T

T T

L





       





                                          

 (27) 

where  and  are constants (with subscripts L or T) that control the degree of anisotropy, and 0    /2 

and   0. When  = /2 and  = 0, the isotropic strength response is obtained. Equations (26) and (27) 

are defined for one-half of the unit sphere (the top half, as shown in Fig. 1(b) where 0    /2). 

Referring to Figure 3, Equation (26) represents the “polar-cap” or longitudinal L distribution of flaws and 

Equation (27) represents an “equatorial-belt” or transverse T distribution of flaws. The polar-cap distribution describes crack planes that are distributed symmetrically (centered) about a plane (in this case, the 

yz plane), and the equatorial-belt distribution describes crack planes that are distributed symmetrically 

(centered) along a line (in this case the x axis). For the polar-cap distribution, kBV is evaluated for a 

uniaxial stress along the x direction; and for the equatorial-belt distribution, it is evaluated in the y

(or optionally the z) direction. 

The separate polar-cap and equatorial-belt distributions are introduced to describe individual and 

distinctly different failure modes. For a unidirectional fiber-reinforced composite, the polar cap can be 

used to represent the fiber strength distribution and the equatorial belt can be used to represent the matrixfiber interface. The equatorial-belt and the polar-cap distributions can be considered to equivalently 

represent global failure planes, which are referred to as “action planes” in the Puck multiaxial strength 

model for composites (see Lutz (2006) for a description). The angular width of the belt or cap with regard 
NASA/TM—2013-217810 12 

Figure 3.—Unit sphere with probability distribution functions (PDFs) describing the anisotropy of flaw 

orientation, where  and are angular coordinates and L and T are constants in the flaw-orientation 

anisotropy function representing longitudinal and transverse distributions. Orientation is described with the 

normal to the crack plane. In this figure, two orientation functions are described: (1) a polar-cap distribution 

describing crack planes distributed symmetrically (centered) about a plane (in this case the y–z plane) 

and (2) an equatorial-belt distribution where crack planes are distributed symmetrically (centered) along a 

line (in this case, the x axis). 
NASA/TM—2013-217810 13 

to angle  indicates the maximum extent of scatter of a flaw orientation. Flaws or action planes oriented 

outside of the scatter bands of the cap and belt are assumed to not exist or not contribute to the likelihood 

of failure. An example of this methodology is provided in Nemeth (2013) for a unidirectional polymer matrix 

composite (PMC) under biaxial loading. 

2.3.2 Critical Strength or Stress Intensity Factor Anisotropy 

Batdorf (1973) approached strength anisotropy using the Ic strength ellipsoid approach (describing 

an ellipsoid rather than a unit sphere). This was applied to graphite—a mildly anisotropic material. In this 

report, strength anisotropy is used with the critical stress-intensity factor (fracture toughness) KIc varying 

with orientation angle on the unit sphere. This is functionally equivalent to Ic varying with orientation, 

and it is also functionally equivalent to the size of the flaw changing with the orientation angle. With 

regard to a ceramic matrix composite (CMC), where failure from loading in the fiber direction is by 

matrix cracking with large-scale fiber bridging, fracture toughness cannot be defined on the global scale 

of the structure. In that case, one has to consider critical strength as a metric. However, it is acceptable to 

use fracture toughness on the local scale at the crack tip, where micromechanics can account for the 

bridging explicitly. In this report, fiber bridging is not considered directly. Of first-order importance 

herein is that local fracture toughness could change with the orientation of the flaw plane (or action plane) 

and the applied loading. The specific micromechanics of how this might occur is not considered here. 

The modeling approach taken here is similar to that described previously for flaw-orientation anisotropy (Section 2.3.1). The critical strength Ic is defined as the fracture strength of the crack in mode I 

loading and is proportional to KIc. Therefore, for anisotropic KIc(, ) = (c Ic(, )) =   c f   I ,max I c c  ,   , 

where c is a constant ( c  Y ac with crack-shape geometry factor Y and critical crack length ac), Ic,max

is the maximum value of Ic over the unit sphere (for all  and ), and fIc , is a normalized function 

expressing the degree of this anisotropy. The unit-sphere formulation described previously is used, except 

the Heaviside step function in Equations (20) and (23) is modified as follows: 

   

    

















 

 



  

 

 



  

I ,max I

I

I I

I ,max I

I

I I

, , 0 

, , 1 

eqc c

eq

eq eqc

eqc c

eq

eq eqc

f

H

f

H

(28)

where Ieqc is substituted for Ic to indicate generalized mixed-mode loading. 

For a transversely isotropic response where anisotropy is only a function of angle , the function 

fIc    over the top half of the unit sphere (0    /2) is arbitrarily defined for the L distribution as 

 























     

    













  

















  



 





  



2

1

0 1cos 1

2

1

2

1 ( ) 1

I

I

c L

L

L L c

f

r

f

L

 (29) 

and for the T distribution 
NASA/TM—2013-217810 14 

 























      

















  

















 



 



 



   

  

  



 





  

      



2 2

1 1

2

cos

2

1

2

1 ( ) 1

2

1 0

I

I

T

T

T

T

c

c T

r

f

f

T

 (30) 

where L, T, L, T, rL,, and rT are constants. The L subscript relates to the polar-cap strength anisotropy 

distribution for crack planes symmetrically distributed (centered) about a plane (in this case, the yz

plane), and the T subscript relates to the equatorial-belt strength anisotropy distribution for crack planes 

symmetrically distributed (centered) along a line (in this case the x axis) (see also Fig. 3 for a reference 

frame). Figure 4 shows a schematic of Equations (29) and (30) in stress-strength space. 

Figure 4.—Schematic of Equations (29) and (30) in normalized 

strength space. Equation (29) correlates to the polar-cap 

distribution, and Equation (30) correlates to the equatorial-belt 

distribution, as shown in Figure 3. Here, L, T, L, and T are 

constants in the critical mode I stress intensity anisotropy 

function representing longitudinal and transverse distributions; 

f

Ic    is the normalized anisotropy function of the critical 

mode I stress-intensity factor, KIc, or critical mode I strength, 

Ic, as a function of angle  and I ,x I ,y and I ,z , , c c c are the 

normalized (by Ic,max) orthogonal critical strength components. 

 
NASA/TM—2013-217810 15 

2.3.3 Generalized Unit-Sphere Failure Probability Formulation for Transverse Isotropy 

Equations (26), (27), (29), and (30) are arbitrary functions chosen for their flexibility to fit to data. 

Other functions could be easily substituted. For an isotropic distribution of the orientation of the flaws, 

(, ) is 1/2, which is consistent with Equation (21). Although not detailed herein, the Weibull 

modulus can also be made anisotropic, varying with angles  and by using, for example, the 

normalizing Equations (29) and (30) to modify m (dividing m by c fI ()). Therefore, in its most 

generalized form for transverse isotropy (considering flaw and critical strength orientation anisotropy), 

failure probability can be expressed as 

 

















      



         

V eqc

eqc

V eqc PfV V eq ( , ) sin d d d d

d

d

1 exp 2

0

2

0 I 0 I

I ,max I (31) 

and, if a power law is assumed for the crack-density function, 

   

 

 

 

( , )

I oV

I

( , )

I

I

I , ,

x, y,z, ,

,

x, y,z, , ( , )

   





















    

    



















 

  

    

V mV

c

eq

BV

m

c

eq

V eq BV

f

k

f

k (32) 

where critical strength or KIc orientation anisotropy is accounted for with the addition of the term c fI . The 

Weibull modulus mV is allowed to be a function of  and . Note that the Weibull scale parameter oV is 

used in the crack-density formulation and that it is also allowed to be a function of  and , but that it is 

assumed to be constant with spatial location x, y, z. Equation (31) can then be more simply expressed as 

  



















      



















    

         

 

V

m

c oV

eq

fV BV V

f

x y z

P k

V

( , ) sin d d d , ( , )

( , , , , ) 1 exp 2

0

2

0

( , )

I

I (33) 

2.4 Orientation Distribution of Critical Flaws Under Multiaxial Load 

The unique feature of the unit-sphere methodology is the assumption that the flaws that are inherently 

present in a brittle material will drive the failure response and that these flaws are microcracks with an 

assumed orientation and geometry. In this manner, the stochastic strength response of the material can be 

predicted for an arbitrary multiaxial stress state. An interesting consequence of this methodology is that the 

orientation of the critical flaws (the failure-initiating flaws) relative to the applied stresses can be predicted 

as a probabilistic distribution. The PDF of the orientation of the critical fracture-causing flaws can be 

obtained from Equations (10), (13), (16), and (19) for an isotropic strength material; from Equations (22), 

(23), and (24) for flaw orientation anisotropy; and from Equation (28) for critical strength (fracturetoughness) orientation anisotropy. How this is obtained is described subsequently for transverse isotropy. 

For a given flaw (a flaw that is assumed to exist at some orientation  and ) and considering all 

possible flaw orientations, the likelihood must be 1 (in other words, 100-percent certainty) that the flaw 

exists at some orientation α and β. In other words, 

 Probability that a flaw exists at some  and  



 



            



2

0

2

0

1 c ( , )sin d d (34) 

where c is a normalizing constant and the function () describes the degree of anisotropy of flaw 

orientation at angular coordinates  and  (where the normal direction to the flaw plane is given by angles 
NASA/TM—2013-217810 16 

and ). The function  was defined in Equations (26) and (27) for a transversely isotropic strength response 

(which is expressed only as a function of angle ). Therefore, in Equation (34) c  must equal 

          

 2

0

2

0 ( , )sin d d

1

c (35) 

The probability that the flaw is oriented in the range between  and ( + d) and between  and 

( + d) is given by (, ) ddin Equation (24). In this case, an infinitesimal area element dA = d

d on the unit sphere is represented as shown in Figure 5. In contrast to Figure 1(b), the sine term is not 

used here to describe dA. Figure 1(b) maintains a convention previously used in other references such as 

Nemeth et al. (2003, 2005). The sine term is only necessary when integration over the whole unit sphere 

is required. In order to determine the PDF of the orientation of the critical flaws, dA must be kept constant 

versus angular coordinates  and . If dA = sin  d dwere used, the nonsensical condition would be 

obtained that the probability of a flaw being oriented at would always be zero regardless of the 

stress state and the magnitude of the stress. 

Equation (24) is the PDF of the degree of anisotropy of the flaw orientation. Integration of the 

numerator of Equation (24) over the entire unit sphere for dA = sin  d dyields a value identical to the 

denominator, satisfying the condition for a PDF that integration over all possible values must equal 1 (for 

the cumulative distribution function). 

Figure 5.—Infinitesimal area element, dA, on the unit sphere. 

In this case, dA = d d for a given  and  orientation of 

the normal to the flaw plane, where Ieq is the effective 

stress. 

The probability that a flaw of strength Ieqc or less that is within the incremental volume element V, 

where Ieq(, , ) Ieqc, and that is oriented between  and ( + d) and between  and ( + d), will 

fail is (see also Eq. (3)) 
NASA/TM—2013-217810 17 

 

         

                       2

0

2

0 ( , )sin d d

( , ) d d ( , ) ( , , , ) ( , ) d d ( , ) V P V V V

fV V (36) 

where PfV(, , , V) is expressed for a constant infinitesimal area element dA = d d on the unit 

sphere as shown in Figure 5 regardless of the specific values of angles  and . 

The failure probability for the whole unit sphere for incremental volume element V is 

 

   

 

     





 



         

   2

0

2

0

2

0

2

0

( , )sin d d

( , ) ( , )sin d d

( , )

V

P V

V

fV (37) 

Reexpressing Equation (37) using the crack-density function of Equation (32), using Equation (35) for 

c , and considering critical strength or KIc orientation anisotropy yields 

V

f

P V c k

mV

c oV

eq

fV BV 





















   



















    

  

         

 



2

0

2

0

( , )

I

I sin d d

( , ) ( , )

( , ) ( , ) ( , ) (38) 

The probability that the critical flaw (the failure-initiating flaw) is oriented in the range between 

and ( + d) and between  and ( + d) that is under a given stress state  is found by dividing 

PfV(, , , V) in Equation (36) by PfV(, V) in Equation (37): 

  

            

             2

0

2

0 ( , ) ( , )sin d d

( , ) d d ( , ) ( , , )

V

V c (39) 

The numerator in Equation (39) refers to an infinitesimal area element dA = dα dβ on the unit sphere as 

shown in Figure 5. The PDF of Equation (39) is therefore 

 

             

           2

0

2

0 ( , ) ( , )sin d d

( , ) ( , ) ( , , )

V

V (40) 

Integration of the numerator of Equation (40) over the entire unit sphere for dA = sin  d dyields 

a value identical to the denominator, satisfying the condition for a PDF that integration over all possible 

values must equal 1 (for the cumulative distribution function). Substituting Equation (32) into 

Equation (40) yields 

   

 

 



   















 

     

















 

     

    

2

0

2

0

( , )

I

I

( , )

I

I

sin d d

( , )

( , ) ( , )

( , )

( , ) ( , )

( , , ) V

V

m

c

eq

BV

m

c

eq

BV

f

k

f

k

 (41) 
NASA/TM—2013-217810 18 

However, for computational convenience, the effective stress Ieq should be normalized by the scale 

parameter oV and the normalized Batdorf crack-density coefficient kBV . If it is assumed that oV and mV

are independent of angles  and  

   



   















  

     

















  

     

    

2

0

2

0

I

I

I

I

sin d d

( , )

( , ) ( , )

( , )

( , ) ( , )

( , , ) V

V

m

c oV

eq

m

c oV

eq

f

f

 (42) 

Equation (42) describes the PDF of the orientation of critical flaw normals from the applied stress 

state . However, these equations only describe the probability of the orientation of the critical flaw, they 

do not describe the angle of subsequent crack propagation. For noncoplanar crack growth, the action of a 

simultaneous normal stress and a shear stress will cause the crack to propagate at an angle  relative to the 

plane of the initial flaw or microcrack. The anisotropy of material properties can also affect the angle of 

propagation. This direction of growth will be the direction of cracking that would be observed in rupture 

experiments (or loading experiments to the point of reaching the proportional limit or some level of 

damage threshold). Actually observing the orientation of the critical flaw would be difficult. This report 

does not consider the topic of crack propagation in an anisotropic medium. Nonetheless, the methodology 

described in this report provides a mechanistic-based framework for which to describe orientation 

distribution for critical flaws. 

3.0 Examples 

3.1 Isotropic Material Critical Flaw Orientation Probability Density Distribution 

Fortran was used to implement a gaussian quadrature numerical integration algorithm to compute the 

generalized probability of failure equation (Eq. (33)) and the PDF equation for the orientation of critical 

flaws (Eq. (42)). This algorithm procedure is the same as that described in Nemeth et al. (2005), except 

that the integration is performed sequentially over smaller discrete sectors of the unit sphere. The numerical output from the Fortran code was used with MATLAB (MathWorks) software to plot the PDD (which 

herein are the results from the PDF) over the unit sphere. Values for parameters related to critical strength 

anisotropy (Eqs. (29) and (30)) are sensitive to the number of sampling points (number of gauss points, 

NGP) used in the numerical integration. When a greater number of points are used, the numerical solution 

asymptotically approaches a single value. In the following plots of the unit sphere (Figs. 6 to 11), a 

smaller number of sampling points was used (NGP = 5) for computational efficiency. This somewhat 

affected the numerical values of parameter rT but not the trends seen in the figures. Results for Table I and 

Figures 12 to 15 were generated for NGP = 30, which represents fully converged solutions. 

 
NASA/TM—2013-217810 19 

TABLE I.—VALUES OF PARAMETER rT (SEE EQ. (30)) FOR VARIOUS STRENGTH RATIOS 

OF THE STRONGEST TO WEAKEST MATERIAL STRENGTH DIRECTIONS 

UNDER UNIAXIAL TENSILE LOADING 

Strength ratio Parameter rT

Parameter, T, deg 

10 90 

Shetty shear sensitivity, C

1.1236 100.0 1.1236 100.0 

1.25 1.57 1.47 1.39 1.29 

1.5 1.95 1.83 1.82 1.59 

2.0 2.68 2.52 2.83 2.21 

3.0 4.15 3.91 5.44 3.51 

4.0 5.64 5.32 8.83 4.87 

6.0 8.68 8.18 17.90 7.68 

8.0 11.77 11.10 30.01 10.59 

10.0 14.91 14.06 45.18 13.58 

Figure 6 shows an example of the PDD from Equation (42) of critical flaws on the unit sphere for an 

isotropic material with a 100-MPa uniaxial tensile load applied in various directions. Figures 6(a), (b), 

and (c) show results for a uniaxial load applied in the x, y, and z material axes, respectively. The plots 

show the direction normal to the flaw plane as depicted schematically in Figure 1. As would be expected, 

the “spot” is centered on the direction of applied load, with the highest probability of orientation 

perpendicular to the applied load (dark red (center of “bullseye” shapes) indicates the highest probability 

and blue (background color) indicates the lowest probability). Figures 6(a), (b), and (c) are shown to 

demonstrate that the numerical algorithm and plotting software correctly reproduce the same results 

regardless of the direction of the applied load. Figure 6(d) shows the PDD numerical distribution versus 

angles  and  for an applied uniaxial tensile load in the z direction. Figure 6 is generated for an 

incremental unit volume with a load level corresponding to PfV = 0.6321 (for oV = 100 MPam1/10), 

mV = 10, and C = 1.5 for mild shear sensitivity (using Eq. (11)). All subsequent plots were generated for 

PfV = 0.6321 and oV = 100 MPam1/10. 

The size (angular extent) of the “spot” (or distribution) shown in Figure 6 is a function of mV and C . 

This is demonstrated in Figure 7 for a uniaxial tensile load applied in the y direction. The plots show that 

the higher mV is, the more localized (more confined or highly aligned) the “spot” with the applied load is 

(flaws are more perpendicular to the applied load) as shown in Figures 7(a), (c), and (e). The same is also 

true with C , where for the shear-insensitive condition (C = 100.0), the spot is more tightly distributed as 

shown in Figures 7(b), (d), and (f) and, as C gets smaller (shown as C = 1.5 for mild shear sensitivity), 

the spot is less tightly distributed as shown in Figures 7(a), (c), and (e). 

Figure 8 shows the effect of shear-sensitivity parameter C for an applied shear stress yz for mV = 10. 

Results for increasing shear sensitivity for C = 100.0 (shear-insensitive), C = 1.5 (mildly shearsensitive), and 0.82 (highly shear-sensitive) are shown in Figures 8(a), (c), and (e), respectively. In 

Figures 8(a) and (b), only the normal tensile stress component contributes to failure, and the spot of 

highest probability is centered on  = 90 and  = 45. In Figures 8(c) and (d), for the case of mild shear 

sensitivity, the spot of the highest probability bifurcates and the two spots are symmetrically offset at a 

small  angle from  = 90 and  = 45. This indicates that the effect of shear sensitivity on the microcrack alters the angle of most probable fracture for an applied shear stress. This is further shown in 

Figures 8(e) and (f) where a high shear sensitivity shows two distinctly different fracture planes, but their 

symmetry is oriented about  = 90 and  = 45. An additional note about Figures 8(e) and (f) is that, 

when the normal stress on the flaw plane becomes compressive, the failure probability associated with 

that orientation for the tensile failure mode is set to zero in the numerical method; hence, there is some 

truncation that can be seen in the figure, where the color fringes drop out (left side of the left spot and 

right side of the right spot). 
NASA/TM—2013-217810 20 

Figure 6.—Probability density distribution (PDD) plots of the orientation of critical flaws on the unit sphere 

for a 100-MPa uniaxial tensile load for an isotropic material. The center of the “spot” indicates the highest 

likelihood of orientation. This tests the consistency of calculated results for loading in various directions. 

Here the Weibull modulus, mV = 10.0, and the Shetty shear sensitivity, C = 1.5. (a) Load in the x direction. (b) Load in the y direction. (c) Load in the z direction. (d) PDD in terms of magnitude normalized to 

radians for angles  and  of the load in the z direction. 
NASA/TM—2013-217810 21 

Figure 7.—Probability density distribution (PDD) plots of the orientation of critical flaws on the unit 

sphere for an isotropic material with a 100-MPa uniaxial tensile load applied in the y direction, 

showing the effect of Weibull modulus, mV, and Shetty shear sensitivity, C . (a) mV = 5.0 and 

C = 1.5. (b) mV = 5.0 and C = 100.0. (c) mV = 10.0 and C = 1.5. (d) mV = 10.0 and 

C = 100.0. (e) mV = 20.0 and C = 1.5. (f) mV = 20.0 and C = 100.0. 
NASA/TM—2013-217810 22 

Figure 8.—Probability density distribution (PDD) plots of the orientation of critical flaws on the unit sphere for 

an isotropic material for a yz = 100-MPa shear stress showing the effect of Shetty shear sensitivity, C . 

(a) Weibull modulus, mV = 10.0, and C = 100.0. (b) PDD in terms of magnitude for mV = 10.0 and 

C = 100.0. (c) mV = 10.0 and C = 1.5. (d) PDD in terms of magnitude for mV = 10.0 and C = 1.5. 

(e) mV = 10.0 and C = 0.82. (f) PDD in terms of magnitude for mV = 10.0 and C = 0.82. 
NASA/TM—2013-217810 23 

Figure 9.—Probability density distribution (PDD) plots of the orientation of critical flaws on the unit sphere for 

an isotropic material. (a) Equibiaxial tensile stress applied in the y and z directions, x = y = 100 MPa. 

(b) PDD in terms of magnitude of equibiaxial tensile stress. (c) Uniaxial compressive stress applied in the 

x direction, x = 500 MPa. (d) PDD in terms of magnitude of uniaxial compressive stress. For the tensile 

failure mode, (a) and (b), the Weibull modulus, mV = 10.0; oV = 100 MPa (MPam1/10), C = 1.5 (Eq. (11)), 

 = 0.22; and the probability of failure of material volume, PfV = 0.6321. For the compressive failure mode, 

(c) and (d), using Equation (12), mCV = 10, the Weibull scale parameter for compression, oCV = 500 MPa 

(MPam1/10), and PfV = 0.6321. 

Figures 9(a) and (b) show an example of the unit sphere PDD for critical flaw orientation for an 

equibiaxial tensile stress state applied in the y and z axial directions (x = 0; y = z). This plot was 

generated for a mildly shear sensitive material (C = 1.5). The critical flaws have equal likelihood to be at 

any angle  in the y-z plane. Figures 9(c) and (d) show an example of the unit sphere PDD for the critical 

flaw orientation for a uniaxial compressive stress state for the compressive stress failure criterion using 

Equation (12). The critical flaw planes are oriented at 45 to the direction of applied loading—the 

x direction in this case. 

3.2 Transversely Isotropic Material Critical Flaw Orientation Probability Density 

Distribution for Off-Axis Uniaxial Load 

In this section, some examples of the critical flaw orientation PDD are shown for the critical strength 

anisotropy model of Section 2.3.2. Results for the flaw orientation anisotropy model of Section 2.3.1 are 

not shown (because plotted results would appear similar to that of the critical strength anisotropy model). 
NASA/TM—2013-217810 24 

Figure 10 is an example of the unit sphere PDD for critical flaw orientation for a 2:1 transversely 

isotropic strength response. This model assumes critical strength anisotropy and uses Equation (30) to 

define the weaker directions of strength (the T distribution). In this case, a 50-MPa uniaxial tensile load 

(in the y material direction) is applied in the plane with the weakest material strength (the y-z plane) for 

PfV = 0.6321. Figures 10(a) and (b) show the unit sphere PDD for the critical flaw orientation for a shearinsensitive (C = 100.0) response, and Figures 10(c) and (d) show the results for a shear-sensitive response 

(C = 1.1236, where KIc = (KIIc or KIIIc) for  = 0.22), where the mode I strength of the flaw and the 

mode II (or mode III) of the flaw are equal. The strength response is for a material that has a more gradual 

critical strength transition from the strong to the weaker material strength direction (T = /2 radians (90) 

and T = 1). As would be expected, a more oval or elliptical spot is observed in comparison to the 

isotropic material in Figure 6(b) because of the effect of the weaker y-z material plane. 

Figure 10.—Probability density distribution (PDD) plots of the orientation of critical flaws on the unit sphere for a 

shear-insensitive and shear-sensitive transversely isotropic material. This is for a 50-MPa uniaxial tensile load 

(in the y material direction) applied along the plane with the weakest material strength (the y-z plane). The strength 

ratio is 2:1, where critical mode I stress-intensity factor, KIc, anisotropy functions representing the transverse 

(equatorial-belt) distribution, T = /2 radians (90) and T = 1. (a) Plot on unit sphere. (b) PDD in terms of magnitude for Shetty shear sensitivity, C = 100.0 (shear-insensitive response), and parameter in KIc anisotropy function, 

rT = 2.21. (c) Plot on unit sphere. (d) PDD in terms of magnitude for C = 1.1236 (shear-sensitive response) and 

T = 2.83. In both cases, the Weibull modulus, mV = 10; the Weibull scale parameter, oV = 100 MPa (MPam1/10); 

and the probability of failure of material volume, PfV = 0.6321. 
NASA/TM—2013-217810 25 

Figure 11 is an example of the unit sphere PDD for the critical flaw orientation for a 3:1 transversely 

isotropic strength response for a uniaxial load offset 15 from the strong material direction ( = 15 and 

 = 0). This model also uses Equation (30) to define the weaker directions of strength (the T distribution); however, the strength response is for a material that has a more abrupt critical strength transition 

from the strong to the weaker material strength direction (T = 0.08727 radians (5), T = 0, and rT = 3.56). 

The load is a uniaxial tensile stress of 99.3 MPa for a PfV = 0.6321. A relatively mild shear sensitivity 

(C = 1.4) is assumed. The parameter T = 0 corresponds to a uniform distribution over the equatorial 

belt where KIc is low relative to the rest of the unit sphere. This modeling scenario is arbitrary but may be 

considered to model the effect of a weak interface, such as between a fiber and matrix or within the 

matrix constituent, or possibly represent anisotropy in thin-film coatings where interfacial layers and 

anisotropically oriented microstructures are present. In addition to the PDD plot on the unit sphere 

(Fig. 11(a)), a numerical plot of the PDD versus angle of orientation of the flaw normal is shown 

(Fig. 11(b)). Contrasting the unit-sphere plots of Figures 10(a) and (c) with Figure 11(a), shows the 

effect of parameter T as a narrower strip (the smaller extent of angle ), where the critical flaw has a 

high likelihood of being centered on  = 90 and  = 0. Figure 11 also shows the relative sensitivity 

of strength anisotropy with critical flaw angle. When the offset uniaxial load is only 15, the probability 

that the critical flaw will be oriented about  = 90 and  = 0 becomes significant. This leads to the 

question of how anisotropy affects the orientation of the most critical flaw that is predicted with this 

model. 

Figures 12 to 15 were prepared to help investigate that question or to at least demonstrate that 

consequence with the unit-sphere model. These graphs plot the angle of highest probability (the 

orientation where the critical flaw PDD value is maximum) versus the angle of a uniaxial tensile load 

applied at some angle . These figures were prepared using mV = 100 MPam1/10 and PfV = 0.6321. The 

strongest material direction is at  = 0, and the weakest material direction is at  = 90. It bears 

Figure 11.—Probability density distribution (PDD) of critical flaw normals for a 3:1 strength response for a transversely isotropic material. The load is a uniaxial tensile stress of 99.3 MPa applied at a 15 offset from the x axis 

( = 15 and  = 0). A relatively mild shear sensitivity (C = 1.4) is assumed. (a) Plot on unit sphere. (b) PDD in 

terms of magnitude for Weibull modulus, mV = 10; Weibull scale parameter, oV = 100 MPam1/10; Shetty shear 

sensitivity, C = 1.4; Poisson’s ratio,  = 0.22; parameters in the KIc anisotropy function representing the transverse 

(equatorial-belt) distribution: T = 0.08727 radians (5), T = 0, and rT = 3.563; probability of failure of material 

volume, PfV = 0.6321. 

 
NASA/TM—2013-217810 26 

repeating that this is a prediction of the critically oriented flaw and not the direction of the subsequent 

crack propagation. The figures illustrate the effect of shear sensitivity C , the abruptness of the critical 

strength transition from strong to weak response using parameter T, and the strong direction to weak 

direction strength ratio modeled with parameter rT (see Table I). This line of inquiry was motivated by 

the observation that in CMCs the angle of microcracking for offset uniaxial tensile loads tends to be 

perpendicular to the direction of applied loading (e.g., Cady et al. (1995) and Lynch and Evans (1996)). 

However, this is not universally the case, as shown in Figure 16 for a 0/90 CMC laminate with a 

30-offset loading. In that case, the matrix cracks tended to be misaligned 12 to 15 with the perpendicular to the loading direction and approximately aligned symmetrically 45 from either fiber axis. 

Figure 16 shows that adjacent plies can influence the microcracking pattern. This suggests that it may 

be plausible to consider the influence of the laminate as a whole and to construct anisotropic unit-sphere 

matrix failure criteria appropriately. However, at the microscale level, CMC failure behavior involves 

more complicated mechanisms including fiber interface debonding and sliding resistances. Therefore, for 

a CMC, the anisotropic unit-sphere material failure criterion may prove useful at the macroscale of the 

composite material or alternatively at the microscale of the individual constituents involving the matrix 

with the fiber where multiscale modeling would be requisite. At the macroscopic scale, Figures 12 to 15 

would be appropriate for a unidirectional laminate where all the fibers are aligned in the same direction. 

When the individual plies are not aligned with one another, Equations (29) and (30) would have to be 

suitably modified. 

Figure 12.—Most probable orientation of the critical flaw relative 

to the offset angle of an applied uniaxial tensile load (at angle 

 on the unit sphere) for various anisotropic strength ratios. 

This example is for a mild critical strength transition (parameters in the critical mode I stress-intensity factor, KIc, anisotropy function representing the equatorial-belt distribution, 

T = /2 radians (90) and T = 1) and high shear sensitivity 

C = 1.1236 (mode I and mode II strength are equal). The 

strongest material direction is at  = 0, and the weakest 

material direction is at  = 90. PDD, probability density 

distribution. 
NASA/TM—2013-217810 27 

Figure 13.—Most probable orientation of the critical flaw relative 

to the offset angle of an applied uniaxial tensile load (at 

angle  on the unit sphere) for various anisotropic strength 

ratios. This example is for a mild critical strength transition 

(parameters in the critical mode I stress-intensity factor, KIc, 

anisotropy function representing the equatorial-belt distribution, T = /2 radians (90) and T = 1) and shear insensitivity C = 100.0. The strongest material direction is at  = 0, 

and the weakest material direction is at  = 90. PDD, 

probability density distribution. 

Figure 12 shows the most probable orientation of the critical flaw relative to the offset angle for an 

applied uniaxial load for a mild mode I critical strength transition (T = /2 radians (90) and T = 1) and 

high shear sensitivity C = 1.1236 (where mode I and mode II strength are equal). Figure 13 shows the 

results for the case of shear insensitivity C = 100.0 (with all other parameters the same). The mild critical 

strength transition is representative of materials such as graphite. Figures 12 and 13 indicate that the 

critical flaw angle deviates from the angle of offset load depending on the level of shear sensitivity and 

the degree of strength anisotropy. Figure 12 indicates that, when the uniaxial tensile load is applied in the 

strong material direction ( = 0), the critical flaw is already oriented at a significant angle of  away 

from the direction of applied load. This contrasts with Figure 13, where this angle is smaller—particularly 

for lower values of strength ratio, and tends to zero at  = 0. However, this angle of deviation increases 

rapidly as the angle of offset load increases from  = 0 for the higher levels of strength ratio. Figure 13 

shows that the gradual critical strength transition from the strong material direction to the weak material 

direction strongly influences the orientation of the critical flaw for an offset load, and Figure 12 shows the 

additional influence of shear sensitivity combined with the mode I critical strength anisotropy. 
NASA/TM—2013-217810 28 

Figure 14.—Most probable orientation of the critical flaw relative 

to the offset angle of an applied uniaxial tensile load (at angle 

 on the unit sphere) for various anisotropic strength ratios. 

This example is for a strong critical strength transition (parameters in the critical mode I stress-intensity factor, KIc, anisotropy function representing the equatorial-belt distribution, 

T = 0.1745 radians (10) and T = 1) and high shear sensitivity C = 1.1236 (mode I and mode II strength are equal). 

The strongest material direction is at  = 0, and the weakest 

material direction is at  = 90. PDD, probability density 

distribution. 

Figure 14 shows the most probable orientation of the critical flaw relative to the offset angle for an 

applied uniaxial load for an abrupt mode I critical strength transition (T = 0.1745 radians (10) and 

T = 1) and high shear sensitivity, C = 1.1236 (where mode I and mode II strength are equal). Figure 15 

shows the results for the case of shear insensitivity C = 100.0 (with all other parameters the same). The 

abrupt mode I critical strength transition is intended to be representative of materials where an interface or 

inhomogeneity or abrupt shift in microstructure is present such as an interfacial boundary between the 

matrix and a fiber coating or the fiber itself, or between different deposition methods and/or materials on 

thin-film brittle coatings. 

Figures 14 and 15 show the abrupt transition of the orientation of the critical flaw to the weak 

material direction  = 90 as the angle of the offset load increases. The abruptness of this transition can 

also be seen in Figure 11, where the offset angle was chosen near to where this transition occurs. Both 

figures indicate expected isotropic behavior prior to this transition where the most probable angle of the 

critical flaw is coincident with the angle of the applied loading (as also seen in Fig. 6 for the different 

directions of applied loading). In Figure 14, the shear sensitivity shifts the angle of where this transition 

occurs to smaller values of  in comparison to the shear insensitive condition of Figure 15. As the ratio of 

the strength anisotropy increases, it also shifts this transition to progressively smaller values of . 
NASA/TM—2013-217810 29 

A comparison of Figures 12 and 13 with Figures 14 and 15 shows the effect of parameter T on the 

most probable angle of the critical flaw. Even for a milder transition of the mode I critical strength 

(T = /2 radians (90)), the most probable orientation for the critical flaw trends toward  = 90 as the 

strength ratio increases. Overall, Figures 12 to 15 show the strong trend that the most probable orientation 

for the critical flaw is to shift toward the weaker material direction ( = 90) as the strength anisotropy 

ratio increases. 

Figure 15.—Most probable orientation of the critical flaw relative 

to the offset angle of an applied uniaxial tensile load (at angle 

 on the unit sphere) for various anisotropic strength ratios. 

This example is for a strong critical strength transition (parameters in the critical mode I stress-intensity factor, KIc, anisotropy function representing the equatorial-belt distribution, 

T = 0.1745 radians (10) and T = 1) and shear insensitivity 

C = 100.0. The strongest material direction is at  = 0, and 

the weakest material direction is at  = 90. PDD, probability 

density distribution. 
NASA/TM—2013-217810 30 

Figure 16.—Matrix cracking pattern from 30 off-axis loading of 

a 0/90 ceramic matrix composite laminate from Lynch and 

Evans (1996). Matrix cracks appear to be 12 to 15 misaligned with the perpendicular to the loading direction, where 

a is the direction of loading, b is perpendicular to the loading 

direction, c is the fiber axis of the exposed ply, d is the fiber 

axis of the underlying ply, and e is 45 to either fiber axis. 

4.0 Summary and Conclusions 

A generalized equation was developed to describe the probability density function for the distribution 

of critical flaws from the Batdorf unit-sphere model for multiaxial loading. This formulation is for isotropic and anisotropic (transversely isotropic) brittle materials. The scenarios considered for anisotropic 

strength response were for (1) flaw orientation anisotropy, whereby a preexisting microcrack has a higher 

likelihood of being oriented in one direction than another direction, and (2) critical strength or fracture 

toughness anisotropy, where the level of critical strength or fracture toughness for mode I crack propagation changes with regard to the orientation of the microstructure. 

A numerical algorithm was developed, and the results were plotted in MATLAB (MathWorks) for the 

unit-sphere probability density distribution (PDD) for various multiaxial stress states for an isotropic 

material. The results demonstrated that the PDD depends on the stress state, shear sensitivity of the flaws, 

and the Weibull modulus. Lower Weibull modulus and higher shear sensitivity meant higher dispersion of 

the distribution of critical flaws for a uniaxial stress state. For a pure shear stress state, when the shear 

sensitivity of the flaws increased, the orientation distribution of the critical flaws bifurcated into the two 

preferred angles of highest probability. This report also showed PDD results for an equibiaxial stress state 

and a uniaxial compressive stress state. 
NASA/TM—2013-217810 31 

For transversely isotropic strength response, an example was shown of the PDD for a uniaxial tensile 

load applied at an angle offset from the strongest material direction. The angle of highest probability for 

the critical flaw also was investigated versus the offset loading angle for various anisotropic strength 

ratios and abruptness of the strength dropoff for the mode I critical strength anisotropy model. These 

results indicate that when there was an abrupt dropoff in mode I critical strength, there was also an abrupt 

transition of the angle of the critical flaw versus the offset angle of the applied uniaxial load such that the 

critical flaw angle became normal to the weakest material direction. These results also were sensitive to 

the level of shear sensitivity. Even for a more gradual dropoff in mode I critical strength (versus 

orientation) the most probable critical flaw angle tended to align with the weakest material direction as 

the loading offset angle increased and as the anisotropic strength ratio increased. Overall, there was a 

strong tendency for the most probable angle of the critical flaw to quickly align with the weakest material 

direction as the offset loading angle increased and as anisotropic strength ratio increased, even with a 

more gradual transition of mode I critical strength anisotropy versus orientation on the unit sphere. 

If one knows the orientation distribution of critical flaws (and the subsequent direction of crack 

propagation), the anisotropic stiffness degradation (the anisotropic elastic constants associated with the 

damaged material) can be determined. The author anticipates that this will be helpful for a follow-on 

phase of this effort, not described here, of enabling the unit-sphere failure criterion methodology to work 

with NASA’s micromechanics analysis code/generalized method of cells (MAC/GMC, Bednarcyk and 

Arnold (2002)). This incorporation will allow the full exercise of the unit-sphere methodology, including 

incremental time/load steps and fatigue analysis (as described in Nemeth, Noel N.; Gyekenyesi, J.P.; and 

Jadaan, Osama M. (2005): Lifetime Reliability Prediction of Ceramic Structures Under Transient 

Thermomechanical Loads. NASA/TP—2005-212505. http://ntrs.nasa.gov/), to predict the durability 

(strength and lifetime) of composite laminates and woven composite structures. 

The development of the anisotropic unit-sphere methodology was an attempt to provide an improved 

mechanistic basis to the problem of predicting composite strength under multiaxial loading in comparison 

to the phenomenologically based polynomial formulations such as Tsai-Wu, Tsai-Hill, and Hashin. The 

interesting consequence of this work was that the PDD on the unit sphere of the orientation of criticalstrength-controlling flaws could also be predicted. That formulation and the demonstration of that 

capability was the subject of this report. 

This technology has general relevance to the anisotropic strength response of brittle materials. This 

may include thin-film coatings such as environmental barrier coatings where tensile stresses play a 

significant role in coating integrity. The need for additional enhancements is anticipated, such as the 

determination of the angle of noncoplanar crack growth and the development of the anisotropic stiffness 

matrix for a damaged material. 

 

NASA/TM—2013-217810 33 

Appendix—Symbols 

A area 

dA infinitesimal area on the surface of a unit radius sphere 

a crack length or crack radius 

ac critical crack length or crack radius 

C Shetty shear-sensitivity coefficient 

c constant 

c normalizing constant 

(x) discrete density function of random variable X where P(X = x) 

fIc  ,  normalized anisotropy function of KIc or Ic as a function of angles  and 

H Heaviside step function 

KI mode I stress-intensity factor 

KIc critical mode I stress-intensity factor 

kBV Batdorf crack density coefficient 

kBV Batdorf uniaxial stress-state normality constant for the volume-flaw 

 failure mode (also known as the Batdorf normalized crack-density coefficient) 

kwV uniaxial Weibull crack-density coefficient 

mCV Weibull modulus for the volume-flaw failure mode for compressive stress states 

mV Weibull modulus for the volume-flaw failure mode 

n number of incremental links; number of failure modes 

P(X = x) probability that a discrete real-valued random variable X equals a possible value x

Pf probability of failure (Pf = 1– Ps) 

PfV probability of failure of material volume 

PfV probability of failure of a crack with a strength between Ic and (Ic + Ic) in V 

(PfV)i probability of failure of the i

th link 

Ps reliability or probability of survival (Ps = 1– Pf) 

PsV probability of survival of material volume V 

(PsV)i probability of survival of the i

th link 

P1V probability of existence of a crack with strength between Ic and (Ic + Ic) in an 

incremental volume 

P2V probability that a crack of critical strength will be oriented in a particular direction 

such that it will grow and cause failure 

p success probability 
NASA/TM—2013-217810 34 

(, ) probability that a crack of critical stressIeqc is oriented in the range between α and 

( + d) and between  and ( + d) 

c(, , ) probability that the critical flaw (the failure-initiating flaw) is oriented in the range 

between  and ( + d) and between and ( + d) under stress state Σ

(, ) probability density distribution (PDD) function of (, )

(, , ) PDD function of the orientation of critical flaw normals from the applied stress state 

 given by angles  and 

rL constant (ratio) or parameter in KIc anisotropy function for polar-cap distribution 

rT constant (ratio) or parameter in KIc anisotropy function for equatorial-belt distribution 

s slope of a line 

V volume 

V incremental volume 

Vi incremental volume associated with the i

th link 

X discrete real-valued random variable 

x real valued number; possible value of discrete real-valued random variable X

x, y, z location in the body of the structure; Cartesian coordinates 

Y crack-shape geometry factor 

 orientation angles or angular coordinates, where the direction normal to the plane of 

the microcrack is specified by the radial line defined by  and  in stress space 

 constant or parameter in KIc anisotropy function representing the exponent of the sine 

or cosine function 

L constant or parameter in KIc anisotropy function representing the longitudinal (polarcap) distribution 

T constant or parameter in KIc anisotropy function representing the transverse 

(equatorial-belt) distribution 

()  as a function of angle  describing the anisotropy of flaw orientation 

(, ) function describing the anisotropy of the flaw orientation where the normal direction 

to the flaw plane is given by angles  and 

V crack-density function 

V(, ) crack-density function dependent on angles  and 

V() crack-density function for applied uniaxial stress, where  is the number of flaws per 

unit volume with strength equal to or less than 

V(Ieqc) crack-density function for equivalent mode I strength, Ieqc, of a flaw: number of 

flaws per unit volume with strength equal to or less than Ieqc

 angle of crack propagation relative to the initial plane of the microcrack 

 constant or parameter in flaw-orientation anisotropy function representing one-half of 

the total angular extent of the anisotropy distribution 
NASA/TM—2013-217810 35 

L constant or parameter in flaw-orientation anisotropy function representing the 

longitudinal (polar-cap) distribution 

T constant or parameter in flaw-orientation anisotropy function or in KIc anisotropy 

function representing the transverse (equatorial-belt) distribution 

 positive number 

 Poisson’s ratio 

L constant or parameter in KIc anisotropy function representing the longitudinal (polarcap) distribution 

T constant or parameter in KIc anisotropy function representing the transverse 

(equatorial-belt) distribution 

 summation function; applied far-field multiaxial stress state 

 applied uniaxial stress 

Ic critical mode I strength 

I ,x I ,y I ,z , ,  c  c  c orthogonal critical strength components normalized byIc,max

Ieq equivalent, or effective, stress 

Ieq,max maximum value of Ieq over the unit sphere from the applied multiaxial stress 

Ieq(x, y, z, , ) equivalent, or effective, stress as a function of location x, y, z and orientation , 

Ieqc critical equivalent, or effective, stress 

Ieqc(x, y, z, , ) mode I far-field strength of a flaw located at coordinates x, y, and z and oriented at 

angles  and 

i local uniaxial stress 

n applied far-field stress component normal to a crack face 

o Weibull scale parameter 

oCV Weibull scale parameter for compression for the volume-flaw failure mode 

normalized to unit volume 

oV Weibull scale parameter for the volume-flaw failure mode normalized to unit volume 

uV threshold stress parameter 

x, y, z orthogonal stress components expressed relative to a global coordinate system 

 shear stress acting on the oblique plane whose normal is determined by angles  and 

, which represents the applied far-field shear stress on a crack face 

xy, yz, zx shear stress components expressed relative to a global coordinate system 

 constant or parameter in flaw-orientation anisotropy function representing the 

exponent of the sine or cosine function 

L constant or parameter in flaw-orientation anisotropy function representing the 

longitudinal (polar-cap) distribution 
NASA/TM—2013-217810 36 

T constant or parameter in flaw-orientation anisotropy function representing the 

transverse (equatorial-belt) distribution 

(, Ieqc) area of a solid angle projected onto a unit radius sphere in three-dimensional stress 

space for which Ieq  Ieqc from an applied multiaxial stress state 

Superscript 

- normalized 

Subscripts 

I, II, III mode I, II, or III 

i ith value or i

th term 

L longitudinal 

max maximum 

n integer; number of links; normal 

T transverse 

V volume or a volume-based property (e.g., indicates volume-flaw analysis) 

Definitions 

L longitudinal 

mode I crack-opening mode 

mode II crack-sliding mode (in-plane shear) 

mode III crack-tearing mode (out-of-plane shear) 

T transverse 

 
NASA/TM—2013-217810 37 

References 

Aboudi, Jacob (1995): Micromechanical Analysis of Thermo-Inelastic Multiphase Short-Fiber Composites. 

Compos. Eng., vol. 5, no. 7, pp. 839–850. 

Aboudi, Jacob; Pindera, Marek-Jerzy; and Arnold, Steven M. (2003): Higher-Order Theory for Periodic 

Multiphase Materials With Inelastic Phases. Int. J. Plast., vol. 19, pp. 805–847. 

Batdorf, S.B. (1973): A Statistical Theory for the Fracture of Transversely Isotropic Brittle Materials of 

Moderate Anisotropy. Air Force Report SAMSO–TR–73–361, Aerospace Report TR–0074 

(4450–76)–1, NTIS AD–770982. 

Batdorf, S.B. (1978a): Fracture Statistics of Polyaxial Stress States. Fracture Mechanics, Nicholas Perrone, 

et al., eds., University Press of Virginia, Charlottesville, VA, pp. 579591. 

Batdorf, S.B. (1978b): New Light on Weibull Theory. Nucl. Eng. Des., vol. 47, pp. 267–272. 

Batdorf, S.B.; and Crose, J.G. (1974): A Statistical Theory for the Fracture of Brittle Structures Subjected to 

Nonuniform Polyaxial Stresses. J. Appl. Mech., vol. 41, no. 2, pp. 459464. 

Batdorf, S.B.; and Heinisch, H.L., Jr. (1978): Weakest Link Theory Reformulated for Arbitrary Fracture 

Criterion. J. Am. Ceram. Soc., vol. 61, nos. 7‒8, pp. 355358. 

Bednarcyk, Brett A.; and Arnold, Steven M. (2002): MAC/GMC 4.0 User’s Manual—Keywords Manual. 

NASA/TM—212077/VOL2. http://ntrs.nasa.gov/ 

Buch, J.D.; Crose, J.G.; and Robinson, E.Y. (1977): Failure Criteria in Graphite Program. AFML–TR– 

77–16. Available from the Air Force Material Laboratory. 

Cady, Carl; Heredia, Fernando E.; and Evans, Anthony G. (1995): In-Plane Mechanical Properties of 

Several Ceramic-Matrix Composites. J. Am. Ceram. Soc., vol. 78, no. 8, pp. 2065‒2078. 

Erdogan, F.; and Sih, G.C. (1963): On the Crack Extension in Plates Under Plane Loading and Transverse 

Shear. J. Basic Eng., vols. 85–86, pp. 519525. 

Evans, A.G. (1978): A General Approach for the Statistical Analysis of Multiaxial Fracture. J. Amer. 

Ceram. Soc., vol. 61, nos. 7–8, pp. 302–308. 

Hellen, T.K.; and Blackburn, W.S. (1975): The Calculation of Stress Intensity Factors for Combined Tensile 

and Shear Loading. Int. J. Fract., vol. 11, no. 4, pp. 605617. 

Hoel, Paul G.; Port, Sydney C.; and Stone, C.J. (1971): Introduction to Probability Theory. Houghton 

Mifflin, Boston, MA, pp. 65 and 131. 

Ichikawa, M. (1991): Proposal of an Approximate Analytical Expression of Maximum Energy Release Rate 

of a Mixed Mode Crack in Relation to Reliability Evaluation of Ceramic Components. J. Soc. Mat. Sci. 

Jpn., vol. 40, pp. 224227. 

Lutz, G. (2006): The Puck Theory of Failure in Laminates in the Context of the New Guideline VDI2014 

Part 3. Conference on Damage in Composite Materials 2006, Stuttgart, Germany. 

Lynch, Christopher S.; and Evans, Anthony G. (1996): Effects of Off-Axis Loading on the Tensile Behavior 

of a Ceramic-Matrix Composite. J. Am. Ceram. Soc., vol. 79, no. 12, pp. 3113‒3123. 

Matsuo, Yotaro (1981): A Probabilistic Analysis of the Brittle Fracture Loci Under Bi-Axial Stress State: 

1st Report; In the Case of Tension Being Dominant. Bull. JSME, vol. 24, no. 188, pp. 290–294. 

Nemeth, Noel N. (2013): Unit-Sphere Multiaxial Stochastic-Strength Model Applied to Anisotropic and 

Composite Materials. NASA/TP—2013-217749, 2013. http://ntrs.nasa.gov/ 

Nemeth, Noel N.; and Bratton, Robert L. (2011): Statistical Models of Fracture Relevant to Nuclear-Grade 

Graphite: Review and Recommendations. NASA/TM—2011-215805. http://ntrs.nasa.gov/ 

Nemeth, Noel N., et al. (2003): CARES/Life Ceramics Analysis and Reliability Evaluation of Structures 

Life Prediction Program. NASA/TM—2003-106316. http://ntrs.nasa.gov/ 

Nemeth, Noel N.; Gyekenyesi, J.P.; and Jadaan, Osama M. (2005): Lifetime Reliability Prediction of 

Ceramic Structures Under Transient Thermomechanical Loads. NASA/TP—2005-212505. 

http://ntrs.nasa.gov/ 

Nemeth, Noel N.; Manderscheid, Jane M.; and Gyekenyesi, John P. (1990): Ceramics Analysis and 

Reliability Evaluation of Structures (CARES). Users and Programmers Manual. NASA TP2916. 

http://ntrs.nasa.gov/ 
NASA/TM—2013-217810 38 

Shetty, D.K. (1987): Mixed-Mode Fracture Criteria for Reliability Analysis and Design With Structural 

Ceramics. J. Eng. Gas Turbines Power, vol. 109, no. 3, pp. 282289. 

Weibull, Waloddi (1939): A Statistical Theory of the Strength of Materials. Ingeniorsvetenskapsakademiens 

Handlingar, no. 151.

REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188 

The 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 of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this 

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)

01-09-2013 

2. REPORT TYPE

Technical Memorandum 

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

Unit-Sphere Anisotropic Multiaxial Stochastic-Strength Model Probability Density 

Distribution for the Orientation of Critical Flaws 

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Nemeth, Noel, N. 

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

WBS 984754.02.07.03.16.03.02 

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

National Aeronautics and Space Administration 

John H. Glenn Research Center at Lewis Field 

Cleveland, Ohio 44135-3191 

8. PERFORMING ORGANIZATION

 REPORT NUMBER

E-18546-2 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

Washington, DC 20546-0001 

10. SPONSORING/MONITOR'S

 ACRONYM(S)

NASA 

11. SPONSORING/MONITORING

 REPORT NUMBER

NASA/TM-2013-217810 

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited 

Subject Category: 37 

Available electronically at http://www.sti.nasa.gov 

This publication is available from the NASA Center for AeroSpace Information, 443-757-5802

13. SUPPLEMENTARY NOTES

An Erratum was added to this report March 2014. 

14. ABSTRACT

Models that predict the failure probability of monolithic glass and ceramic components under multiaxial loading have been developed by 

authors such as Batdorf, Evans, and Matsuo. These “unit-sphere” failure models assume that the strength-controlling flaws are randomly 

oriented, noninteracting planar microcracks of specified geometry but of variable size. This report develops a formulation to describe the 

probability density distribution of the orientation of critical strength-controlling flaws that results from an applied load. This distribution is a 

function of the multiaxial stress state, the shear sensitivity of the flaws, the Weibull modulus, and the strength anisotropy. Examples are 

provided showing the predicted response on the unit sphere for various stress states for isotropic and transversely isotropic (anisotropic) 

materials--including the most probable orientation of critical flaws for offset uniaxial loads with strength anisotropy. The author anticipates 

that this information could be used to determine anisotropic stiffness degradation or anisotropic damage evolution for individual brittle (or 

quasi-brittle) composite material constituents within finite element or micromechanics-based software.

15. SUBJECT TERMS

Batdorf; Failure probability; Multiaxial; Strength; Probability density function; Flaw orientation; Unit sphere; Anisotropy; Transverse 

isotropy 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF

 ABSTRACT

UU 

18. NUMBER

 OF 

 PAGES

46 

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email:help@sti.nasa.gov) 

a. REPORT

U 

b. ABSTRACT

U 

c. THIS 

PAGE

U 

19b. TELEPHONE NUMBER (include area code)

443-757-5802 

Standard Form 298 (Rev. 8-98)

Prescribed by ANSI Std. Z39-18