National Aeronautics and

Space Administration

IS04

George C. Marshall Space Flight Center

Marshall Space Flight Center, Alabama

35812
The NASA STI Program Office...in Profile

Since its founding, NASA has been dedicated to

the advancement of aeronautics and space

science. The NASA Scientific and Technical 

Information (STI) Program Office plays a key

part in helping NASA maintain this important

role.

The NASA STI Program Office is operated by Langley 

Research Center, the lead center for NASA’s scientific 

and technical information. The NASA STI Program 

Office provides access to the NASA STI Database, the 

largest collection of aeronautical and space science 

STI in the world. The Program Office is also NASA’s 

institutional mechanism for disseminating the results 

of its research and development activities. These 

results are published by NASA in the NASA STI 

Report Series, which includes the following report 

types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major significant 

phase of research that present the results of 

NASA programs and include extensive data 

or theoretical analysis. Includes compilations 

of significant scientific and technical data 

and information deemed to be of continuing 

reference value. NASA’s counterpart of peerreviewed formal professional papers but has 

less stringent limitations on manuscript length 

and extent of graphic presentations.

• TECHNICAL MEMORANDUM. Scientific 

and technical findings 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. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientific and technical 

conferences, symposia, seminars, or other 

meetings sponsored or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientific, 

technical, or historical information from 

NASA programs, projects, and mission, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. 

English-language translations of foreign 

scientific and technical material pertinent to 

NASA’s mission.

Specialized services that complement the STI Program 

Office’s diverse offerings include creating custom 

thesauri, building customized databases, organizing 

and publishing research results...even providing 

videos.

For more information about the NASA STI Program 

Office, see the following:

• Access the NASA STI Program Home Page at 

http://www.sti.nasa.gov

• E-mail your question via the Internet to 

help@sti.nasa.gov

• Fax your question to the NASA Access Help 

Desk at 301–621–0134

• Telephone the NASA Access Help Desk at 

301–621–0390

• Write to:

NASA Access Help Desk

NASA Center for AeroSpace Information

7121 Standard Drive

Hanover, MD 21076–1320

301–621–0390


NASA/CR—2005–214006

The M-Integral for Computing 

Stress Intensity Factors in Generally 

Anisotropic Materials

P.A. Wawrzynek and B.J. Carter

Fracture Analysis Consultants, Inc., Ithaca, New York

L. Banks-Sills, Consultant

Tel Aviv University, Tel Aviv, Israel

July 2005

Prepared for Marshall Space Flight Center

under Contract NAS8–36801 

and sponsored by the 

Space Shuttle Main Engine Program

managed at the Marshall Space Flight Center

National Aeronautics and

Space Administration

Marshall Space Flight Center • MSFC, Alabama 35812
ii

Available from:

NASA Center for AeroSpace Information National Technical Information Service

7121 Standard Drive 5285 Port Royal Road

Hanover, MD 21076–1320 Springfield, VA 22161

301–621–0390 703–487–4650
iii

TABLE OF CONTENTS

1. INTRODUCTION ........................................................................................................................ 1

1.1 Background ............................................................................................................................ 1

1.2 Anisotropic Material Properties ............................................................................................. 2

2. CRACK-FRONT STRESS AND DISPLACEMENT FIELDS .................................................... 3

2.1 General Anisotropy ................................................................................................................ 3

2.2 Special Anisotropies ............................................................................................................... 5

3. M-INTEGRAL FORMULATION ................................................................................................ 7

3.1 The Relationship Between J and K ........................................................................................ 7

3.2 The M-Integral Expressions ................................................................................................... 8

4. M-INTEGRAL IMPLEMENTATION .......................................................................................... 13

4.1 Integration Domains and q-Functions .................................................................................... 13

4.2 Detailed Implementation Steps............................................................................................... 17

 4.2.1 Compute material constitutive matrices ..................................................................... 17

 4.2.2 Compute the crack-front coordinate system .............................................................. 18

 4.2.3 Transform nodal coordinates, displacements, and material matrices ........................ 19

 4.2.4 Determine the nature of the isotropy/anisotropy ....................................................... 19

 4.2.5 Compute anisotropic parameters ................................................................................ 20

 4.2.6 Compute the auxiliary nodal displacements .............................................................. 20

 4.2.7 Compute displacement derivatives ............................................................................ 21

 4.2.8 Compute local stresses ............................................................................................... 21

 4.2.9 Compute the contribution to the M-integral ............................................................... 22

 4.2.10 Compute the area of the virtual extension ................................................................. 23

 4.2.11 Compute the stress intensity factors .......................................................................... 23

5. PRESCRIBED DISPLACEMENT TESTS .................................................................................. 24

6. THICK PLATE ANALYSIS ......................................................................................................... 25

7. BRAZILIAN DISK ANALYSIS .................................................................................................. 29

8. SUMMARY .................................................................................................................................. 33

REFERENCES .................................................................................................................................. 34

APPENDIX A: BRAZILIAN DISK ANALYSIS RESULTS ............................................................ 35
iv

LIST OF FIGURES

1. Local crack front coordinate system used throughout this report ..................................... 3

2. The integration contour for the 2-D J-integral .................................................................. 7

3. Domain of integration for computing the 3-D M-integral ................................................. 10

4. The q function used in evaluating the 3-D M-integral ....................................................... 10

5. The domain of integration and the virtual crack extension for q(1) using one ring 

of elements ......................................................................................................................... 14

6. The domain of integration and the virtual crack extension for q(2) using one ring 

of elements ......................................................................................................................... 14

7. The domain of integration and the virtual crack extension for q(3) using one ring 

of elements ......................................................................................................................... 15

8. The domain of integration and the virtual crack extension for q(4) using one ring 

of elements ......................................................................................................................... 15

9. The domain of integration and the virtual crack extension for q(1) using two rings 

of elements ......................................................................................................................... 15

10. The domain of integration and the virtual crack extension for q(2) using two rings 

of elements ......................................................................................................................... 16

11. The domain of integration and the virtual crack extension for q(3) using two rings 

of elements ......................................................................................................................... 16

12. The domain of integration and the virtual crack extension for q(3) using two rings 

of elements ......................................................................................................................... 16

13. The definition of the local crack-front coordinates for case q(1), based on specified 

nodal coordinates ............................................................................................................... 18

14. A schematic of how the coordinates r and are determined ............................................. 21

15. The geometry used for the thick plate analysis ................................................................. 26

θ
LIST OF FIGURES (Continued)

16. The simply supported boundary conditions (a), and the mode I (b) and II 

(c) load cases used in the analyses...................................................................................... 26

17. The mesh used for the thick plate verification analyses..................................................... 27

18. A detail of the thick plate mesh showing the crack region ................................................ 28

19. The relative errors between the theoretical and numerical solutions for the stress 

 intensity factors for a thick center-cracked plate ............................................................... 28

20. The geometry of the Brazilian disk test specimen ............................................................. 29

21. A graphical illustration of the cubic material orientation relative to the local 

crack-front axes ................................................................................................................. 30

22. A typical mesh used in the disk analyses (with greatly magnified displacements) ........... 31

23. Typical Brazilian disk analysis results, where F K Wt

P a a

n

n =

( ) +

π

2π 1 2

, with n=I, II .......... 32
The NASA STI Program Office...in Profile

Since its founding, NASA has been dedicated to

the advancement of aeronautics and space

science. The NASA Scientific and Technical 

Information (STI) Program Office plays a key

part in helping NASA maintain this important

role.

The NASA STI Program Office is operated by Langley 

Research Center, the lead center for NASA’s scientific 

and technical information. The NASA STI Program 

Office provides access to the NASA STI Database, the 

largest collection of aeronautical and space science 

STI in the world. The Program Office is also NASA’s 

institutional mechanism for disseminating the results 

of its research and development activities. These 

results are published by NASA in the NASA STI 

Report Series, which includes the following report 

types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major significant 

phase of research that present the results of 

NASA programs and include extensive data 

or theoretical analysis. Includes compilations 

of significant scientific and technical data 

and information deemed to be of continuing 

reference value. NASA’s counterpart of peerreviewed formal professional papers but has 

less stringent limitations on manuscript length 

and extent of graphic presentations.

• TECHNICAL MEMORANDUM. Scientific 

and technical findings 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. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientific and technical 

conferences, symposia, seminars, or other 

meetings sponsored or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientific, 

technical, or historical information from 

NASA programs, projects, and mission, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. 

English-language translations of foreign 

scientific and technical material pertinent to 

NASA’s mission.

Specialized services that complement the STI Program 

Office’s diverse offerings include creating custom 

thesauri, building customized databases, organizing 

and publishing research results...even providing 

videos.

For more information about the NASA STI Program 

Office, see the following:

• Access the NASA STI Program Home Page at 

http://www.sti.nasa.gov

• E-mail your question via the Internet to 

help@sti.nasa.gov

• Fax your question to the NASA Access Help 

Desk at 301–621–0134

• Telephone the NASA Access Help Desk at 

301–621–0390

• Write to:

NASA Access Help Desk

NASA Center for AeroSpace Information

7121 Standard Drive

Hanover, MD 21076–1320

301–621–0390
 

CONTRACTOR REPORT 

THE M-INTEGRAL FOR COMPUTING STRESS INTENSITY FACTORS 

IN GENERALLY ANISOTROPIC MATERIALS 

1. INTRODUCTION 

The objective of this project is to develop and demonstrate a capability for computing stress 

intensity factors in generally anisotropic materials. These objectives have been met. 

The primary deliverable of this project is this report and the information it contains. In addition, 

we have delivered the source code for a subroutine that will compute stress intensity factors for 

anisotropic materials encoded in both the C and Python programming languages and made 

available a version of the FRANC3D program that incorporates this subroutine. 

Single crystal super alloys are commonly used for components in the hot sections of 

contemporary jet and rocket engines. Because these components have a uniform atomic lattice 

orientation throughout, they exhibit anisotropic material behavior. This means that stress 

intensity solutions developed for isotropic materials are not appropriate for the analysis of crack 

growth in these materials. Until now, a general numerical technique did not exist for computing 

stress intensity factors of cracks in anisotropic materials and cubic materials in particular. Such a 

capability was developed during the project and is described and demonstrated herein. 

1.1 Background 

Expressions for the stress and displacement fields near the tip of a crack in an anisotropic 

material were first published by Sih, Paris, and Irwin [1965]. Their results were used by a 

number of authors to create numerical techniques for computing stress intensity factors in such 

materials (e.g., Atluri et al. [1975], Boone et al. [1987], and Tohgo et al. [1993]). However, a 

restriction not strongly emphasized in the original or any of the subsequent papers is that these 

results are only valid for the case where the plane of analysis (typically a plane perpendicular to 

the crack front) is a plane of material symmetry. 

The expressions for the stress and displacement field near a crack front in a generally anisotropic 

material were first published by Hoenig [1982]. These results do not appear to be widely 

referenced in the literature and it does not appear that a numerical technique for computing stress 

intensity factors based on these results has been developed previously. The Hoenig paper 

contains a significant number of typographical errors. We believe that these have all been 

identified and corrected in the work reported herein. 
2 

1.2 Anisotropic Material Properties 

The primary materials of concern for the present study are single-crystal alloys. These materials 

exhibit cubic behavior, which means that they can be characterized by three properties (usually 

E, G, and v) and their orientation is given by three angles. The current implementation supports 

the more general case of arbitrarily oriented orthotropic materials, which are characterized by 

nine properties (E1, E2, E3, G12, G23, G13, ν12, ν23, and ν13) and three angles. The fully general 

case of a material with 21 independent material properties could be accommodated easily in the 

current software, but as it is rare to assume such general behavior in engineering practice, this is 

not currently supported in order to simplify the interface. 

As mentioned, three angles are required to specify the orientation of the material axes relative to 

the Cartesian axes. There are a number of conventions for this. Here we adopted the use of Euler 

angles, which are commonly used in many finite element applications. Our convention conforms 

to that used by the ANSYS software, where the first rotation, θz, is about the original z-axis. The 

second rotation, θx, is about the rotated x-axis, and the third rotation, θy, is about the rotated yaxis. In all cases, a "right-handed" rule is used to determine the sign of the rotation. 

The remainder of this report is divided into seven sections. In the following section, the 

expressions for the stress and displacement fields near crack fronts in anisotropic materials are 

developed. In Section 3, the M-integral is formulated for anisotropic materials. In Section 4, the 

details of the present implementation of the M-integral are described. In Section 5, the results of 

verification studies with prescribed nodal displacements are presented. In Section 6, the results 

of verification studies of cracks in thick center-crack plates are presented. In Section 7, the 

results of a large number of analyses of cracks in single-crystal Brazilian disk specimens are 

summarized. Finally, Section 8 summarizes the report. 
 3 

2. CRACK-FRONT STRESS AND DISPLACEMENT FIELDS 

2.1 General Anisotropy 

The stress and displacement fields in the vicinity of a crack-front in a general anisotropic 

material were presented by Hoenig [1982] based on the Lekhnitskii [1963] formalism. These are 

given by 

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎭

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎬

⎫

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎩

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎨

⎧

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ −

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ −

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

=

⎪

⎪

⎪

⎭

⎪

⎪

⎪

⎬

⎫

⎪

⎪

⎪

⎩

⎪

⎪

⎪

⎨

⎧

∑

∑

∑

∑

∑

=

−

=

−

=

−

=

−

=

−

3

1

1

3

1

1

3

1

1

3

1

1

3

1

2 1

Re

Re

Re

Re

Re

2

1

i i

i ij j

i i

i i ij j

i i

i ij j

i i

ij j

i i

i ij j

yz

zx

xy

yy

xx

Q

N K

Q

N K

Q

N K

Q

N K

Q

N K

r

λ

μ λ

μ

μ

π

σ

σ

σ

σ

σ

 (1) 

and 

 ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ = ∑=

− 3

1

1 Re 2

j

i mijNjl Kl

Qj

r

u

π

, (2) 

where Re represents the real part of the quantity in parenthesis, two repeated indices obey the 

summation convention from 1 to 3, the coordinates x, y, z, r, and θ refer to the crack coordinates 

illustrated in Figure 1, and Kj are the stress intensity factors KI, KII, and KIII. It is assumed that the 

stresses and displacements are functions only of the in-plane coordinates, r and θ. 

Figure 1. Local crack front coordinate system used throughout this report. 
4 

The quantities pi, λi, mij, Nij, and Qi, are derived from the material compliance tensor, ijkl S (where 

ij ijkl kl ε = S σ ). The overbar indicates quantities specified in the Cartesian coordinate system. 

Before extracting the quantities, it is necessary to rotate the tensor to the crack-front coordinate 

system, Sijkl, and contract it to a 6x6 matrix, Sij, where i = 1..6. It is common in many branches of 

solid mechanics to work with the contracted form of the constitutive matrix. However, there are 

two different conventions for the ordering of the elements in the matrix. Both are used in this 

work. In the theoretical Sections (2 and 3), indices of the contracted matrix are taken in the order 

{x, y, z, yz, xz, xy}. This is consistent with the Hoenig and Sih papers and much of the literature 

for orthotropic materials (principally composites). However, in the actual code and those 

portions of the report concerning the implementation, the contracted matrix indices are taken in 

the order {x, y, z, xy, yz, zx}. This is consistent with most finite element software and much of the 

finite element literature. 

For plane stress conditions, the full compliance matrix is employed. For generalized plane strain 

conditions (which are assumed in the implementation and throughout this report), a reduced 

compliance matrix is obtained from 

33

3 3

S

S S

S S i j

ij′ = ij − (3) 

where i,j = 1..6, S'ij is symmetric, and S'i3 = S'3i = 0. 

The parameters μi, i = 1..3, are the eigenvalues of the compatibility equations with positive 

imaginary parts. These are found from the characteristic sixth order polynomial equation 

( ) ( ) ( ) 0 2 l4 μ l2 μ − l3 μ = (4) 

where 

( )

() ( ) ( )

() ( ) 26 22

2

12 66

3

16

4

4 11

25 46 24

2

14 56

3

3 15

45 44

2

2 55

2 2 2

2

l S S S S S S

l S S S S S S

l S S S

= ′ − ′ + ′ + ′ − ′ + ′

= ′ − ′ + ′ + ′ + ′ − ′

= ′ − ′ + ′

μ μ μ μ μ

μ μ μ μ

μ μ μ

 (5) 

Solution of eqn. 4 leads to three pairs of complex conjugate roots. The three with positive 

imaginary parts are chosen. The parameters λi are given by 

( )

( )i

i

i l

l

μ

μ λ

2

3 = − , i = 1..3. (6) 

The parameters mij are given by 
 5 

( )

( )

( ) i i i i i

i i i i i

i i i i i

m S S S S S

m S S S S S

m S S S S S

μ μ λ μ

μ μ λ μ

μ μ λ μ

3 41 46 42 45 44

2 21 26 22 25 24

16 12 15 14

2

1 11

= ′ − ′ + ′ + ′ − ′

= ′ − ′ + ′ + ′ − ′

= ′ − ′ + ′ + ′ − ′

. (7) 

The matrix Nij is defined as 

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎣

⎡

− − −

= − − −

1 2 3

1 2 3

1 1 1

λ λ λ

Nij μ μ μ , (8) 

which has an inverse 

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎣

⎡

− − −

− − −

− − −

= −

1 2 2 1 2 1 1 2

3 1 1 3 1 3 3 1

2 3 3 2 3 2 2 3

1 1

μ λ μ λ λ λ μ μ

μ λ μ λ λ λ μ μ

μ λ μ λ λ λ μ μ

N

Nij (9) 

and |N| is the determinate of Nij. Finally, 

Qi = cosθ + μi sinθ (10) 

2.2 Special Anisotropies 

There are special anisotropies for which the above development is mathematically degenerate. 

These include monoclinic and cubic materials oriented such that z = 0 is a plane of material 

symmetry. Hoenig [1982] gives expressions for the stresses and displacements for the degenerate 

cases. However, the implementation follows the more widely referenced expressions given 

earlier by Sih et. al. [1965]. These stress and displacement expressions are 

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎭

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎬

⎫

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎩

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎨

⎧

⎥

⎦

⎤ ⎢

⎣

⎡

⎥

⎦

⎤ ⎢

⎣

⎡ −

⎥

⎦

⎤ ⎢

⎣

⎡

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ − ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

− ⎥ +

⎦

⎤ ⎢

⎣

⎡

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ − ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

−

⎥

⎦

⎤ ⎢

⎣

⎡

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ − ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

− ⎥ +

⎦

⎤ ⎢

⎣

⎡

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ − ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

−

⎥

⎦

⎤ ⎢

⎣

⎡

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ − ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

− ⎥ +

⎦

⎤ ⎢

⎣

⎡

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ − ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

−

=

⎪

⎪

⎪

⎭

⎪

⎪

⎪

⎬

⎫

⎪

⎪

⎪

⎩

⎪

⎪

⎪

⎨

⎧

3

3

3

2

2

1

1

1 2 1 2 1 2

1 2

1 1 2 2 1

2

2

1

1 2

1

2

1

2

2

2

1 1 2

1

2

2

1 2

1 2

1 Re

Re

1 Re 1 1 Re

1 1 1 Re 1 Re

1 Re Re

2

1

Q

K

Q

K

Q Q

K

Q Q

K

Q Q

K

Q Q

K

Q Q

K

Q Q

K

r

III

III

I II

I II

I II

yz

zx

xy

yy

xx

μ

μ μ

μ μ μ μ

μ μ

μ μ

μ μ

μ μ

μ μ

μ μ

μ μ

μ μ

μ μ

π

σ

σ

σ

σ

σ

 (11) 
6 

and 

( ) ( )

( ) ( )

⎪

⎪

⎪

⎪

⎭

⎪

⎪

⎪

⎪

⎬

⎫

⎪

⎪

⎪

⎪

⎩

⎪

⎪

⎪

⎪

⎨

⎧

⎥

⎦

⎤ ⎢

⎣

⎡

′ + ′

⎥

⎦

⎤ ⎢

⎣

⎡ − ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

− ⎥ +

⎦

⎤ ⎢

⎣

⎡ − ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

−

⎥

⎦

⎤ ⎢

⎣

⎡ − ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

− ⎥ +

⎦

⎤ ⎢

⎣

⎡ − ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

−

=

⎪

⎭

⎪

⎬

⎫

⎪

⎩

⎪

⎨

⎧

45 3 44

3

2 2 1 1

1 2

1 2 2 2 1 1

1 2

2 2 1 1

1 2

1 2 2 2 1 1

1 2

3

2

1

Re

1 Re 1 Re

1 Re 1 Re

2

S S

Q K

K q Q q Q K q Q q Q

K p Q p Q K p Q p Q

r

u

u

u

III

I II

I II

μ

μ μ

μ μ

μ μ

μ μ

μ μ

μ μ

π

 (12) 

In these expressions, μ1 and μ2 are the roots of the equation 

2 (2 ) 2 26 22 0 2

12 66

3

16

4 S11′ μ − S′ μ + S′ + S′ μ − S′ μ + S′ = (13) 

with positive imaginary parts. 

μ3 is the root of the equation 

2 0. 45 55

2 S44′ μ + S′ μ + S′ = (14) 

with a positive imaginary part. 

The parameters pj and qj are defined as 

j j

j j j

q S S S

p S S S

μ μ

μ μ

12 26 22

16 12

2

11

= ′ − ′ + ′

= ′ − ′ + ′

 i = 1,2 (15) 

Finally, 

Qi = cosθ + μi sinθ . (16) 
 7 

3. M-INTEGRAL FORMULATION 

3.1 The Relationship Between J and K

As a first step in formulating an M-integral for anisotropic materials, the relationship between the 

J-integral and the stress intensity factors must be established. 

The J-integral developed by Rice [1968] is given by 

ds

x

u J Wn T i

∫ x i Γ ⎟

⎠

⎞ ⎜

⎝

⎛

∂

∂ = − (17) 

where, for 2-D, the integration path Γ is illustrated in Figure 2, W is the strain energy density 

W ij ij σ ε 2

1 = (18) 

ni is the outer normal to Γ, Ti = σijnj, ui is the displacement vector, and ds is the differential arc 

length along Γ. 

Figure 2. The integration contour for the 2-D J-integral. 

By means of the crack closure integral, Hoenig [1982] developed the relationship between the Jintegral and the stress intensity factors as 

( [ ] [ ] [ ]) KI m iNij K j KII m iNij K j KIII m iNij K j J 1

3

1

1

1

2 2

1 Im Im Im − − − = − + + (19) 

where Im is the imaginary part of the bracketed quantities and summation is applied to repeated 

indices. 

For monoclinic or cubic materials oriented such that z = 0 is a symmetry plane, the relationship 

between J and K is 
8 

( ) [ ]

2 2

2 44 55 45

1

22

1 2

1 2 11

2

1 2 11 16

2

22

1 2

1 2

2 26

1

1

Im

I II III

I II

S S K K S S S K

J S S K S S K

+ ′ ′ − ′ ⎥

⎦

⎤

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

+ ′ − ′

⎢

⎢

⎣

⎡

+ + ′ − ′ ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ ′ ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ + = ′ −

μ μ

μ μ

μ μ

μ μ

μ μ

 (20) 

3.2 The M-Integral Expressions 

Yau et. al. [1980] developed the M-integral from the J-integral as a way to extract the stress 

intensity factors for the three fracture modes from the global energy release rate. To obtain the 

M-integral from eqn. 19, two solutions are assumed and superposed. This is possible if the 

material is assumed to be linear elastic. Define 

(1) (2)

(1) (2)

(1) (2)

(1) (2)

i i i

i i i

ij ij ij

ij ij ij

K K K

u u u

= +

= +

= +

= +

ε ε ε

σ σ σ

. (21) 

Substitution of the stress, strain, and displacement fields into the J-integral expression (eqn. 17) 

leads to 

 (1) (2) (1,2) J = J + J + M (22) 

where 

∫

∫

∫

Γ

Γ

Γ

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

∂

∂ − ∂∂ = −

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

∂

∂ = −

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

∂

∂ = −

ds

x

u T

x

u M W n T

ds

x

u J W n T

ds

x

u J W n T

i

i

i

i

i

i

i

i

1

(1)

(2)

1

(2)

(1)

1

(1,2) (1,2)

1

(2)

(2)

1

(2) (2)

1

(1)

(1)

1

(1) (1)

 (23) 

The interaction strain energy density is given by 

(1,2) (1) (2) (2) (1) W ij ij ij ij =σ ε =σ ε (24) 

and M(1,2) is the M-integral. 

The superposed stress intensity factors in eqn. 21 can be substituted into eqn. 19 to yield 
 9 

[ ( ) ( )

( )()

( )( )] 1 (2)

3

1 (2) (2)

3

(1)

1 (1)

1

1 (2) (2)

1

(1)

1 (1)

2

1 (2) (2)

2

(1)

2

(1) (2) 1

Im Im

Im Im

Im Im

III i ij j III i ij j

II i ij j II i ij j

I i ij j I i ij j

K m N K K m N K

K m N K K m N K

J J J K m N K K m N K

− −

− −

− −

+ +

+ +

= + − +

 (25) 

where 

[ ( ) ( ) ( )]

[ ] ( )( )( ) 1 (2)

3

1 (2) (2)

1

1 (2) (2)

2

(2)

2

(2) 1

1 (1)

3

1 (1) (1)

1

1 (1) (1)

2

(1)

2

(1) 1

Im Im Im

Im Im Im

I i ij j II i ij j III i ij j

I i ij j II i ij j III i ij j

J K m N K K m N K K m N K

J K m N K K m N K K m N K

− − −

− − −

= − + +

= − + +

. (26) 

Equating eqn. 22 and eqn. 25, and noting eqn. 23, leads to the definition of the M-integral as 

[ ( ) ( )

( )()

( )( )] 1 (2)

3

1 (2) (2)

3

(1)

1 (1)

1

1 (2) (2)

1

(1)

1 (1)

2

1 (2) (2)

2

(1)

2

1

Im Im

Im Im

Im Im

III i ij j III i ij j

II i ij j II i ij j

I i ij j I i ij j

K m N K K m N K

K m N K K m N K

M K m N K K m N K

− −

− −

− −

+ +

+ +

= − +

. (27) 

To this point, the J- and M-integrals have been presented as contour integrals. It is well known 

that in the context of finite elements the evaluation of area and volume integrals lead to more 

accurate and stable results than direct integration along a contour. Li et. al. (1985) introduced the 

2-D equivalent domain integral (EDI) form of the J-integral that transforms eqn. 17 to an area 

integral. An analogous technique can be used to formulate a 3-D equivalent domain version of 

the M-Integral where the integration takes place over a volume. The resulting expressions are 

∫

∫

∫

=

= = ∂

∂

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ − ∂

∂

+

∂

∂ =

= = ∂

∂

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ − ∂

∂ =

L q t

V

j

j

i

ij

i

ij

j

V j

k

k

k i

ij

k

A q ds

dV i j

x

q W

x

u

x

u M

dV i k j

x

q W

x

u J

.

1,2,3 1,2

1,2,3 , 1,2

1

(1,2)

1

(1)

(2)

1

(2)

(1,2) (1)

1

( )

1

( )

( ) ( )

σ σ δ

σ δ

 (28) 

where δij is the Kronecker delta and q is a function that is 1 at the stress intensity evaluation point 

and zero on the outer boundary of the domain of integration; qt is the value of the q function 

along the crack front. The domain of integration is a cylinder that encloses a portion of the crack 

front. This is illustrated in Figure 3. The q function can be thought of as a virtual crack 

extension. This is illustrated in Figure 4. 
10 

Figure 3. Domain of integration for computing the 3-D M-integral. 

In order to evaluation the M-integral numerically, two solutions are chosen. Solution (1) is the 

finite element solution for the problem of interest. Solution (2) is an auxiliary solution; three 

solutions are chosen in order to obtain KI

(1), KII(1), and KIII(1). These are denoted (2a), (2b), and 

(2c) and are chosen as follows 

Solution KI KII KIII

(2a) 1.0 0.0 0.0 

(2b) 0.0 1.0 0.0 

(2c) 0.0 0.0 1.0 

Figure 4. The q function used in evaluating the 3-D M-integral. 

Substituting these into eqn. 27 yields the following expressions for M
 11 

[ ( ) ( ) ( )]

[ ] ( ) () ( )

[ ] ( ) ( ) ()1

3 3

1 (1)

3 2

1

1 3

1 (1)

3 1

1

2 3

(1)

2

(1,2 ) 1

1

1 3

1

3 2

1 (1)

1 2

1 (1)

1 1

1

2 2

(1)

2

(1,2 ) 1

1

2 3

1

3 1

1 (1)

2 2

1

1 1

1 (1)

2 1

(1)

2

(1,2 ) 1

Im Im 2 Im

Im 2 Im Im

2 Im Im Im

− − − − −

− − − − −

− − − − −

= − + + + +

= − + + + +

= − + + + +

I i i i i II i i i i III i i

c

I i i i i II i i III i i i i

b

I i i II i i i i III i i i i

a

M K m N m N K m N m N K m N

M K m N m N K m N K m N m N

M K m N K m N m N K m N m N

 (29) 

and 

∫

∫

∫

∂

∂

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ − ∂

∂

+

∂

∂ =

∂

∂

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ − ∂

∂

+

∂

∂ =

∂

∂

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ − ∂

∂

+

∂

∂ =

A

j

j

c i c

ij

c

i

ij

c

A

j

j

b i b

ij

b

i

ij

b

A

j

j

a i a

ij

a

i

ij

a

dA

x

q W

x

u

x

u M

dA

x

q W

x

u

x

u M

dA

x

q W

x

u

x

u M

1

(1,2 )

1

(1)

(2 )

1

(2 )

(1,2 ) (1)

1

(1,2 )

1

(1)

(2 )

1

(2 )

(1,2 ) (1)

1

(1,2 )

1

(1)

(2 )

1

(2 )

(1,2 ) (1)

σ σ δ

σ σ δ

σ σ δ

 (30) 

These can be combined into one system of equations, 

() ( ) ( )

( ) () ( )

( )( ) ( ) ⎪

⎭

⎪

⎬

⎫

⎪

⎩

⎪

⎨

⎧

=

⎪

⎭

⎪

⎬

⎫

⎪

⎩

⎪

⎨

⎧

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎣

⎡

+ +

+ +

+ +

− − − − −

− − − − −

− − − − −

q

c

q

b

q

a

III

II

I

i i i i i i i i i i

i i i i i i i i i i

i i i i i i i i i i

M A

M A

M A

K

K

K

m N m N m N m N m N

m N m N m N m N m N

m N m N m N m N m N

(1,2 )

(1,2 )

(1,2 )

(1)

(1)

(1)

1

3 3

1

3 2

1

1 3

1

3 1

1

2 3

1

1 3

1

3 2

1

1 2

1

1 1

1

2 2

1

2 3

1

3 1

1

2 2

1

1 1

1

2 1

Im Im 2Im

Im 2Im Im

2Im Im Im

 

 (31) 

that can be solved for the stress intensity factors. The terms in the coefficient matrix are obtained 

from eqs. 7 and 9 and depend on the material properties and orientation. 

Terms in the right-hand-side vector (eqn. 30) are obtained from the finite element results and 

eqs. 1 and 2. 

For the special case where z = 0 is a plane of material symmetry, the corresponding expressions 

for the M-integral are 

( ) ( )

.

1

2

1

Im

1

2

1

Im

2 (1)

44 55 45

(1,2 )

(1)

1 2 11 16

(1)

22

1 2

1 2 11

(1,2 )

(1)

22

1 2

1 2 11

(1)

22

1 2

1 2

26

(1,2 )

III

c

I II

b

I II

a

M S S S K

M S S K S S K

M S S K S S K

= ′ ′ − ′

⎥

⎥

⎦

⎤

⎢

⎢

⎣

⎡

+ + ′ − ′ ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ = ′ − ′

⎥

⎥

⎦

⎤

⎢

⎢

⎣

⎡

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

+ ′ − ′ ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ ′ ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ + = ′ −

μ μ

μ μ

μ μ

μ μ

μ μ

μ μ

μ μ

 (32) 

These can be arranged into the system of equations 
12 

( ) ( )

⎪

⎪

⎪

⎭

⎪

⎪

⎪

⎬

⎫

⎪

⎪

⎪

⎩

⎪

⎪

⎪

⎨

⎧

=

⎪

⎭

⎪

⎬

⎫

⎪

⎩

⎪

⎨

⎧

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

⎡

′ ′ − ′

+ ′ − ′ ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ ′ − ′

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ ′ − ′ ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ ′ ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ + ′ −

q

c

q

b

q

a

III

II

I

A

M

A

M

A

M

K

K

K

S S S

S S S S

S S S S

(1,2 )

(1,2 )

(1,2 )

(1)

(1)

(1)

2

44 55 45

22 1 2 11 16

1 2

2 1 2 11

1

22

1 2

2 1 2 11

1

22

1 2

1 2

26

0 0

Im 0 1 Im

0 1 Im Im

μ μ

μ μ

μ μ

μ μ

μ μ

μ μ

μ μ

. 

 (33) 

where the terms in the right-hand-side vector are determined from eqn. 30 using eqs. 11 and 12 

to determine the stresses and displacements. 
 13 

4. M-INTEGRAL IMPLEMENTATION 

4.1 Integration Domains and q-Functions 

There are two important aspects of the implementation of the 3-D M-integral that do not follow 

directly from the formulation. These are the selection of the elements that will participate in the 

integration domain, and the related issue of selecting the virtual crack extension, the q-function. 

Our experience with 2-D implementations of the J-integral indicates that accurate stress intensity 

factors can be computed using eight triangular quarter-point elements around the crack tip. We 

obtain accurate results integrating only in these crack-tip elements. Other authors, however, 

report more accurate results when the integration is performed in "rings" of elements that are 

remote from, but completely encircle the crack tip (e.g., Banks-Sills and Sherman [1992]). In 3-

D, unlike in 2-D, the domain of integration cannot be remote from the crack-front. That is, the 

crack-front elements must be included in the integration domain, if not, it cannot be 

demonstrated that the J-integral (or M-integral) is path independent. 

Considering this experience and constraints, we have implemented the two options of evaluating 

the M-integral: either in the crack-front elements only (Fig. 5-8), or in the crack-front elements 

plus the first ring of elements surrounding the crack-front elements (Fig 9-12). The crack-front 

elements are assumed to be 15-noded wedge elements with the appropriate side nodes moved to 

the quarter-points. The elements in the next ring out from the crack front are assumed to be 20-

noded brick elements. 

Related to the selection of elements in the integration domain is the selection of an appropriate 

virtual crack extension or q-function. Banks-Sills and Sherman studied this issue in some detail 

[1988]. They determined that the most accurate results are obtained if the integral is evaluated at 

the corner node of an element with a linear variation of the q-function along the crack front in 

both sets of elements adjacent to the evaluation point. This q-variation is illustrated in the right 

of Figures 5 and 9. For the purposes of this report, this is denoted as q(1). This type of q-function 

has the attractive feature that it can be used for nodes where the crack-front meets a free surface, 

as shown in Figures 6, 7, 10, and 11, which are denoted as q(2) and q(3). 

An alternative q-function that extends only the mid-side node is shown in Figures 8 and 12 

(denoted q(4)). Banks-Sills and Sherman found that this extension was less accurate than the 

corner node extensions. However, this extension does overcome complications that arise for 

curved crack fronts. 

Curved crack fronts are usually approximated by piecewise linear segments. For anisotropic 

materials, one must consider the orientation of the material coordinate system relative to the 

crack front coordinate system. For cases q(2), q(3), and q(4) the relative orientation of the material 

coordinate system can be determined from the tangent of the participating crack front segment 

and the local normal to the crack plane. For case q(1), two crack front segments participate in the 
14 

integration domain, and, in general, their tangents will not be equal. This means that the local 

material coordinate system orientation is ambiguous. How this ambiguity is resolved in the 

current implementation is presented in Section 4.2.2. However, the accuracy of this resolution 

has not been assessed. 

 Figure 5. The domain of integration and the virtual crack extension for q(1)

using one ring of elements. 

 Figure 6. The domain of integration and the virtual crack extension for q(2)

using one ring of elements. 
 15 

 Figure 7. The domain of integration and the virtual crack extension for q(3)

using one ring of elements. 

 Figure 8. The domain of integration and the virtual crack extension for q(4)

using one ring of elements. 

 Figure 9. The domain of integration and the virtual crack extension for q(1)

using two rings of elements. 
16 

 Figure 10. The domain of integration and the virtual crack extension for q(2)

using two rings of elements. 

 

Figure 11. The domain of integration and the virtual crack extension for q(3)

using two rings of elements. 

 Figure 12. The domain of integration and the virtual crack extension for q(3)

using two rings of elements. 
 17 

A study of the relative accuracy of the combinations of element domains and q-functions has 

been performed. The results are reported in Section 5. 

4.2 Detailed Implementation Steps 

The following sub-sections detail the steps in the current implementation of the M-integral for 

anisotropic materials. 

4.2.1 Compute material constitutive matrices 

Material properties are specified by nine material properties: E1, E2, E3, ν12, ν23, ν13, G12, G23, 

and G31. The orientation of the material coordinate system relative to the Cartesian coordinate 

system of the model is specified by three Euler angles, θz, θx, and θy. 

The material compliance matrix for an orthotropic material in the local material coordinate 

system is 

[ ] .

1/

1/ 0

1/ 0 0

1/ 0 0 0

1/ / 0 0 0

1/ / / 0 0 0

ˆ

23

23

12

3

2 23 3

1 12 2 13 3

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

⎡

−

− −

=

sym G

G

G

E

E v E

E v E v E

S (34) 

This is inverted to obtain the local material stiffness matrix 

ˆ ˆ −1 C = S . (35) 

A 3x3 transformation matrix that will rotate the local material coordinate axis to the material 

coordinate axes in the Cartesian space is determined from the Euler angles as 

[ ]

⎥

⎥

⎥

⎦

⎤

−

+

⎢

⎢

⎢

⎣

⎡

−

+

− −

=

cos( )sin( )

sin( ) cos( ) cos( )sin( ) cos( )

cos( )sin( ) sin( )sin( ) cos( )

cos( )sin( ) sin( )

sin( ) cos( ) cos( )sin( )sin( ) cos( ) cos( )

cos( ) cos( ) sin( )sin( )sin( ) sin( ) cos( )

x y

z y z x y

z y z x y

x y x

z y z x y z y

z y z x y z x

t

θ θ

θ θ θ θ θ

θ θ θ θ θ

θ θ θ

θ θ θ θ θ θ θ

θ θ θ θ θ θ θ

. (36) 

Components of the 3x3 rotation matrices are used to formulate 6x6 matrices that rotate the 

material stiffness and compliance matrices. These two matrices are different, and are given by 
18 

[ ]

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

⎡

+ + +

+ + +

+ + + =

31 11 32 12 33 13 31 12 32 11 32 13 33 12 33 11 31 13

21 31 22 32 23 33 21 32 22 31 22 33 23 32 23 31 21 33

11 21 12 22 13 23 11 22 12 21 12 23 13 22 13 21 11 23

31 32 32 33 33 31

2

33

2

32

2

31

21 22 22 23 23 21

2

23

2

22

2

21

11 12 12 13 13 11

2

13

2

12

2

11

2 2 2

2 2 2

2 2 2

t t t t t t t t t t t t t t t t t t

t t t t t t t t t t t t t t t t t t

t t t t t t t t t t t t t t t t t t

t t t t t t t t t

t t t t t t t t t

t t t t t t t t t

Rs (37) 

and 

[ ]

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

⎡

+ + +

+ + +

+ + + =

31 11 32 12 33 13 31 12 32 11 32 13 33 12 33 11 31 13

21 31 22 32 23 33 21 32 22 31 22 33 23 32 23 31 21 33

11 21 12 22 13 23 11 22 12 21 12 23 13 22 13 21 11 23

31 32 32 33 33 31

2

33

2

32

2

31

21 22 22 23 23 21

2

23

2

22

2

21

11 12 12 13 13 11

2

13

2

12

2

11

2 2 2

2 2 2

2 2 2

t t t t t t t t t t t t t t t t t t

t t t t t t t t t t t t t t t t t t

t t t t t t t t t t t t t t t t t t

t t t t t t t t t

t t t t t t t t t

t t t t t t t t t

Rc . (38) 

The material matrices are then rotated into the Cartesian system by 

[ ] [ ][ ][ ]

T C Rc C Rc ˆ = and [ ] [ ][ ][ ]T S Rs S Rs ˆ = . (39) 

4.2.2 Compute the crack-front coordinate system 

For cases q(2), q(3), and q(4), the local z-axis is equal to the tangent to the crack-front segment that 

is participating in the integration. The local y-axis is parallel to the local normal to the crack 

surface. The local x-axis is perpendicular to the local y- and z-axes. 

As mentioned above, the local crack front coordinate system for q(1) is ambiguous. In the current 

implementation, the local crack front axes are determined from the coordinates of the nodes that 

are shared by the crack-front elements on either side of the evaluation location. This is illustrated 

in Figure 13. 

Figure 13. The definition of the local crack-front coordinates for case q(1), based 

on specified nodal coordinates 
 19 

4.2.3 Transform nodal coordinates, displacements, and material matrices 

The M-integral implementation requires that all input nodal coordinates, finite element 

displacements, and the material stiffness and compliance matrices be rotated to the local crackfront coordinate system. Furthermore, it requires that the nodal coordinates be translated so that 

the M-evaluation point is at the local coordinate origin. 

If the components of the normalized crack-front axes specified in the model Cartesian coordinate 

system are given by x', y', and z', then the 3x3 matrix that will transform nodal coordinates and 

displacements in the Cartesian system into the crack-front system is given by 

[ ]

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎣

⎡

′ ′ ′

′ ′ ′

′ ′ ′

′ =

z y z

y y y

x x x

x y z

x y z

x y z

t . (40) 

The corresponding transformations are 

 {x }= [ ] t′ {x − xorigin } ∗ and {u }= [t′]{u} ∗ . (41) 

The material rotation matrices are similar to those given in eqs. 37 and 38, but use components 

of the t' matrix. The corresponding transformations are 

[] [ ][ ][ ]

T C Rc C Rc = ′ ′ and [ ] [ ][ ][ ]T

s Rs S = R′ S ′ . (42) 

4.2.4 Determine the nature of the isotropy/anisotropy 

Components of the compliance matrix are examined to determine if the material is 1) isotropic, 

2) the local z = 0 plane is a plane of material symmetry, or 3) the material is generally 

anisotropic. 

For the isotropic and z = 0 symmetry cases, the compliance matrix will have distinct patterns of 

zero and non-zero entries. z = 0 symmetry will have the pattern 

[ ]

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

⎡

∗ ∗

∗ ∗

∗ ∗ ∗ ∗

∗ ∗ ∗ ∗

∗ ∗ ∗ ∗

∗ ∗ ∗ ∗

=

0 0 0 0

0 0 0 0

0 0

0 0

0 0

0 0

S . (43) 

The first check is to see if any of the upper right 6 entries of the compliance matrix are non-zero. 

If so, the material is assumed to be generally anisotropic. 
20 

An isotropic material has the pattern 

[ ]

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

⎡

∗

∗

∗

∗ ∗ ∗

∗ ∗ ∗

=

0 0 0 0 0

0 0 0 0 0

0 0 0 0 0

0 0 0

0 0 0

* * * 0 0 0

S . (44) 

The second check is to see if the 14, 24, 34, or 56 elements are zero. If so, and the following 

conditions hold 

 S11 = S22 = S33, S12 = S13 = S23, S44 = S55 = S66, and S44 = S11/2(1-S12/S11) (45) 

the material is assumed to be isotropic. Otherwise, z = 0 symmetry is assumed. 

4.2.5 Compute anisotropic parameters 

If the material has been determined to be z = 0 symmetric, then the μ's are determined by finding 

the roots to eqs. 13 and 14. 

If the material has been determined to be generally anisotropic, then the μ's, λ's, m's, and N-1 are 

determined from eqs. 4 through 9. 

4.2.6 Compute the auxiliary nodal displacements 

For all elements in the domain of integration, we determine the nodal displacements associated 

with the three auxiliary solutions, u(2a), u(2b), u(2c). The generalized plane strain expressions for 

these values are given by eqs. 2 and 12 for the general anisotropy case and the z = 0 symmetry 

case, respectively. For isotropic materials in plane strain, the expression is 

⎪

⎭

⎪

⎬

⎫

⎪

⎩

⎪

⎨

⎧

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

⎡

⎟

⎠

⎞ ⎜

⎝

⎛ ⎟ − + + ⎠⎞ ⎜

⎝

⎛ − −

⎟

⎠

⎞ ⎜

⎝

⎛ ⎟ − + ⎠⎞ ⎜

⎝

⎛ − +

=

⎪

⎭

⎪

⎬

⎫

⎪

⎩

⎪

⎨

⎧

III

II

I

z

y

x

K

K

K

r

G

u

u

u

2

0 0 2sin

0

2

1 2 sin

2

cos

2

2 2 cos

2

sin

0

2

2 2 cos

2

sin

2

1 2 sin

2

cos

2

1 2 2

2 2

θ

θ

ν

θ θ

ν

θ

θ

ν

θ θ

ν

θ

π

.(46) 

These expressions give the displacements as a function of the in-plane coordinates r and θ. A 

general approach for determining these must account for curved crack fronts where the "ends" of 

the element may not be perpendicular to the crack front. The technique for doing this in the 
 21 

current implementation is illustrated in Figure 14. The evaluation point x has element parametric 

coordinates ξ, η, and ζ. The associated point on the crack plane, x', is determined from the 

parametric coordinates ξ and η in the quadrilateral region that is the intersection of the 

evaluation ring with the crack plane. The associated point on the crack front, x'', is determined 

from the parametric coordinate ξ and the associated crack-front segment. r is the magnitude of 

the vector from x'' to x, and θ is the angle between this vector and the vector from x'' to x'. 

Figure 14. A schematic of how the coordinates r and θ are determined. 

4.2.7 Compute displacement derivatives 

For each integration point within an element, we compute 

x

u

∂

∂ (1)

, x

u a

∂

∂ (2 )

, x

u

b

∂

∂ (2 )

, and 

x

u c

∂

∂ (2 )

. 

This is done using standard finite element interpolation. The expression for the displacement 

derivatives is 

[ ] ∂u ∂x = [J ] [∂N ∂ξ ]{u} −1 (47) 

where J is the Jacobian matrix and N represents the element shape functions. The Jacobian is 

computed from the product of the shape function derivatives with respect to the parametric 

coordinates and the element's nodal coordinates: 

[J ] = [∂N ∂ξ ]{x}. (48) 

4.2.8 Compute local stresses 

For each integration point, we compute the local element stresses. For the finite element stresses, 

σ

(1), the expression is 
22 

[ ]

⎪

⎪

⎪

⎪

⎭

⎪⎪

⎪

⎪

⎬

⎫

⎪

⎪

⎪

⎪

⎩

⎪⎪

⎪

⎪

⎨

⎧

∂ ∂ + ∂ ∂

∂ ∂ + ∂ ∂

∂ ∂ + ∂ ∂

∂ ∂

∂ ∂

∂ ∂

=

⎪

⎪

⎪

⎪

⎭

⎪⎪

⎪

⎪

⎬

⎫

⎪

⎪

⎪

⎪

⎩

⎪⎪

⎪

⎪

⎨

⎧

u x u z

u z u y

u y u x

u z

u y

u x

C

z x

y z

x y

z

y

x

zx

yz

xy

z

y

x

σ

σ

σ

σ

σ

σ

. (49) 

For the auxiliary solutions (σ

(2a), σ(2b), and σ(2c)), stresses are given by eqn. 1 for the general 

anisotropy case, eqn. 11 for z = 0 symmetry, and for the case of isotropy by 

⎪

⎭

⎪

⎬

⎫

⎪

⎩

⎪

⎨

⎧

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

⎡

−

⎟

⎠

⎞ ⎜

⎝

⎛ −

⎟

⎠

⎞ ⎜

⎝

⎛ +

⎟

⎠

⎞ ⎜

⎝

⎛ ⎟ − + ⎠⎞ ⎜

⎝

⎛ −

=

⎪

⎪

⎪

⎭

⎪

⎪

⎪

⎬

⎫

⎪

⎪

⎪

⎩

⎪

⎪

⎪

⎨

⎧

III

II

I

zx

yz

xy

y

x

K

K

K

r

2

0 0 sin

2

0 0 cos

0

2

3 sin

2

1 sin

2

cos

2

3

cos

2

cos

2

sin

0

2

3

cos

2

cos

2

sin

3

3 sin

2

1 sin

2

cos

0

2

3

cos

2

2 cos

2

sin

2

3 sin

2

1 sin

2

cos

2

1

θ

θ

θ θ θ θ θ θ

θ θ θ θ θ θ

θ θ θ θ θ θ

π

σ

σ

σ

σ

σ

.(50) 

For all cases, σz comes from the plane strain condition 

 

( ) σ z = S31σ x + S32σ y + S34σ xy + S35σ yz + S36σ zx S33 . (51) 

4.2.9 Compute the contribution to the M-integral 

The M-integral is evaluated numerically using Gauss integration. It has contributions from all 

Gauss points for all elements in the integration domain. The expressions for the numerical 

integration is 

∑ ∑

∑ ∑

∑ ∑

= =

= =

= =

∂

∂

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ − ∂

∂

+

∂

∂ =

∂

∂

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ − ∂

∂

+

∂

∂ =

∂

∂

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ − ∂

∂

+

∂

∂ =

elem

elem

gp

gp j

j

c i c

ij

c

i

gp gp ij

c

elem

elem

gp

gp j

j

b i b

ij

b

i

gp gp ij

b

elem

elem

gp

gp j

j

a i a

ij

a

i

gp gp ij

a

x

q W

x

u

x

u M w J

x

q W

x

u

x

u M w J

x

q W

x

u

x

u M w J

#

1

#

1

1

(1,2 )

1

(1)

(2 )

1

(2 )

(1,2 ) (1)

#

1

#

1

1

(1,2 )

1

(1)

(2 )

1

(2 )

(1,2 ) (1)

#

1

#

1

1

(1,2 )

1

(1)

(2 )

1

(2 )

(1,2 ) (1)

σ σ δ

σ σ δ

σ σ δ

 (52) 

where wgp is the Gauss integration weight and |Jgp| is the determinate of the Jacobian matrix 

evaluated at the integration point. Summation convention is assumed with i = 1, 2, 3 and j = 1, 2, 3.
 23 

The current implementation uses 27 integration points for brick elements and 28 integration 

points for wedge (crack-front) elements. 

4.2.10 Compute the area of the virtual extension 

Aq is the area of the virtual extension along the crack front. For cases q(2), q(3), and q(4) only one 

crack-front segment participates in the integration. If this is denoted l, and, noting that the virtual 

extension at the evaluation point is 1, 

A A l q q 2

(2) (3) 1 = = and A l q 3(4) 2 = (53) 

For case q(1), two crack-front segments participate in the integration. If these are denoted l1 and 

l2, then 

( ) 2 1 2

(1) 1 A l l q = + . (54) 

4.2.11 Compute the stress intensity factors 

Finally, the stress intensity factors are computed. Eqn. 31 is used for general anisotropy, eqn 33 

for z = 0 symmetry, and for isotropy, the expression is 

( )

( )

( ) ⎪

⎭

⎪

⎬

⎫

⎪

⎩

⎪

⎨

⎧

=

⎪

⎭

⎪

⎬

⎫

⎪

⎩

⎪

⎨

⎧

⎥

⎥

⎥

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎢

⎢

⎢

⎣

⎡

+

−

−

q

c

q

b

q

a

III

II

I

M A

M A

M A

K

K

K

E

E

E

(1,2 )

(1,2 )

(1,2 )

2

2

2 1 0 0

0 2 1 0

0 0 2 1

ν

ν

ν

. (55) 
24 

5. PRESCRIBED DISPLACEMENT TESTS 

An initial series of tests were performed to verify the implementation. In these tests, no finite 

element analysis was performed. Rather, nodal displacements for all elements participating in the 

M-integral evaluation were determined from eqs. 2, 12, and 46, assuming unit stress intensity 

values. The M-integral is evaluated, and the computed stress intensity values are compared to the 

prescribed values. 

The purpose of these tests is two-fold. First, as verification; if the M-integral evaluation is 

encoded correctly, then the computed stress intensity factor should be close to the prescribed 

value. Second, the absolute difference between the prescribed and computed values give a 

quantitative assessment of the "intrinsic" error in the computations and an indication of the 

relative performance of differing q functions and numbers of element rings. By intrinsic error, 

we mean those errors that arise due to the approximate numerical techniques used to evaluate the 

integrals. 

Table 1 shows the computed stress intensity factors for isotropy, z = 0 symmetry, and general 

anisotropy. For each, all three fracture modes are examined for both 1 and 2 element rings 

around the crack front for four different q functions. 

Of all the tests, the largest error was 4.7%. This was for general anisotropy, mode I, one ring of 

elements, and q(4). In general, all the q(4) tests with one ring of elements performed poorly. If 

these are discounted, then the largest error for the one-ring tests was 0.89%. Among the tworing tests the largest error was 0.46%. 

Based on these tests, the following conclusions are drawn: 

1. All combinations of element rings and q functions produce accurate results with the 

exception of q(4) with one ring of elements. This configuration should be avoided. 

2. In general the two-ring configurations are slightly more accurate than the one-ring 

configurations, but the difference is not significant. 
 25 

6. THICK PLATE ANALYSES 

A series of finite element analyses were performed to verify the performance of the M-integral 

calculations. The models were essentially very thick center-cracked plates. The objective was to 

make the plates thick enough that plane-strain conditions would be approximated near the center 

of the plate. 

The geometry of the model is shown in Figure 15. The plate is simply supported as shown in 

Figure 16a. Two load cases, mode I (Figure 16b) and mode II (Figure 16c) were analyzed. All 

length and force units are normalized with respect to those used for the material properties. 

Table 1. Computed stress intensity factors for prescribed nodal displacements. 

 prescribed KI prescribed KII prescribed KIII

rings q(n) KI KII KIII KI KII KIII KI KII KIII

 Isotropic 

1 1.0089 0.0 0.0 0.0 1.0062 0.0 0.0 0.0 1.0070 

2 1.0089 0.0 0.0 0.0 1.0062 0.0 0.0 0.0 1.0070 

3 1.0089 0.0 0.0 0.0 1.0062 0.0 0.0 0.0 1.0070 

1 

4 0.9591 0.0 0.0 0.0 0.9647 0.0 0.0 0.0 0.9562 

 

1 0.9984 0.0 0.0 0.0 0.9987 0.0 0.0 0.0 0.9992 

2 0.9984 0.0 0.0 0.0 0.9987 0.0 0.0 0.0 0.9992 

3 0.9984 0.0 0.0 0.0 0.9987 0.0 0.0 0.0 0.9992 

2 

4 0.9984 0.0 0.0 0.0 0.9987 0.0 0.0 0.0 0.9992 

 Z symmetry 

1 1.0080 0.0 0.0 0.0 1.0065 0.0 0.0 0.0 1.0069 

2 1.0080 0.0 0.0 0.0 1.0065 0.0 0.0 0.0 1.0069 

3 1.0080 0.0 0.0 0.0 1.0065 0.0 0.0 0.0 1.0069 

1 

4 0.9534 0.0 0.0 0.0 0.9624 0.0 0.0 0.0 0.9543 

 

1 1.0034 0.0 0.0 0.0 1.0009 0.0 0.0 0.0 1.0001 

2 1.0034 0.0 0.0 0.0 1.0009 0.0 0.0 0.0 1.0001 

3 1.0034 0.0 0.0 0.0 1.0009 0.0 0.0 0.0 1.0001 

2 

4 1.0034 0.0 0.0 0.0 1.0009 0.0 0.0 0.0 1.0001 

 general anisotropy 

1 1.0078 0.0 -0.0007 0.0 1.0064 0.0 -0.0013 0.0 1.0073 

2 1.0078 0.0 -0.0007 0.0 1.0064 0.0 -0.0013 0.0 1.0073 

3 1.0078 0.0 -0.0007 0.0 1.0064 0.0 -0.0013 0.0 1.0073 

1 

4 0.9529 0.0 -0.0024 0.0 0.9702 0.0 -0.0062 0.0 0.9591 

 

1 1.0046 0.0 0.0018 0.0 0.9968 0.0 0.0045 0.0 0.9980 

2 1.0046 0.0 0.0018 0.0 0.9968 0.0 0.0045 0.0 0.9980 

3 1.0046 0.0 0.0018 0.0 0.9968 0.0 0.0045 0.0 0.9980 

2 

4 1.0046 0.0 0.0018 0.0 0.9968 0.0 0.0045 0.0 0.9980 
26 

Figure 15. The geometry used for the thick plate analysis. 

Figure 16. The simply supported boundary conditions (a), and the mode I (b) 

and II (c) load cases used in the analyses. 

Three sets of material properties were investigated 

Isotropic: 

E = 1000.0, ν = 0.25 

Cubic material oriented with z = 0 symmetry: 

E1 = E2 = E3 = G12 = G23 = G13 = 1000.0 

ν12 = ν23 = ν13 = 0.25 

θz = 45.0°, θx = θy = 0.0°
 27 

Cubic material oriented with a general orientation: 

E1 = E2 = E3 = G12 = G23 = G13 = 1000.0 

ν12 = ν23 = ν13 = 0.25 

θz = θx = θy = 45.0°

The model and meshes were generated with FRANC3D. The finite element analysis was 

performed with ANSYS. The resulting nodal displacements were read back into FRANC3D 

where the M-integrals were evaluated to determine the stress intensity factors. 

The same mesh was used for all analyses; it is shown in Figure 17, with a detail of the crack 

region shown in Figure 18. The mesh is predominantly 10-noded tetrahedral elements. 15-noded 

quarter point wedge elements are used around the crack front. One ring of 20-noded brick 

elements is used around the crack-front elements. 13-noded pyramid elements are used to 

transition from the brick elements to the tetrahedral elements in the bulk of the mesh. 

Because of the simple loading and the simple supports, the well-known K solution for a centercracked plate is valid for all three material behaviors. The relative errors between the theoretical 

solution and the numerical solutions are shown graphically in Figure 19. It is expected that there 

will be a relatively large error near the surfaces of the plate where plane strain conditions are not 

satisfied. It can be seen that near the center of the plate the M-integral solutions are very 

accurate. 

Figure 17. The mesh used for the thick plate verification analyses. 
28 

Figure 18. A detail of the thick plate mesh showing the crack region. 

-5

-3

-1

1

3

5

0 0.2 0.4 0.6 0.8 1

normalized distance along crack front

% error in stress intensity factor

isotropic Mode I

z=0 sym Mode I

general Mode I

isotropic Mode II

z=0 sym Mode II

general Mode II

Figure 19. The relative errors between the theoretical and numerical solutions for the stress 

intensity factors for a thick center-cracked plate. 
 29 

7. BRAZILIAN DISK ANALYSIS 

A large number of analyses were performed to determine a stress intensity solution for a 

Brazilian disk test specimen. The specimen geometry is shown in Figure 20. The material 

properties considered are those for a typical single crystal nickel alloy used in an engine hot 

section. An isotropic material was considered as well. In testing, it is observed that the two ends 

of the crack do not usually grow at the same rate, therefore, a range of a1 and a2 lengths are 

considered. 

W

t

y

x

P

θ

a1

a2

W

t

y

x

P

θ

a1

a2

Figure 20. The geometry of the Brazilian disk test specimen. 

The single crystal material has cubic symmetry. This means that E = E001 = E010 = E100, ν = 

ν001/010 = ν001/100 = ν010/100, and G = G001/010 = G001/100 = G010/100. The properties used here are E/G

= 0.952 and ν = 0.4. The properties used for the isotropic material are E/G = 2.8 and ν = 0.4. 

Two different material orientations were considered for the anisotropic material, <111/110> and 

<111/112> where <y/x> defines the orientation of the local crack front axes relative to the 

material axes. These orientations are shown graphically in Figure 21, where the red axes indicate 

the local crack front-axes. The corresponding Euler angles are θz = 50.7685, θx = −24.0948, and 

θy = 26.5650, for the <111/110> case, and θz = −35.2644, θx = −45.0, and θy = −90.0, for the 

<111/112> case. 
30 

Figure 21. A graphical illustration of the cubic material orientation relative to the local crackfront axes. 

Twenty-eight crack configurations were considered. These are shown in Table 2. Eleven load 

orientations, θ's, were considered. The isotropic and <111/110> material configurations give 

symmetric results for negative and positive θ angles. For these cases, 0°, 4°, 8°, 12°, 16° and 20°

were analyzed. For the <111/112> material orientation, −4°, −8°, −12°, −16°, and −20° were 

analyzed in addition to the positive angles. 

Table 2. Analyzed crack configuration. 

2a1/W 2a2/W

 

0.2 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8 

0.3 0.3, 0.4, 0.5, 0.6, 0.7, 0.8 

0.4 0.4, 0.5, 0.6, 0.7, 0.8 

0.5 0.5, 0.6, 0.7, 0.8 

0.6 0.6, 0.7, 0.8 

0.7 0.7, 0.8 

0.8 0.8 

 

A typical mesh used in these analyses is shown in Figure 22 (with greatly magnified 

displacements). The meshes were two elements thick and periodic boundary conditions were 

applied to the front and back surfaces to enforce plane strain behavior. The meshes were 

predominantly 15-noded wedge elements, with quarter-point wedge elements used at the crack 

fronts and 20-noded brick elements used in the first ring of elements surrounding the crack-front 

elements. The M-integral was evaluated using q(1). 
 31 

Figure 22. A typical mesh used in the disk analyses (with greatly 

magnified displacements). 

Typical analysis results for are shown graphically in Figure 23. Complete analysis results are 

given in graphical and tabular form in the appendices. 
32 

0.000

0.500

1.000

1.500

2.000

2.500

3.000

3.500

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

0.000

0.500

1.000

1.500

2.000

2.500

3.000

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

Mode I: Crack Tip a2

Mode I: Crack Tip a1

2a2/W

2a2/W

0.000

0.500

1.000

1.500

2.000

2.500

3.000

3.500

4.000

4.500

0.2 0.4 0.6 0.8

2a1/W

FII

0 .2 0.3 0.4 0 .5 0.6

0 .7 0.8

Mode II: Crack Tip a1

2a2/W

0.000

0.500

1.000

1.500

2.000

2.500

3.000

0.2 0.4 0.6 0.8

2a1/W

FII

0.2 0.3 0 .4 0.5 0 .6

0.7 0.8

Mode II: Crack Tip a2

2a2/W

Figure 23. Typical Brazilian disk analysis results, where 

( ) 1 2 P 2 a a

K Wt F n

n + = π

π , with n = I, II. 
 33 

8. SUMMARY 

This report describes the successful effort to develop and verify a numerical technique for 

extracting stress intensity factors from finite element results for anisotropic materials. The report 

can be summarized as follows: 

• Expressions for the stress and displacement fields near a crack front in anisotropic material 

are presented. 

• The expressions for the M-integral that can be used to determine the mode-wise stress 

intensity factors in such materials are developed. 

• The details of an implementation of the anisotropic M-integral are presented. 

• The M-integral is shown to give consistent and accurate results for test cases where nodal 

displacements are prescribed to be the theoretical near crack-front values. 

• The M-integral is shown to give accurate results for a simple model with a known stress 

intensity factor solution. 

• The new capability is used to develop the K-solution for crack growth in an anisotropic 

Brazilian disk specimen. 
34 

REFERENCES 

Atluri S.N., Kobayashi A.S., Nakagaki M., "A finite-element program for fracture mechanics 

analysis of composite material." In: Sendeckyj G.P., editor. Fracture Mechanics of Composites, 

ASTM-STP 593. American Society for Testing and Materials, 1975. p.~86-98. 

Banks-Sill, L. and Sherman, D., "On quarter-point three-dimensional finite elements in linear 

elastic fracture mechanics," Int J Fracture 1989; 41:177-196. 

Banks-Sill, L. and Sherman, D., "On the computation of stress intensity factors for threedimensional geometries by means of the stiffness derivative and J-integral methods," Int J 

Fracture 1992; 53:1-20. 

Boone T.J., Wawrzynek P.A., and Ingraffea A.R., "Finite element modeling of fracture 

propagation in orthotropic materials," Engng Fract Mech. 1987; 26:185-201. 

Hoenig A., "Near-tip behavior of a crack in a plane anisotropic elastic body," Engng Fract Mech. 

1982; 16(3):393-403. 

Lekhnitskii S.G., Theory of Elasticity of an Anisotropic Body. San Francisco: Holden-Day, 1950, 

in Russian, 1963, in English, translated by Fern P. 

Rice J.R., "A path independent integral and the approximate analysis of strain concentration by 

notches and cracks," J Appl Mech. 1968; 35:379-386. 

Sih G.C. Paris P.C. Irwin G.R., "On cracks in rectilinearly anisotropic bodies," Int J Fract Mech.

1965; 1:189--302. 

Tohgo K., Wang A.S.D., and Chou T., "A criterion for splitting crack initiation in unidirectional 

fiber-reinforced composites," J Compos Ma. 1993; 27:1054-1076. 

Yau J.F., Wang S.S., Corten H.T., "A mixed-mode crack analysis of isotropic solids using 

conservation laws of elasticity," J Appl Mech. 1980; 47:335-341. 
 35 

APPENDIX A: BRAZILIAN DISK ANALYSIS RESULTS 

This appendix contains stress intensity factor results for the Brazilian disk specimens. 

Analyses were performed for a cubic material with properties E/G = 0.952 and ν = 0.4 and for an 

isotropic material with properties E/G = 2.8 and ν = 0.4. Two material orientations relative to the 

crack-font were considered for the cubic material, <111/110> and <111/112> where <y/x> 

defines the orientation of the local crack front axes relative to the material axes. 

For the isotropic material and the cubic material in the <111/110> orientation, θ angles of 0°, 4°, 

8°, 12°, 16° and 20° were analyzed. Due to material symmetries, these results are symmetric 

about θ = 0°. The cubic material in the <111/112> orientation is not symmetric about θ = 0°. For 

that case, −4°, −8°, −12°, −16°, and −20° were analyzed in addition to the positive angles. 

A total of 644 analyses were performed. 

The first part of this appendix contains the stress intensity factor values in tabular form. Each 

page contains the data for one load orientation for one material type/orientation. 

The second part of this appendix contains plots of normalized stress intensity factors for Modes I 

and II. The normalized stress intensity factors are defined as 

( ) 1 2 P 2 a a

K Wt F n

n + = π

π where n = I, II. 

 
36 

<111/110>, θ = 20°

Mode I: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.128 0.127 0.107 0.058 -0.035 -0.200 -0.495 

 tip 2 0.128 0.147 0.163 0.174 0.175 0.163 0.124 

0.3 tip 1 0.147 0.127 0.078 -0.015 -0.182 -0.480 

 tip 2 0.147 0.163 0.176 0.180 0.170 0.135 

0.4 tip 1 0.145 0.097 0.003 -0.164 -0.466 

 tip 2 0.145 0.158 0.164 0.156 0.125 

0.5 tip 1 0.111 0.018 -0.151 -0.455 

 tip 2 0.111 0.118 0.112 0.083 

0.6 tip 1 0.025 -0.144 -0.452 

 tip 2 0.025 0.020 -0.007 

0.7 tip 1 -0.150 -0.460 

 tip 2 -0.150 -0.177 

0.8 tip 1 -0.488 

 tip 2 -0.488 

Mode II: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 -0.520 -0.612 -0.717 -0.834 -0.967 -1.108 -1.256 

 tip 2 -0.519 -0.593 -0.669 -0.749 -0.837 -0.936 -1.050 

0.3 tip 1 -0.685 -0.790 -0.911 -1.050 -1.199 -1.361 

 tip 2 -0.683 -0.759 -0.841 -0.932 -1.034 -1.153 

0.4 tip 1 -0.868 -0.992 -1.138 -1.298 -1.475 

 tip 2 -0.866 -0.951 -1.046 -1.153 -1.279 

0.5 tip 1 -1.082 -1.235 -1.405 -1.600 

 tip 2 -1.080 -1.180 -1.295 -1.430 

0.6 tip 1 -1.344 -1.527 -1.743 

 tip 2 -1.342 -1.467 -1.616 

0.7 tip 1 -1.668 -1.908 

 tip 2 -1.665 -1.834 

0.8 tip 1 -2.108 

 tip 2 -2.105 

Mode III: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 -0.020 -0.031 -0.046 -0.065 -0.081 -0.083 -0.053 

 tip 2 0.020 0.023 0.027 0.032 0.037 0.044 0.051 

0.3 tip 1 -0.034 -0.049 -0.068 -0.084 -0.086 -0.056 

 tip 2 0.034 0.037 0.042 0.047 0.053 0.060 

0.4 tip 1 -0.053 -0.071 -0.088 -0.090 -0.061 

 tip 2 0.053 0.057 0.062 0.068 0.075 

0.5 tip 1 -0.076 -0.093 -0.095 -0.067 

 tip 2 0.076 0.082 0.088 0.094 

0.6 tip 1 -0.099 -0.103 -0.076 

 tip 2 0.099 0.106 0.113 

0.7 tip 1 -0.111 -0.086 

 tip 2 0.111 0.120 

0.8 tip 1 -0.097 

 tip 2 0.097 
 37 

<111/110>, θ = 16°

Mode I: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.186 0.199 0.196 0.168 0.098 -0.042 -0.312 

 tip 2 0.186 0.214 0.241 0.263 0.279 0.284 0.264 

0.3 tip 1 0.228 0.226 0.199 0.131 -0.009 -0.280 

 tip 2 0.228 0.255 0.280 0.299 0.306 0.291 

0.4 tip 1 0.255 0.230 0.162 0.023 -0.248 

 tip 2 0.255 0.282 0.302 0.313 0.302 

0.5 tip 1 0.258 0.192 0.054 -0.218 

 tip 2 0.258 0.281 0.293 0.286 

0.6 tip 1 0.216 0.079 -0.193 

 tip 2 0.216 0.231 0.227 

0.7 tip 1 0.095 -0.177 

 tip 2 0.095 0.092 

0.8 tip 1 -0.180 

 tip 2 -0.179 

Mode II: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 -0.448 -0.535 -0.640 -0.769 -0.928 -1.110 -1.312 

 tip 2 -0.446 -0.511 -0.578 -0.650 -0.734 -0.831 -0.947 

0.3 tip 1 -0.598 -0.704 -0.836 -1.001 -1.191 -1.407 

 tip 2 -0.595 -0.663 -0.737 -0.822 -0.921 -1.041 

0.4 tip 1 -0.773 -0.908 -1.079 -1.279 -1.510 

 tip 2 -0.770 -0.847 -0.935 -1.039 -1.166 

0.5 tip 1 -0.989 -1.167 -1.377 -1.625 

 tip 2 -0.987 -1.080 -1.190 -1.325 

0.6 tip 1 -1.268 -1.490 -1.758 

 tip 2 -1.265 -1.385 -1.533 

0.7 tip 1 -1.625 -1.918 

 tip 2 -1.623 -1.790 

0.8 tip 1 -2.115 

 tip 2 -2.114 

Mode III: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.010 0.003 -0.012 -0.034 -0.063 -0.088 -0.084 

 tip 2 -0.010 -0.011 -0.011 -0.009 -0.006 0.000 0.007 

0.3 tip 1 0.004 -0.010 -0.033 -0.062 -0.087 -0.083 

 tip 2 -0.004 -0.004 -0.002 0.000 0.005 0.011 

0.4 tip 1 -0.011 -0.033 -0.062 -0.087 -0.083 

 tip 2 0.011 0.012 0.014 0.018 0.024 

0.5 tip 1 -0.034 -0.063 -0.088 -0.084 

 tip 2 0.034 0.037 0.040 0.046 

0.6 tip 1 -0.066 -0.092 -0.089 

 tip 2 0.066 0.071 0.076 

0.7 tip 1 -0.097 -0.096 

 tip 2 0.097 0.104 

0.8 tip 1 -0.105 

 tip 2 0.106 
38 

<111/110>, θ = 12°

Mode I: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.237 0.264 0.282 0.281 0.247 0.142 -0.098 

 tip 2 0.237 0.274 0.310 0.345 0.376 0.399 0.403 

0.3 tip 1 0.302 0.321 0.323 0.290 0.188 -0.050 

 tip 2 0.302 0.340 0.377 0.410 0.437 0.446 

0.4 tip 1 0.361 0.365 0.335 0.236 0.001 

 tip 2 0.361 0.400 0.436 0.466 0.480 

0.5 tip 1 0.406 0.379 0.284 0.052 

 tip 2 0.406 0.446 0.479 0.496 

0.6 tip 1 0.421 0.329 0.099 

 tip 2 0.422 0.458 0.479 

0.7 tip 1 0.367 0.140 

 tip 2 0.367 0.391 

0.8 tip 1 0.166 

 tip 2 0.166 

Mode II: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 -0.356 -0.431 -0.527 -0.656 -0.834 -1.067 -1.353 

 tip 2 -0.354 -0.406 -0.460 -0.521 -0.594 -0.684 -0.800 

0.3 tip 1 -0.481 -0.579 -0.710 -0.893 -1.134 -1.433 

 tip 2 -0.478 -0.533 -0.595 -0.669 -0.760 -0.879 

0.4 tip 1 -0.634 -0.768 -0.957 -1.206 -1.521 

 tip 2 -0.631 -0.695 -0.772 -0.866 -0.990 

0.5 tip 1 -0.836 -1.030 -1.289 -1.619 

 tip 2 -0.833 -0.913 -1.014 -1.144 

0.6 tip 1 -1.117 -1.387 -1.736 

 tip 2 -1.114 -1.223 -1.365 

0.7 tip 1 -1.509 -1.881 

 tip 2 -1.507 -1.666 

0.8 tip 1 -2.070 

 tip 2 -2.067 

Mode III: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.041 0.039 0.031 0.012 -0.021 -0.068 -0.108 

 tip 2 -0.041 -0.046 -0.049 -0.052 -0.053 -0.050 -0.044 

0.3 tip 1 0.043 0.035 0.017 -0.016 -0.063 -0.102 

 tip 2 -0.043 -0.048 -0.050 -0.052 -0.051 -0.046 

0.4 tip 1 0.039 0.021 -0.011 -0.058 -0.097 

 tip 2 -0.039 -0.043 -0.044 -0.044 -0.040 

0.5 tip 1 0.024 -0.008 -0.054 -0.093 

 tip 2 -0.024 -0.026 -0.027 -0.023 

0.6 tip 1 -0.006 -0.052 -0.092 

 tip 2 0.006 0.006 0.009 

0.7 tip 1 -0.052 -0.093 

 tip 2 0.053 0.057 

0.8 tip 1 -0.098 

 tip 2 0.098 
 39 

<111/110>, θ = 8°

Mode I: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.278 0.317 0.355 0.386 0.399 0.364 0.192 

 tip 2 0.278 0.321 0.365 0.410 0.456 0.500 0.534 

0.3 tip 1 0.362 0.401 0.435 0.453 0.423 0.257 

 tip 2 0.362 0.408 0.455 0.504 0.551 0.590 

0.4 tip 1 0.449 0.486 0.508 0.483 0.324 

 tip 2 0.449 0.500 0.551 0.602 0.645 

0.5 tip 1 0.540 0.566 0.546 0.394 

 tip 2 0.540 0.595 0.649 0.697 

0.6 tip 1 0.625 0.610 0.464 

 tip 2 0.625 0.684 0.735 

0.7 tip 1 0.672 0.533 

 tip 2 0.674 0.730 

0.8 tip 1 0.593 

 tip 2 0.593 

Mode II: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 -0.247 -0.302 -0.377 -0.484 -0.650 -0.910 -1.314 

 tip 2 -0.245 -0.281 -0.320 -0.364 -0.419 -0.493 -0.598 

0.3 tip 1 -0.337 -0.413 -0.522 -0.691 -0.958 -1.374 

 tip 2 -0.335 -0.373 -0.418 -0.474 -0.548 -0.654 

0.4 tip 1 -0.452 -0.563 -0.737 -1.011 -1.439 

 tip 2 -0.449 -0.495 -0.552 -0.628 -0.738 

0.5 tip 1 -0.612 -0.790 -1.071 -1.512 

 tip 2 -0.608 -0.668 -0.749 -0.863 

0.6 tip 1 -0.855 -1.145 -1.602 

 tip 2 -0.851 -0.937 -1.061 

0.7 tip 1 -1.242 -1.719 

 tip 2 -1.237 -1.374 

0.8 tip 1 -1.881 

 tip 2 -1.880 

Mode III: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.067 0.072 0.072 0.065 0.045 -0.002 -0.089 

 tip 2 -0.067 -0.075 -0.082 -0.089 -0.095 -0.099 -0.098 

0.3 tip 1 0.079 0.080 0.074 0.054 0.007 -0.079 

 tip 2 -0.079 -0.087 -0.094 -0.101 -0.106 -0.106 

0.4 tip 1 0.088 0.082 0.062 0.017 -0.068 

 tip 2 -0.088 -0.096 -0.103 -0.109 -0.110 

0.5 tip 1 0.090 0.071 0.027 -0.056 

 tip 2 -0.090 -0.098 -0.104 -0.107 

0.6 tip 1 0.080 0.035 -0.047 

 tip 2 -0.080 -0.087 -0.090 

0.7 tip 1 0.043 -0.040 

 tip 2 -0.044 -0.047 

0.8 tip 1 -0.036 

 tip 2 0.037 
40 

<111/110>, θ = 4°

Mode I: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.304 0.352 0.404 0.460 0.522 0.583 0.595 

 tip 2 0.304 0.351 0.401 0.453 0.510 0.573 0.642 

0.3 tip 1 0.401 0.455 0.515 0.582 0.650 0.674 

 tip 2 0.401 0.452 0.507 0.566 0.632 0.706 

0.4 tip 1 0.509 0.573 0.645 0.720 0.754 

 tip 2 0.509 0.567 0.629 0.699 0.776 

0.5 tip 1 0.634 0.711 0.793 0.838 

 tip 2 0.634 0.701 0.774 0.855 

0.6 tip 1 0.783 0.871 0.927 

 tip 2 0.783 0.861 0.947 

0.7 tip 1 0.955 1.019 

 tip 2 0.955 1.046 

0.8 tip 1 1.118 

 tip 2 1.120 

Mode II: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 -0.127 -0.157 -0.198 -0.259 -0.362 -0.550 -0.956 

 tip 2 -0.125 -0.144 -0.164 -0.187 -0.217 -0.260 -0.329 

0.3 tip 1 -0.175 -0.216 -0.279 -0.383 -0.575 -0.989 

 tip 2 -0.172 -0.192 -0.215 -0.246 -0.288 -0.358 

0.4 tip 1 -0.236 -0.300 -0.407 -0.602 -1.025 

 tip 2 -0.233 -0.257 -0.288 -0.331 -0.403 

0.5 tip 1 -0.325 -0.435 -0.635 -1.065 

 tip 2 -0.321 -0.353 -0.399 -0.473 

0.6 tip 1 -0.470 -0.675 -1.113 

 tip 2 -0.465 -0.514 -0.592 

0.7 tip 1 -0.729 -1.179 

 tip 2 -0.723 -0.810 

0.8 tip 1 -1.283 

 tip 2 -1.274 

Mode III: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.084 0.094 0.103 0.110 0.113 0.103 0.048 

 tip 2 -0.084 -0.095 -0.105 -0.115 -0.126 -0.137 -0.148 

0.3 tip 1 0.104 0.114 0.121 0.125 0.116 0.064 

 tip 2 -0.104 -0.115 -0.125 -0.137 -0.149 -0.161 

0.4 tip 1 0.124 0.132 0.137 0.130 0.079 

 tip 2 -0.124 -0.135 -0.147 -0.160 -0.172 

0.5 tip 1 0.143 0.149 0.144 0.097 

 tip 2 -0.143 -0.156 -0.170 -0.184 

0.6 tip 1 0.163 0.159 0.113 

 tip 2 -0.163 -0.177 -0.192 

0.7 tip 1 0.176 0.133 

 tip 2 -0.176 -0.193 

0.8 tip 1 0.151 

 tip 2 -0.153 
 41 

<111/110>, θ = 0°

Mode I: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.312 0.364 0.421 0.488 0.570 0.681 0.852 

 tip 2 0.312 0.362 0.413 0.468 0.529 0.600 0.688 

0.3 tip 1 0.415 0.474 0.544 0.632 0.750 0.936 

 tip 2 0.415 0.468 0.525 0.589 0.662 0.754 

0.4 tip 1 0.530 0.604 0.698 0.824 1.022 

 tip 2 0.530 0.590 0.657 0.734 0.829 

0.5 tip 1 0.668 0.767 0.901 1.112 

 tip 2 0.668 0.739 0.820 0.919 

0.6 tip 1 0.844 0.986 1.208 

 tip 2 0.844 0.930 1.033 

0.7 tip 1 1.078 1.310 

 tip 2 1.077 1.185 

0.8 tip 1 1.425 

 tip 2 1.426 

Mode II: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 -0.001 -0.002 -0.002 -0.003 -0.004 -0.004 -0.012 

 tip 2 0.001 0.000 0.000 0.000 -0.001 0.000 -0.001 

0.3 tip 1 -0.002 -0.002 -0.003 -0.004 -0.004 -0.010 

 tip 2 0.001 0.001 0.001 0.000 0.000 -0.001 

0.4 tip 1 -0.002 -0.003 -0.003 -0.004 -0.010 

 tip 2 0.001 0.001 0.001 0.001 0.000 

0.5 tip 1 -0.002 -0.003 -0.004 -0.009 

 tip 2 0.002 0.001 0.002 -0.001 

0.6 tip 1 -0.004 -0.005 -0.012 

 tip 2 0.001 0.001 0.001 

0.7 tip 1 -0.005 -0.009 

 tip 2 0.001 0.001 

0.8 tip 1 -0.009 

 tip 2 0.001 

Mode III: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.090 0.103 0.115 0.128 0.143 0.162 0.193 

 tip 2 -0.090 -0.102 -0.113 -0.125 -0.137 -0.152 -0.171 

0.3 tip 1 0.113 0.126 0.139 0.156 0.177 0.211 

 tip 2 -0.113 -0.125 -0.136 -0.150 -0.166 -0.185 

0.4 tip 1 0.137 0.151 0.168 0.191 0.228 

 tip 2 -0.137 -0.150 -0.163 -0.180 -0.200 

0.5 tip 1 0.165 0.183 0.207 0.247 

 tip 2 -0.165 -0.179 -0.196 -0.218 

0.6 tip 1 0.199 0.225 0.268 

 tip 2 -0.199 -0.216 -0.239 

0.7 tip 1 0.247 0.293 

 tip 2 -0.247 -0.271 

0.8 tip 1 0.320 

 tip 2 -0.320 
42 

<111/112>, θ = 20°

Mode I: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.044 0.029 -0.007 -0.072 -0.180 -0.349 -0.610 

 tip 2 0.044 0.050 0.053 0.048 0.035 0.006 -0.047 

0.3 tip 1 0.035 -0.001 -0.068 -0.178 -0.352 -0.619 

 tip 2 0.035 0.037 0.034 0.021 -0.007 -0.059 

0.4 tip 1 0.001 -0.067 -0.180 -0.358 -0.633 

 tip 2 0.001 -0.003 -0.017 -0.046 -0.098 

0.5 tip 1 -0.072 -0.188 -0.371 -0.654 

 tip 2 -0.072 -0.088 -0.117 -0.171 

0.6 tip 1 -0.207 -0.394 -0.687 

 tip 2 -0.207 -0.240 -0.297 

0.7 tip 1 -0.433 -0.735 

 tip 2 -0.434 -0.497 

0.8 tip 1 -0.807 

 tip 2 -0.807 

Mode II: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 -0.548 -0.648 -0.759 -0.882 -1.014 -1.142 -1.260 

 tip 2 -0.548 -0.628 -0.709 -0.795 -0.889 -0.993 -1.109 

0.3 tip 1 -0.726 -0.839 -0.964 -1.102 -1.240 -1.373 

 tip 2 -0.726 -0.808 -0.896 -0.992 -1.099 -1.221 

0.4 tip 1 -0.923 -1.052 -1.197 -1.345 -1.495 

 tip 2 -0.923 -1.014 -1.115 -1.227 -1.358 

0.5 tip 1 -1.148 -1.301 -1.461 -1.629 

 tip 2 -1.149 -1.255 -1.376 -1.517 

0.6 tip 1 -1.418 -1.590 -1.781 

 tip 2 -1.418 -1.550 -1.705 

0.7 tip 1 -1.739 -1.956 

 tip 2 -1.739 -1.916 

0.8 tip 1 -2.166 

 tip 2 -2.166 

Mode III: all zero 
 43 

<111/112>, θ = 16°

Mode I: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.123 0.124 0.107 0.061 -0.031 -0.191 -0.460 

 tip 2 0.123 0.142 0.157 0.168 0.171 0.160 0.125 

0.3 tip 1 0.143 0.126 0.080 -0.012 -0.175 -0.449 

 tip 2 0.143 0.160 0.171 0.175 0.166 0.133 

0.4 tip 1 0.143 0.097 0.004 -0.161 -0.440 

 tip 2 0.143 0.155 0.160 0.151 0.119 

0.5 tip 1 0.109 0.015 -0.153 -0.436 

 tip 2 0.109 0.113 0.104 0.072 

0.6 tip 1 0.017 -0.154 -0.443 

 tip 2 0.017 0.006 -0.028 

0.7 tip 1 -0.169 -0.465 

 tip 2 -0.169 -0.208 

0.8 tip 1 -0.512 

 tip 2 -0.511 

Mode II: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 -0.453 -0.548 -0.664 -0.806 -0.976 -1.158 -1.339 

 tip 2 -0.453 -0.520 -0.590 -0.668 -0.756 -0.858 -0.979 

0.3 tip 1 -0.614 -0.731 -0.876 -1.051 -1.243 -1.439 

 tip 2 -0.614 -0.684 -0.762 -0.852 -0.957 -1.082 

0.4 tip 1 -0.804 -0.951 -1.134 -1.336 -1.548 

 tip 2 -0.804 -0.885 -0.978 -1.086 -1.218 

0.5 tip 1 -1.037 -1.226 -1.439 -1.669 

 tip 2 -1.038 -1.136 -1.251 -1.392 

0.6 tip 1 -1.333 -1.558 -1.810 

 tip 2 -1.334 -1.459 -1.613 

0.7 tip 1 -1.700 -1.978 

 tip 2 -1.701 -1.875 

0.8 tip 1 -2.185 

 tip 2 -2.186 

Mode III: all zero 
44 

<111/112>, θ = 12°

Mode I: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.196 0.215 0.223 0.207 0.149 0.012 -0.261 

 tip 2 0.196 0.226 0.255 0.282 0.304 0.314 0.303 

0.3 tip 1 0.247 0.255 0.240 0.183 0.046 -0.228 

 tip 2 0.247 0.277 0.305 0.329 0.341 0.333 

0.4 tip 1 0.286 0.273 0.217 0.081 -0.195 

 tip 2 0.286 0.316 0.341 0.355 0.349 

0.5 tip 1 0.304 0.248 0.112 -0.164 

 tip 2 0.304 0.330 0.345 0.340 

0.6 tip 1 0.274 0.137 -0.142 

 tip 2 0.274 0.289 0.283 

0.7 tip 1 0.151 -0.133 

 tip 2 0.150 0.142 

0.8 tip 1 -0.146 

 tip 2 -0.146 

Mode II: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 -0.329 -0.409 -0.517 -0.664 -0.866 -1.121 -1.410 

 tip 2 -0.329 -0.379 -0.433 -0.496 -0.572 -0.667 -0.789 

0.3 tip 1 -0.458 -0.567 -0.716 -0.924 -1.188 -1.491 

 tip 2 -0.458 -0.512 -0.575 -0.651 -0.747 -0.871 

0.4 tip 1 -0.622 -0.774 -0.987 -1.260 -1.580 

 tip 2 -0.622 -0.686 -0.765 -0.863 -0.992 

0.5 tip 1 -0.842 -1.061 -1.343 -1.680 

 tip 2 -0.843 -0.924 -1.028 -1.164 

0.6 tip 1 -1.150 -1.444 -1.800 

 tip 2 -1.151 -1.263 -1.409 

0.7 tip 1 -1.570 -1.950 

 tip 2 -1.572 -1.737 

0.8 tip 1 -2.147 

 tip 2 -2.147 

Mode III: all zero 
 45 

<111/112>, θ = 8°

Mode I: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.255 0.292 0.325 0.349 0.347 0.280 0.048 

 tip 2 0.255 0.295 0.336 0.377 0.418 0.455 0.477 

0.3 tip 1 0.333 0.368 0.394 0.396 0.332 0.104 

 tip 2 0.333 0.375 0.419 0.462 0.501 0.527 

0.4 tip 1 0.412 0.441 0.445 0.385 0.161 

 tip 2 0.412 0.458 0.503 0.545 0.573 

0.5 tip 1 0.488 0.496 0.438 0.218 

 tip 2 0.489 0.536 0.580 0.610 

0.6 tip 1 0.546 0.491 0.273 

 tip 2 0.546 0.591 0.623 

0.7 tip 1 0.537 0.321 

 tip 2 0.539 0.571 

0.8 tip 1 0.350 

 tip 2 0.351 

Mode II: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 -0.187 -0.240 -0.319 -0.439 -0.631 -0.936 -1.391 

 tip 2 -0.187 -0.217 -0.250 -0.291 -0.344 -0.419 -0.530 

0.3 tip 1 -0.270 -0.349 -0.470 -0.666 -0.978 -1.447 

 tip 2 -0.270 -0.303 -0.343 -0.396 -0.470 -0.582 

0.4 tip 1 -0.383 -0.506 -0.706 -1.024 -1.506 

 tip 2 -0.383 -0.424 -0.477 -0.552 -0.666 

0.5 tip 1 -0.549 -0.753 -1.079 -1.573 

 tip 2 -0.549 -0.604 -0.682 -0.799 

0.6 tip 1 -0.813 -1.147 -1.657 

 tip 2 -0.812 -0.896 -1.021 

0.7 tip 1 -1.240 -1.773 

 tip 2 -1.240 -1.378 

0.8 tip 1 -1.935 

 tip 2 -1.939 

Mode III: all zero 
46 

<111/112>, θ = 4°

Mode I: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.296 0.345 0.397 0.454 0.514 0.565 0.540 

 tip 2 0.296 0.342 0.391 0.442 0.499 0.560 0.626 

0.3 tip 1 0.392 0.447 0.507 0.572 0.630 0.616 

 tip 2 0.392 0.443 0.497 0.555 0.620 0.689 

0.4 tip 1 0.500 0.564 0.633 0.698 0.693 

 tip 2 0.500 0.557 0.618 0.686 0.758 

0.5 tip 1 0.624 0.698 0.769 0.773 

 tip 2 0.624 0.689 0.759 0.834 

0.6 tip 1 0.767 0.843 0.856 

 tip 2 0.768 0.842 0.919 

0.7 tip 1 0.922 0.941 

 tip 2 0.923 1.004 

0.8 tip 1 1.026 

 tip 2 1.030 

Mode II: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 -0.041 -0.061 -0.096 -0.154 -0.260 -0.471 -0.948 

 tip 2 -0.041 -0.050 -0.060 -0.075 -0.096 -0.131 -0.198 

0.3 tip 1 -0.070 -0.104 -0.163 -0.271 -0.485 -0.971 

 tip 2 -0.070 -0.080 -0.094 -0.115 -0.149 -0.213 

0.4 tip 1 -0.115 -0.174 -0.283 -0.500 -0.993 

 tip 2 -0.115 -0.128 -0.148 -0.182 -0.246 

0.5 tip 1 -0.188 -0.299 -0.518 -1.019 

 tip 2 -0.188 -0.208 -0.242 -0.307 

0.6 tip 1 -0.320 -0.543 -1.049 

 tip 2 -0.319 -0.354 -0.422 

0.7 tip 1 -0.582 -1.097 

 tip 2 -0.578 -0.652 

0.8 tip 1 -1.186 

 tip 2 -1.178 

Mode III: all zero 
 47 

<111/112>, θ = 0°

Mode I: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.316 0.370 0.429 0.499 0.585 0.698 0.871 

 tip 2 0.316 0.366 0.418 0.474 0.537 0.609 0.699 

0.3 tip 1 0.421 0.483 0.556 0.648 0.769 0.956 

 tip 2 0.421 0.475 0.534 0.599 0.675 0.769 

0.4 tip 1 0.540 0.617 0.715 0.844 1.045 

 tip 2 0.540 0.602 0.671 0.750 0.848 

0.5 tip 1 0.683 0.786 0.924 1.138 

 tip 2 0.684 0.757 0.840 0.942 

0.6 tip 1 0.865 1.011 1.236 

 tip 2 0.865 0.954 1.059 

0.7 tip 1 1.106 1.342 

 tip 2 1.106 1.217 

0.8 tip 1 1.461 

 tip 2 1.461 

Mode II: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.095 0.107 0.117 0.128 0.141 0.157 0.179 

 tip 2 0.095 0.107 0.118 0.129 0.140 0.153 0.168 

0.3 tip 1 0.118 0.129 0.140 0.154 0.171 0.197 

 tip 2 0.118 0.129 0.141 0.154 0.168 0.184 

0.4 tip 1 0.141 0.153 0.168 0.186 0.214 

 tip 2 0.141 0.154 0.168 0.184 0.202 

0.5 tip 1 0.167 0.184 0.204 0.235 

 tip 2 0.168 0.183 0.201 0.222 

0.6 tip 1 0.201 0.224 0.255 

 tip 2 0.202 0.222 0.246 

0.7 tip 1 0.247 0.287 

 tip 2 0.248 0.276 

0.8 tip 1 0.319 

 tip 2 0.324 

Mode III: all zero 
48 

<111/112>, θ = −4°

Mode I: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.319 0.371 0.427 0.489 0.559 0.634 0.693 

 tip 2 0.319 0.369 0.421 0.476 0.537 0.604 0.681 

0.3 tip 1 0.422 0.481 0.547 0.622 0.705 0.777 

 tip 2 0.422 0.476 0.535 0.598 0.670 0.751 

0.4 tip 1 0.538 0.608 0.689 0.781 0.865 

 tip 2 0.538 0.600 0.668 0.743 0.830 

0.5 tip 1 0.674 0.761 0.861 0.960 

 tip 2 0.674 0.747 0.827 0.920 

0.6 tip 1 0.840 0.948 1.057 

 tip 2 0.840 0.926 1.025 

0.7 tip 1 1.042 1.164 

 tip 2 1.042 1.149 

0.8 tip 1 1.281 

 tip 2 1.280 

Mode II: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.213 0.250 0.296 0.359 0.455 0.615 0.937 

 tip 2 0.213 0.241 0.271 0.303 0.342 0.391 0.463 

0.3 tip 1 0.278 0.325 0.389 0.487 0.652 0.983 

 tip 2 0.278 0.308 0.341 0.381 0.432 0.506 

0.4 tip 1 0.356 0.421 0.524 0.694 1.032 

 tip 2 0.356 0.391 0.433 0.486 0.563 

0.5 tip 1 0.459 0.564 0.740 1.092 

 tip 2 0.459 0.504 0.561 0.643 

0.6 tip 1 0.614 0.797 1.156 

 tip 2 0.614 0.676 0.766 

0.7 tip 1 0.867 1.241 

 tip 2 0.868 0.968 

0.8 tip 1 1.361 

 tip 2 1.361 

Mode III: all zero 
 49 

<111/112>, θ = −8°

Mode I: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.308 0.354 0.400 0.444 0.479 0.485 0.386 

 tip 2 0.308 0.356 0.405 0.456 0.509 0.564 0.615 

0.3 tip 1 0.404 0.452 0.499 0.539 0.552 0.464 

 tip 2 0.404 0.455 0.509 0.566 0.625 0.682 

0.4 tip 1 0.506 0.557 0.603 0.623 0.544 

 tip 2 0.506 0.564 0.625 0.689 0.752 

0.5 tip 1 0.619 0.671 0.698 0.630 

 tip 2 0.619 0.685 0.754 0.825 

0.6 tip 1 0.743 0.777 0.720 

 tip 2 0.742 0.818 0.895 

0.7 tip 1 0.861 0.815 

 tip 2 0.861 0.946 

0.8 tip 1 0.910 

 tip 2 0.909 

Mode II: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.309 0.365 0.434 0.526 0.661 0.868 1.220 

 tip 2 0.309 0.352 0.396 0.444 0.501 0.573 0.671 

0.3 tip 1 0.407 0.477 0.571 0.710 0.923 1.286 

 tip 2 0.407 0.451 0.501 0.559 0.633 0.733 

0.4 tip 1 0.523 0.619 0.763 0.984 1.359 

 tip 2 0.523 0.575 0.636 0.712 0.818 

0.5 tip 1 0.674 0.823 1.052 1.441 

 tip 2 0.675 0.740 0.822 0.934 

0.6 tip 1 0.895 1.133 1.536 

 tip 2 0.895 0.984 1.106 

0.7 tip 1 1.234 1.656 

 tip 2 1.234 1.371 

0.8 tip 1 1.820 

 tip 2 1.820 

Mode III: all zero 
50 

<111/112>, θ = −12°

Mode I: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.287 0.325 0.358 0.378 0.375 0.311 0.086 

 tip 2 0.287 0.331 0.376 0.420 0.464 0.504 0.527 

0.3 tip 1 0.371 0.405 0.429 0.430 0.371 0.151 

 tip 2 0.371 0.417 0.465 0.512 0.556 0.587 

0.4 tip 1 0.455 0.482 0.488 0.435 0.223 

 tip 2 0.455 0.506 0.557 0.607 0.644 

0.5 tip 1 0.537 0.548 0.501 0.297 

 tip 2 0.537 0.594 0.648 0.693 

0.6 tip 1 0.610 0.569 0.373 

 tip 2 0.610 0.671 0.723 

0.7 tip 1 0.637 0.449 

 tip 2 0.638 0.698 

0.8 tip 1 0.520 

 tip 2 0.520 

Mode II: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.386 0.455 0.538 0.644 0.793 1.005 1.319 

 tip 2 0.386 0.439 0.495 0.555 0.624 0.709 0.819 

0.3 tip 1 0.507 0.592 0.701 0.854 1.073 1.400 

 tip 2 0.507 0.563 0.624 0.696 0.783 0.896 

0.4 tip 1 0.649 0.761 0.920 1.148 1.488 

 tip 2 0.649 0.713 0.787 0.878 0.997 

0.5 tip 1 0.829 0.994 1.231 1.587 

 tip 2 0.830 0.909 1.005 1.131 

0.6 tip 1 1.081 1.328 1.704 

 tip 2 1.081 1.186 1.324 

0.7 tip 1 1.447 1.845 

 tip 2 1.447 1.602 

0.8 tip 1 2.030 

 tip 2 2.029 

Mode III: all zero 
 51 

<111/112>, θ = −16°

Mode I: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.259 0.288 0.304 0.301 0.258 0.128 -0.194 

 tip 2 0.259 0.299 0.337 0.374 0.406 0.429 0.426 

0.3 tip 1 0.328 0.347 0.346 0.306 0.180 -0.139 

 tip 2 0.328 0.369 0.408 0.444 0.472 0.476 

0.4 tip 1 0.390 0.392 0.355 0.234 -0.080 

 tip 2 0.390 0.433 0.473 0.506 0.518 

0.5 tip 1 0.438 0.406 0.289 -0.020 

 tip 2 0.438 0.483 0.521 0.541 

0.6 tip 1 0.456 0.345 0.042 

 tip 2 0.456 0.500 0.527 

0.7 tip 1 0.394 0.098 

 tip 2 0.394 0.428 

0.8 tip 1 0.139 

 tip 2 0.138 

Mode II: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.446 0.523 0.615 0.728 0.878 1.078 1.340 

 tip 2 0.446 0.508 0.572 0.641 0.719 0.811 0.926 

0.3 tip 1 0.584 0.677 0.793 0.949 1.157 1.432 

 tip 2 0.585 0.649 0.719 0.799 0.894 1.012 

0.4 tip 1 0.743 0.862 1.024 1.241 1.532 

 tip 2 0.744 0.816 0.900 0.999 1.124 

0.5 tip 1 0.940 1.108 1.336 1.643 

 tip 2 0.941 1.029 1.135 1.268 

0.6 tip 1 1.205 1.445 1.772 

 tip 2 1.205 1.321 1.467 

0.7 tip 1 1.574 1.926 

 tip 2 1.574 1.739 

0.8 tip 1 2.121 

 tip 2 2.120 

Mode III: all zero 
52 

<111/112>, θ = −20°

Mode I: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.225 0.243 0.243 0.213 0.129 -0.063 -0.447 

 tip 2 0.225 0.259 0.290 0.318 0.338 0.343 0.314 

0.3 tip 1 0.277 0.279 0.251 0.168 -0.022 -0.405 

 tip 2 0.277 0.311 0.341 0.365 0.374 0.353 

0.4 tip 1 0.315 0.289 0.209 0.022 -0.360 

 tip 2 0.315 0.348 0.376 0.390 0.376 

0.5 tip 1 0.325 0.249 0.065 -0.315 

 tip 2 0.326 0.357 0.377 0.371 

0.6 tip 1 0.285 0.105 -0.272 

 tip 2 0.285 0.310 0.311 

0.7 tip 1 0.134 -0.238 

 tip 2 0.134 0.142 

0.8 tip 1 -0.225 

 tip 2 -0.225 

Mode II: 

 2a2/W 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.493 0.576 0.672 0.787 0.934 1.115 1.316 

 tip 2 0.493 0.562 0.632 0.708 0.791 0.888 1.004 

0.3 tip 1 0.644 0.741 0.859 1.012 1.201 1.417 

 tip 2 0.644 0.715 0.792 0.878 0.978 1.097 

0.4 tip 1 0.814 0.936 1.095 1.293 1.525 

 tip 2 0.815 0.894 0.984 1.088 1.215 

0.5 tip 1 1.020 1.186 1.396 1.644 

 tip 2 1.021 1.116 1.228 1.363 

0.6 tip 1 1.290 1.512 1.781 

 tip 2 1.290 1.412 1.560 

0.7 tip 1 1.649 1.943 

 tip 2 1.649 1.817 

0.8 tip 1 2.143 

 tip 2 2.142 

Mode III: all zero 
 53 

Isotropic, θ = 20°

Mode I: 

 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.175 0.173 0.144 0.071 -0.070 -0.304 -0.626 

 tip 2 0.175 0.200 0.220 0.234 0.235 0.214 0.161 

0.3 tip 1 0.199 0.171 0.098 -0.043 -0.279 -0.607 

 tip 2 0.199 0.221 0.236 0.240 0.224 0.176 

0.4 tip 1 0.194 0.122 -0.019 -0.256 -0.588 

 tip 2 0.194 0.211 0.217 0.204 0.161 

0.5 tip 1 0.141 0.000 -0.239 -0.574 

 tip 2 0.141 0.148 0.138 0.098 

0.6 tip 1 0.008 -0.231 -0.570 

 tip 2 0.008 -0.001 -0.038 

0.7 tip 1 -0.242 -0.584 

 tip 2 -0.242 -0.279 

0.8 tip 1 -0.623 

 tip 2 -0.623 

Mode II: 

 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 -0.475 -0.570 -0.689 -0.835 -1.011 -1.192 -1.338 

 tip 2 -0.475 -0.543 -0.616 -0.696 -0.788 -0.896 -1.020 

0.3 tip 1 -0.638 -0.758 -0.908 -1.091 -1.281 -1.443 

 tip 2 -0.638 -0.711 -0.793 -0.887 -0.997 -1.125 

0.4 tip 1 -0.834 -0.987 -1.177 -1.378 -1.557 

 tip 2 -0.834 -0.918 -1.016 -1.130 -1.265 

0.5 tip 1 -1.077 -1.274 -1.486 -1.684 

 tip 2 -1.077 -1.180 -1.301 -1.446 

0.6 tip 1 -1.386 -1.611 -1.830 

 tip 2 -1.386 -1.518 -1.676 

0.7 tip 1 -1.759 -2.005 

 tip 2 -1.760 -1.938 

0.8 tip 1 -2.217 

 tip 2 -2.217 

Mode III: all zero 
54 

Isotropic, θ = 16°

Mode I: 

 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.239 0.256 0.253 0.215 0.115 -0.092 -0.451 

 tip 2 0.239 0.275 0.307 0.336 0.355 0.358 0.328 

0.3 tip 1 0.292 0.290 0.254 0.156 -0.051 -0.410 

 tip 2 0.292 0.327 0.357 0.380 0.387 0.363 

0.4 tip 1 0.327 0.293 0.196 -0.009 -0.370 

 tip 2 0.327 0.360 0.386 0.397 0.378 

0.5 tip 1 0.328 0.234 0.030 -0.331 

 tip 2 0.328 0.357 0.371 0.357 

0.6 tip 1 0.265 0.061 -0.299 

 tip 2 0.265 0.282 0.272 

0.7 tip 1 0.080 -0.282 

 tip 2 0.079 0.071 

0.8 tip 1 -0.290 

 tip 2 -0.290 

Mode II: 

 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 -0.394 -0.479 -0.589 -0.735 -0.933 -1.176 -1.422 

 tip 2 -0.394 -0.451 -0.513 -0.583 -0.665 -0.768 -0.895 

0.3 tip 1 -0.535 -0.647 -0.796 -1.000 -1.252 -1.514 

 tip 2 -0.535 -0.597 -0.667 -0.752 -0.856 -0.986 

0.4 tip 1 -0.710 -0.863 -1.073 -1.336 -1.614 

 tip 2 -0.710 -0.783 -0.870 -0.978 -1.113 

0.5 tip 1 -0.940 -1.157 -1.430 -1.726 

 tip 2 -0.941 -1.032 -1.146 -1.289 

0.6 tip 1 -1.257 -1.542 -1.858 

 tip 2 -1.257 -1.380 -1.535 

0.7 tip 1 -1.679 -2.021 

 tip 2 -1.680 -1.855 

0.8 tip 1 -2.230 

 tip 2 -2.232 

Mode III: all zero 
 55 

Isotropic, θ = 12°

Mode I: 

 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.291 0.324 0.346 0.346 0.302 0.159 -0.184 

 tip 2 0.292 0.335 0.378 0.420 0.457 0.484 0.485 

0.3 tip 1 0.369 0.392 0.395 0.354 0.214 -0.127 

 tip 2 0.369 0.414 0.458 0.499 0.531 0.538 

0.4 tip 1 0.440 0.446 0.408 0.273 -0.065 

 tip 2 0.440 0.487 0.532 0.568 0.581 

0.5 tip 1 0.496 0.462 0.331 -0.003 

 tip 2 0.496 0.544 0.585 0.603 

0.6 tip 1 0.514 0.386 0.055 

 tip 2 0.515 0.559 0.582 

0.7 tip 1 0.434 0.104 

 tip 2 0.433 0.460 

0.8 tip 1 0.133 

 tip 2 0.134 

Mode II: 

 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 -0.304 -0.373 -0.465 -0.597 -0.793 -1.079 -1.455 

 tip 2 -0.304 -0.349 -0.397 -0.453 -0.521 -0.611 -0.735 

0.3 tip 1 -0.417 -0.510 -0.645 -0.846 -1.140 -1.531 

 tip 2 -0.417 -0.465 -0.521 -0.591 -0.682 -0.807 

0.4 tip 1 -0.560 -0.697 -0.903 -1.205 -1.612 

 tip 2 -0.560 -0.618 -0.690 -0.783 -0.912 

0.5 tip 1 -0.759 -0.970 -1.282 -1.706 

 tip 2 -0.759 -0.834 -0.932 -1.067 

0.6 tip 1 -1.052 -1.374 -1.817 

 tip 2 -1.052 -1.158 -1.303 

0.7 tip 1 -1.493 -1.959 

 tip 2 -1.494 -1.656 

0.8 tip 1 -2.153 

 tip 2 -2.152 

Mode III: all zero 
56 

Isotropic, θ = 8°

Mode I: 

 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.330 0.375 0.417 0.451 0.467 0.426 0.204 

 tip 2 0.330 0.380 0.431 0.482 0.535 0.586 0.625 

0.3 tip 1 0.427 0.471 0.508 0.529 0.494 0.281 

 tip 2 0.427 0.479 0.534 0.590 0.646 0.691 

0.4 tip 1 0.526 0.568 0.593 0.565 0.359 

 tip 2 0.526 0.585 0.645 0.705 0.755 

0.5 tip 1 0.630 0.661 0.638 0.441 

 tip 2 0.630 0.695 0.759 0.814 

0.6 tip 1 0.730 0.713 0.525 

 tip 2 0.730 0.799 0.860 

0.7 tip 1 0.787 0.605 

 tip 2 0.788 0.854 

0.8 tip 1 0.675 

 tip 2 0.675 

Mode II: 

 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 -0.207 -0.255 -0.322 -0.422 -0.583 -0.855 -1.333 

 tip 2 -0.207 -0.237 -0.271 -0.310 -0.359 -0.428 -0.533 

0.3 tip 1 -0.285 -0.353 -0.455 -0.619 -0.898 -1.388 

 tip 2 -0.285 -0.318 -0.358 -0.408 -0.477 -0.582 

0.4 tip 1 -0.388 -0.491 -0.660 -0.944 -1.447 

 tip 2 -0.388 -0.428 -0.479 -0.550 -0.658 

0.5 tip 1 -0.534 -0.706 -0.998 -1.513 

 tip 2 -0.534 -0.588 -0.662 -0.774 

0.6 tip 1 -0.765 -1.066 -1.593 

 tip 2 -0.765 -0.844 -0.964 

0.7 tip 1 -1.154 -1.704 

 tip 2 -1.152 -1.285 

0.8 tip 1 -1.859 

 tip 2 -1.863 

Mode III: all zero 
 57 

Isotropic, θ = 4°

Mode I: 

 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.353 0.407 0.462 0.520 0.583 0.645 0.663 

 tip 2 0.354 0.408 0.463 0.521 0.584 0.653 0.730 

0.3 tip 1 0.462 0.520 0.582 0.651 0.721 0.752 

 tip 2 0.462 0.520 0.581 0.648 0.721 0.803 

0.4 tip 1 0.581 0.648 0.722 0.800 0.844 

 tip 2 0.581 0.646 0.717 0.794 0.881 

0.5 tip 1 0.717 0.798 0.884 0.939 

 tip 2 0.717 0.793 0.875 0.966 

0.6 tip 1 0.879 0.972 1.039 

 tip 2 0.879 0.967 1.063 

0.7 tip 1 1.066 1.140 

 tip 2 1.065 1.167 

0.8 tip 1 1.251 

 tip 2 1.251 

Mode II: 

 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 -0.105 -0.130 -0.165 -0.220 -0.312 -0.484 -0.875 

 tip 2 -0.105 -0.120 -0.137 -0.157 -0.184 -0.221 -0.285 

0.3 tip 1 -0.145 -0.181 -0.236 -0.330 -0.506 -0.904 

 tip 2 -0.145 -0.162 -0.182 -0.209 -0.246 -0.310 

0.4 tip 1 -0.199 -0.255 -0.350 -0.530 -0.935 

 tip 2 -0.198 -0.219 -0.246 -0.285 -0.349 

0.5 tip 1 -0.277 -0.375 -0.558 -0.970 

 tip 2 -0.277 -0.305 -0.345 -0.411 

0.6 tip 1 -0.405 -0.593 -1.014 

 tip 2 -0.405 -0.448 -0.518 

0.7 tip 1 -0.641 -1.071 

 tip 2 -0.639 -0.718 

0.8 tip 1 -1.164 

 tip 2 -1.161 

Mode III: all zero 
58 

Isotropic, θ = 0°

Mode I: 

 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.361 0.418 0.477 0.545 0.626 0.732 0.898 

 tip 2 0.361 0.417 0.474 0.534 0.601 0.677 0.770 

0.3 tip 1 0.474 0.537 0.608 0.695 0.810 0.992 

 tip 2 0.474 0.534 0.597 0.667 0.748 0.845 

0.4 tip 1 0.599 0.675 0.769 0.893 1.088 

 tip 2 0.599 0.667 0.741 0.826 0.928 

0.5 tip 1 0.748 0.847 0.980 1.189 

 tip 2 0.748 0.827 0.916 1.023 

0.6 tip 1 0.933 1.074 1.295 

 tip 2 0.933 1.027 1.138 

0.7 tip 1 1.175 1.406 

 tip 2 1.174 1.291 

0.8 tip 1 1.529 

 tip 2 1.530 

Mode II: 

 

2a1/W 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0.2 tip 1 0.000 0.000 0.000 -0.001 -0.001 -0.001 -0.007 

 tip 2 0.000 0.000 0.000 0.000 0.000 -0.001 -0.001 

0.3 tip 1 0.000 0.000 -0.001 -0.001 -0.001 -0.006 

 tip 2 0.000 0.000 0.000 0.000 0.000 -0.001 

0.4 tip 1 0.000 -0.001 -0.001 -0.001 -0.007 

 tip 2 0.000 0.000 0.000 0.000 -0.001 

0.5 tip 1 -0.001 -0.001 -0.001 -0.003 

 tip 2 0.000 0.000 0.000 -0.001 

0.6 tip 1 -0.001 -0.002 -0.008 

 tip 2 -0.001 0.000 -0.001 

0.7 tip 1 -0.002 -0.004 

 tip 2 -0.001 -0.001 

0.8 tip 1 -0.005 

 tip 2 -0.002 

Mode III: all zero 
 59 

<111/110>, θ = 20°

-2.500

-2.000

-1.500

-1.000

-0.500

0.000

0.500

1.000

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

-2.000

-1.500

-1.000

-0.500

0.000

0.500

1.000

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

Mode I: Crack Tip a2

Mode I: Crack Tip a1

2a2/W

2a2/W

-6.000

-5.000

-4.000

-3.000

-2.000

-1.000

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0 .2 0.3 0.4 0 .5 0.6

0 .7 0.8

Mode II: Crack Tip a1

2a2/W

-5.000

-4.500

-4.000

-3.500

-3.000

-2.500

-2.000

-1.500

-1.000

-0.500

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0.2 0.3 0 .4 0.5 0 .6

0.7 0.8

Mode II: Crack Tip a2

2a2/W
60 

<111/110>, θ = 16°

-1.500

-1.000

-0.500

0.000

0.500

1.000

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

-1.500

-1.000

-0.500

0.000

0.500

1.000

1.500

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

Mode I: Crack Tip a2

Mode I: Crack Tip a1

2a2/W

2a2/W

-6.000

-5.000

-4.000

-3.000

-2.000

-1.000

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0 .2 0.3 0.4 0 .5 0.6

0 .7 0.8

Mode II: Crack Tip a1

2a2/W

-5.000

-4.500

-4.000

-3.500

-3.000

-2.500

-2.000

-1.500

-1.000

-0.500

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0.2 0.3 0 .4 0.5 0 .6

0.7 0.8

Mode II: Crack Tip a2

2a2/W
 61 

<111/110>, θ = 12°

-0.600

-0.400

-0.200

0.000

0.200

0.400

0.600

0.800

1.000

1.200

1.400

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

-0.500

0.000

0.500

1.000

1.500

2.000

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

Mode I: Crack Tip a2

Mode I: Crack Tip a1

2a2/W

2a2/W

-6.000

-5.000

-4.000

-3.000

-2.000

-1.000

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0 .2 0.3 0.4 0 .5 0.6

0 .7 0.8

Mode II: Crack Tip a1

2a2/W

-5.000

-4.500

-4.000

-3.500

-3.000

-2.500

-2.000

-1.500

-1.000

-0.500

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0.2 0.3 0 .4 0.5 0 .6

0.7 0.8

Mode II: Crack Tip a2

2a2/W
62 

<111/110>, θ = 8°

0.000

0.200

0.400

0.600

0.800

1.000

1.200

1.400

1.600

1.800

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

0.000

0.500

1.000

1.500

2.000

2.500

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

Mode I: Crack Tip a2

Mode I: Crack Tip a1

2a2/W

2a2/W

-6.000

-5.000

-4.000

-3.000

-2.000

-1.000

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0 .2 0.3 0.4 0 .5 0.6

0 .7 0.8

Mode II: Crack Tip a1

2a2/W

-4.500

-4.000

-3.500

-3.000

-2.500

-2.000

-1.500

-1.000

-0.500

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0.2 0.3 0 .4 0.5 0 .6

0.7 0.8

Mode II: Crack Tip a2

2a2/W
 63 

<111/110>, θ = 4°

0.000

0.500

1.000

1.500

2.000

2.500

3.000

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

0.000

0.500

1.000

1.500

2.000

2.500

3.000

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

Mode I: Crack Tip a2

Mode I: Crack Tip a1

2a2/W

2a2/W

-4.500

-4.000

-3.500

-3.000

-2.500

-2.000

-1.500

-1.000

-0.500

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0 .2 0.3 0.4 0 .5 0.6

0 .7 0.8

Mode II: Crack Tip a1

2a2/W

-3.000

-2.500

-2.000

-1.500

-1.000

-0.500

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0.2 0.3 0 .4 0.5 0 .6

0.7 0.8

Mode II: Crack Tip a2

2a2/W
64 

<111/110>, θ =0°

0.000

0.500

1.000

1.500

2.000

2.500

3.000

3.500

4.000

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

0.000

0.500

1.000

1.500

2.000

2.500

3.000

3.500

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

Mode I: Crack Tip a2

Mode I: Crack Tip a1

2a2/W

2a2/W

-0.060

-0.050

-0.040

-0.030

-0.020

-0.010

0.000

0.010

0.2 0.4 0.6 0.8

2a1/W

FII

0 .2 0.3 0.4 0 .5 0.6

0 .7 0.8

Mode II: Crack Tip a1

2a2/W

-0.006

-0.004

-0.002

0.000

0.002

0.004

0.006

0.2 0.4 0.6 0.8

2a1/W

FII

0.2 0.3 0 .4 0.5 0 .6

0.7 0.8

Mode II: Crack Tip a2

2a2/W
 65 

<111/112>, θ = 20°

-3.000

-2.500

-2.000

-1.500

-1.000

-0.500

0.000

0.500

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

-2.500

-2.000

-1.500

-1.000

-0.500

0.000

0.500

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

Mode I: Crack Tip a2

Mode I: Crack Tip a1

2a2/W

2a2/W

-6.000

-5.000

-4.000

-3.000

-2.000

-1.000

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0 .2 0.3 0.4 0 .5 0.6

0 .7 0.8

Mode II: Crack Tip a1

2a2/W

-5.000

-4.500

-4.000

-3.500

-3.000

-2.500

-2.000

-1.500

-1.000

-0.500

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0.2 0.3 0 .4 0.5 0 .6

0.7 0.8

Mode II: Crack Tip a2

2a2/W
66 

<111/112>, θ = 16°

-2.500

-2.000

-1.500

-1.000

-0.500

0.000

0.500

1.000

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

-2.000

-1.500

-1.000

-0.500

0.000

0.500

1.000

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

Mode I: Crack Tip a2

Mode I: Crack Tip a1

2a2/W

-6.000

-5.000

-4.000

-3.000

-2.000

-1.000

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0 .2 0.3 0.4 0 .5 0.6

0 .7 0.8

Mode II: Crack Tip a1

2a2/W

-5.000

-4.500

-4.000

-3.500

-3.000

-2.500

-2.000

-1.500

-1.000

-0.500

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0.2 0.3 0 .4 0.5 0 .6

0.7 0.8

Mode II: Crack Tip a2

2a2/W
 67 

<111/112>, θ =12°

-1.500

-1.000

-0.500

0.000

0.500

1.000

1.500

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

-1.000

-0.500

0.000

0.500

1.000

1.500

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

Mode I: Crack Tip a2

Mode I: Crack Tip a1

2a2/W

2a2/W

-7.000

-6.000

-5.000

-4.000

-3.000

-2.000

-1.000

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0 .2 0.3 0.4 0 .5 0.6

0 .7 0.8

Mode II: Crack Tip a1

2a2/W

-5.000

-4.500

-4.000

-3.500

-3.000

-2.500

-2.000

-1.500

-1.000

-0.500

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0.2 0.3 0 .4 0.5 0 .6

0.7 0.8

Mode II: Crack Tip a2

2a2/W
68 

<111/112>, θ =8°

0.000

0.200

0.400

0.600

0.800

1.000

1.200

1.400

1.600

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

0.000

0.500

1.000

1.500

2.000

2.500

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

Mode I: Crack Tip a2

Mode I: Crack Tip a1

2a2/W

2a2/W

-7.000

-6.000

-5.000

-4.000

-3.000

-2.000

-1.000

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0 .2 0.3 0.4 0 .5 0.6

0 .7 0.8

Mode II: Crack Tip a1

2a2/W

-4.500

-4.000

-3.500

-3.000

-2.500

-2.000

-1.500

-1.000

-0.500

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0.2 0.3 0 .4 0.5 0 .6

0.7 0.8

Mode II: Crack Tip a2

2a2/W
 69 

<111/112>, θ = 4°

0.000

0.500

1.000

1.500

2.000

2.500

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

0.000

0.500

1.000

1.500

2.000

2.500

3.000

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

Mode I: Crack Tip a2

Mode I: Crack Tip a1

2a2/W

2a2/W

-4.500

-4.000

-3.500

-3.000

-2.500

-2.000

-1.500

-1.000

-0.500

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0 .2 0.3 0.4 0 .5 0.6

0 .7 0.8

Mode II: Crack Tip a1

2a2/W

-3.000

-2.500

-2.000

-1.500

-1.000

-0.500

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0.2 0.3 0 .4 0.5 0 .6

0.7 0.8

Mode II: Crack Tip a2

2a2/W
70 

<111/112>, θ = 0°

0.000

0.500

1.000

1.500

2.000

2.500

3.000

3.500

4.000

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

0.000

0.500

1.000

1.500

2.000

2.500

3.000

3.500

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

Mode I: Crack Tip a2

Mode I: Crack Tip a1

2a2/W

2a2/W

0.000

0.100

0.200

0.300

0.400

0.500

0.600

0.700

0.800

0.2 0.4 0.6 0.8

2a1/W

FII

0 .2 0.3 0.4 0 .5 0.6

0 .7 0.8

Mode II: Crack Tip a1

2a2/W

0.000

0.100

0.200

0.300

0.400

0.500

0.600

0.700

0.800

0.2 0.4 0.6 0.8

2a1/W

FII

0.2 0.3 0 .4 0.5 0 .6

0.7 0.8

Mode II: Crack Tip a2

2a2/W
 71 

<111/112>, θ = −4°

0.000

0.500

1.000

1.500

2.000

2.500

3.000

3.500

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

0.000

0.500

1.000

1.500

2.000

2.500

3.000

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

Mode I: Crack Tip a2

Mode I: Crack Tip a1

2a2/W

2a2/W

0.000

0.500

1.000

1.500

2.000

2.500

3.000

3.500

4.000

4.500

0.2 0.4 0.6 0.8

2a1/W

FII

0 .2 0.3 0.4 0 .5 0.6

0 .7 0.8

Mode II: Crack Tip a1

2a2/W

0.000

0.500

1.000

1.500

2.000

2.500

3.000

0.2 0.4 0.6 0.8

2a1/W

FII

0.2 0.3 0 .4 0.5 0 .6

0.7 0.8

Mode II: Crack Tip a2

2a2/W
72 

<111/112>, θ = −8°

0.000

0.500

1.000

1.500

2.000

2.500

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

0.000

0.500

1.000

1.500

2.000

2.500

3.000

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

Mode I: Crack Tip a2

Mode I: Crack Tip a1

2a2/W

2a2/W

0.000

1.000

2.000

3.000

4.000

5.000

6.000

0.2 0.4 0.6 0.8

2a1/W

FII

0 .2 0.3 0.4 0 .5 0.6

0 .7 0.8

Mode II: Crack Tip a1

2a2/W

0.000

0.500

1.000

1.500

2.000

2.500

3.000

3.500

4.000

0.2 0.4 0.6 0.8

2a1/W

FII

0.2 0.3 0 .4 0.5 0 .6

0.7 0.8

Mode II: Crack Tip a2

2a2/W
 73 

<111/112>, θ = −12°

0.000

0.200

0.400

0.600

0.800

1.000

1.200

1.400

1.600

1.800

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

0.000

0.500

1.000

1.500

2.000

2.500

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

Mode I: Crack Tip a2

Mode I: Crack Tip a1

2a2/W

2a2/W

0.000

1.000

2.000

3.000

4.000

5.000

6.000

0.2 0.4 0.6 0.8

2a1/W

FII

0 .2 0.3 0.4 0 .5 0.6

0 .7 0.8

Mode II: Crack Tip a1

2a2/W

0.000

0.500

1.000

1.500

2.000

2.500

3.000

3.500

4.000

4.500

0.2 0.4 0.6 0.8

2a1/W

FII

0.2 0.3 0 .4 0.5 0 .6

0.7 0.8

Mode II: Crack Tip a2

2a2/W
74 

<111/112>, θ = −16°

-1.000

-0.500

0.000

0.500

1.000

1.500

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

-1.000

-0.500

0.000

0.500

1.000

1.500

2.000

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

Mode I: Crack Tip a2

Mode I: Crack Tip a1

2a2/W

2a2/W

0.000

1.000

2.000

3.000

4.000

5.000

6.000

0.2 0.4 0.6 0.8

2a1/W

FII

0 .2 0.3 0.4 0 .5 0.6

0 .7 0.8

Mode II: Crack Tip a1

2a2/W

0.000

0.500

1.000

1.500

2.000

2.500

3.000

3.500

4.000

4.500

5.000

0.2 0.4 0.6 0.8

2a1/W

FII

0.2 0.3 0 .4 0.5 0 .6

0.7 0.8

Mode II: Crack Tip a2

2a2/W
 75 

<111/112>, θ = −20°

-2.500

-2.000

-1.500

-1.000

-0.500

0.000

0.500

1.000

1.500

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

-2.000

-1.500

-1.000

-0.500

0.000

0.500

1.000

1.500

2.000

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

Mode I: Crack Tip a2

Mode I: Crack Tip a1

2a2/W

2a2/W

0.000

1.000

2.000

3.000

4.000

5.000

6.000

0.2 0.4 0.6 0.8

2a1/W

FII

0 .2 0.3 0.4 0 .5 0.6

0 .7 0.8

Mode II: Crack Tip a1

2a2/W

0.000

0.500

1.000

1.500

2.000

2.500

3.000

3.500

4.000

4.500

5.000

0.2 0.4 0.6 0.8

2a1/W

FII

0.2 0.3 0 .4 0.5 0 .6

0.7 0.8

Mode II: Crack Tip a2

2a2/W
76 

Isotropic, θ = 20°

-3.000

-2.500

-2.000

-1.500

-1.000

-0.500

0.000

0.500

1.000

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

-2.500

-2.000

-1.500

-1.000

-0.500

0.000

0.500

1.000

1.500

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

Mode I: Crack Tip a2

Mode I: Crack Tip a1

2a2/W

2a2/W

-6.000

-5.000

-4.000

-3.000

-2.000

-1.000

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0 .2 0.3 0.4 0 .5 0.6

0 .7 0.8

Mode II: Crack Tip a1

2a2/W

-5.000

-4.500

-4.000

-3.500

-3.000

-2.500

-2.000

-1.500

-1.000

-0.500

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0.2 0.3 0 .4 0.5 0 .6

0.7 0.8

Mode II: Crack Tip a2

2a2/W
 77 

Isotropic, θ = 16°

-2.500

-2.000

-1.500

-1.000

-0.500

0.000

0.500

1.000

1.500

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

-2.000

-1.500

-1.000

-0.500

0.000

0.500

1.000

1.500

2.000

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

Mode I: Crack Tip a2

Mode I: Crack Tip a1

2a2/W

2a2/W

-7.000

-6.000

-5.000

-4.000

-3.000

-2.000

-1.000

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0 .2 0.3 0.4 0 .5 0.6

0 .7 0.8

Mode II: Crack Tip a1

2a2/W

-5.000

-4.500

-4.000

-3.500

-3.000

-2.500

-2.000

-1.500

-1.000

-0.500

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0.2 0.3 0 .4 0.5 0 .6

0.7 0.8

Mode II: Crack Tip a2

2a2/W
78 

Isotropic, θ = 12°

-1.000

-0.500

0.000

0.500

1.000

1.500

2.000

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

-1.000

-0.500

0.000

0.500

1.000

1.500

2.000

2.500

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

Mode I: Crack Tip a2

Mode I: Crack Tip a1

2a2/W

2a2/W

-7.000

-6.000

-5.000

-4.000

-3.000

-2.000

-1.000

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0 .2 0.3 0.4 0 .5 0.6

0 .7 0.8

Mode II: Crack Tip a1

2a2/W

-5.000

-4.500

-4.000

-3.500

-3.000

-2.500

-2.000

-1.500

-1.000

-0.500

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0.2 0.3 0 .4 0.5 0 .6

0.7 0.8

Mode II: Crack Tip a2

2a2/W
 79 

Isotropic, θ = 8°

0.000

0.500

1.000

1.500

2.000

2.500

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

0.000

0.500

1.000

1.500

2.000

2.500

3.000

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

Mode I: Crack Tip a2

Mode I: Crack Tip a1

2a2/W

2a2/W

-6.000

-5.000

-4.000

-3.000

-2.000

-1.000

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0 .2 0.3 0.4 0 .5 0.6

0 .7 0.8

Mode II: Crack Tip a1

2a2/W

-4.500

-4.000

-3.500

-3.000

-2.500

-2.000

-1.500

-1.000

-0.500

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0.2 0.3 0 .4 0.5 0 .6

0.7 0.8

Mode II: Crack Tip a2

2a2/W
80 

Isotropic, θ = 4°

0.000

0.500

1.000

1.500

2.000

2.500

3.000

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

0.000

0.500

1.000

1.500

2.000

2.500

3.000

3.500

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

Mode I: Crack Tip a2

Mode I: Crack Tip a1

2a2/W

2a2/W

-4.000

-3.500

-3.000

-2.500

-2.000

-1.500

-1.000

-0.500

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0 .2 0.3 0.4 0 .5 0.6

0 .7 0.8

Mode II: Crack Tip a1

2a2/W

-3.000

-2.500

-2.000

-1.500

-1.000

-0.500

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0.2 0.3 0 .4 0.5 0 .6

0.7 0.8

Mode II: Crack Tip a2

2a2/W
 81 

Isotropic, θ = 0°

0.000

0.500

1.000

1.500

2.000

2.500

3.000

3.500

4.000

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

0.000

0.500

1.000

1.500

2.000

2.500

3.000

3.500

4.000

0.2 0.4 0.6 0.8

2a1/W

FI

0 .2 0.3 0 .4 0.5 0.6

0 .7 0.8

Mode I: Crack Tip a2

Mode I: Crack Tip a1

2a2/W

2a2/W

-0.035

-0.030

-0.025

-0.020

-0.015

-0.010

-0.005

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0 .2 0.3 0.4 0 .5 0.6

0 .7 0.8

Mode II: Crack Tip a1

2a2/W

-0.005

-0.004

-0.004

-0.003

-0.003

-0.002

-0.002

-0.001

-0.001

0.000

0.2 0.4 0.6 0.8

2a1/W

FII

0.2 0.3 0 .4 0.5 0 .6

0.7 0.8

Mode II: Crack Tip a2

2a2/W
vi

REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions 

for reducing this burden, to Washington Headquarters Services, Directorate for Information Operation and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office 

of Management and Budget, Paperwork Reduction Project (0704-0188), Washington, DC 20503

1. AGENCY USE ONLY (Leave Blank) 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED

4. TITLE AND SUBTITLE 5. FUNDING NUMBERS

6. AUTHORS

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

 REPORT NUMBER

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

 AGENCY REPO NUMBER

11. SUPPLEMENTARY NOTES

12a. DISTRIBUTION/AVAILABILITY STATEMENT 12b. DISTRIBUTION CODE

13. ABSTRACT (Maximum 200 words)

14. SUBJECT TERMS 15. NUMBER OF PAGES

16. PRICE CODE

17. SECURITY CLASSIFICATION

 OF REPORT

18. SECURITY CLASSIFICATION

 OF THIS PAGE

19. SECURITY CLASSIFICATION

 OF ABSTRACT

20. LIMITATION OF ABSTRACT

NSN 7540-01-280-5500 Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Std. 239-18

298-102

Unclassified Unclassified Unclassified Unlimited

The M-Integral for Computing Stress Intensity Factors 

in Generally Anisotropic Materials

P.A. Wawrzynek,* B.J. Carter,* and L. Banks-Sills**

*Fracture Analysis Consultants, Inc. **Tel Aviv University

121 Eastern Heights Dr. Tel Aviv, Israel

Ithaca, NY 14850 

NASA’s Space Shuttle Main Engine (SSME) Program

George C. Marshall Space Flight Center

Marshall Space Flight Center, AL 35812

Prepared for NASA’s Space Shuttle Main Engine (SSME) Program by Fracture Analysis Consultants, Inc.

Technical Monitor: G.R. Swanson

Unclassified-Unlimited

Subject Category 39

Availability: NASA CASI 301–621–0390

Single-crystal super alloys are commonly used for components in the hot sections of contemporary jet and rocket engines. Due to the anisotropic nature of single-crystal materials the use of 

existing isotropic fracture mechanics calculations leads to errors in stress intensity factors. The 

difference can be substantial.

Presented in this report is the solution for calculating stress intensity factors in generally anisotropic materials using the M-integral. Included are examples of this solution applied to Brazilian 

disk crack growth specimens. 

88

M–1141

NAS8–36801

July 2005 Contractor Report

NASA/CR—2005–214006

single-crystal, anisotropic, fracture mechanics, Brazilian disk, 

stress intensity fractures
The NASA STI Program Office...in Profile

Since its founding, NASA has been dedicated to

the advancement of aeronautics and space

science. The NASA Scientific and Technical 

Information (STI) Program Office plays a key

part in helping NASA maintain this important

role.

The NASA STI Program Office is operated by Langley 

Research Center, the lead center for NASA’s scientific 

and technical information. The NASA STI Program 

Office provides access to the NASA STI Database, the 

largest collection of aeronautical and space science 

STI in the world. The Program Office is also NASA’s 

institutional mechanism for disseminating the results 

of its research and development activities. These 

results are published by NASA in the NASA STI 

Report Series, which includes the following report 

types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major significant 

phase of research that present the results of 

NASA programs and include extensive data 

or theoretical analysis. Includes compilations 

of significant scientific and technical data 

and information deemed to be of continuing 

reference value. NASA’s counterpart of peerreviewed formal professional papers but has 

less stringent limitations on manuscript length 

and extent of graphic presentations.

• TECHNICAL MEMORANDUM. Scientific 

and technical findings 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. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientific and technical 

conferences, symposia, seminars, or other 

meetings sponsored or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientific, 

technical, or historical information from 

NASA programs, projects, and mission, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. 

English-language translations of foreign 

scientific and technical material pertinent to 

NASA’s mission.

Specialized services that complement the STI Program 

Office’s diverse offerings include creating custom 

thesauri, building customized databases, organizing 

and publishing research results...even providing 

videos.

For more information about the NASA STI Program 

Office, see the following:

• Access the NASA STI Program Home Page at 

http://www.sti.nasa.gov

• E-mail your question via the Internet to 

help@sti.nasa.gov

• Fax your question to the NASA Access Help 

Desk at 301–621–0134

• Telephone the NASA Access Help Desk at 

301–621–0390

• Write to:

NASA Access Help Desk

NASA Center for AeroSpace Information

7121 Standard Drive

Hanover, MD 21076–1320

301–621–0390
National Aeronautics and

Space Administration

IS04

George C. Marshall Space Flight Center

Marshall Space Flight Center, Alabama

35812