Thesis Melson

Fatigue Crack Growth Analysis
with Finite Element Methods

and a Monte Carlo Simulation
Joshua H. Melson
Thesis submitted to the Faculty of the

Virginia Polytechnic Institute and State University
in partial fulllment of the requirements for the degree of
Master of Science
In
Mechanical Engineering
Robert L. West, Chair

Norman E. Dowling, Co-Chair
John M. Kennedy
May 2, 2014
Blacksburg, VA
Keywords: Fracture, eXtended Finite Element Method,

Phantom Nodes, Monte Carlo
Copyright 2014, Joshua H. Melson
Fatigue Crack Growth Analysis

with Finite Element Methods
and a Monte Carlo Simulation
Joshua H. Melson
Abstract
Fatigue crack growth in engineered structures reduces the structures load carrying capacity
and will eventually lead to failure. Cycles required to grow a crack from an initial length to
the critical length is called the fatigue fracture life. In this thesis, ve dierent methods for
analyzing the fatigue fracture life of a center cracked plate were compared to experimental
data previously collected by C.M. Hudson in a 1969 NASA report studying the
R-ratio eects
on crack growth in 7075-T6 aluminum alloy. The Paris, Walker, and Forman fatigue crack
growth models were t the experimental data. The Walker equation best t the data since it
incorporated
R-ratio
eects and had a similar Root Mean Square Error (RMSE) compared
to the other models. There was insucient data in the unstable region of crack growth to
adequately t the Forman equation.
Analytical models were used as a baseline for all fatigue fracture life comparisons.
Life
estimates from AFGROW and nite elements with mid-side nodes moved to their quarter
point location compared very with the analytical model with errors less than 3%. The Virtual
Crack Closure Technique (VCCT) was selected as a method for crack propagation along a
predened path. Stress intensity factors (SIFs) for shorter crack lengths were found to be
low, resulting in an overestimated life of about 8%. The eXtended Finite Element Method
with Phantom Nodes (XFEM-PN) was used, allowing crack propagation along a solution
dependent path, independent of the mesh.
Low SIFs throughout growth resulted in life
estimates 20% too large. All nite element analyses were performed in Abaqus 6-13.3. An
integrated polynomial method was developed for calculating life based on Abaqus' results,
leading to coarser meshes with answers closer to the analytical estimate. None of the ve
methods for estimating life compared well with the experimental data, with analytical errors
on life ranging from 10-20%. These errors were attributed to the limited number of crack
growth experiments run at each
R-ratio,
and the large variability typically seen in growth
rates.
Monte Carlo simulations were run to estimate the distribution on life. It was shown that
material constants in the Walker model must be selected based on their interrelation with
a multivariate normal probability density function. Both analytical and XFEM-PN simulations had similar coecients of variation on life of approximately 3% with similar normal
distributions.
It was concluded that Abaqus' XFEM-PN is a reasonable means of estimating fatigue fracture
life and its variation, and this method could be extended to other geometries and threedimensional analyses.
Acknowledgments
I would like to thank my adviser and committee chair, Dr.
Bob West, for his technical
guidance, availability to answer questions, and patience throughout my graduate career.

Without his support, this would not have been possible.
I would also like to thank my committee co-chair, Dr.
Norman Dowling, who's technical
expertise in fracture mechanics was invaluable. He was always willing to meet and clear up
any questions encountered along the way.
Thank you to my committee member, Dr.
John Kennedy, who helped guide me through
undergraduate and graduate school.

Thanks to The MathWorks, the College of Engineering, and the Mechanical Engineering
Department of Virginia Tech for providing funding to pursue my master's degree.
iii
Contents
Introduction
1.1
Overview and Statement of Need
. . . . . . . . . . . . . . . . . . . . . . . .
1.2
Research Goals and Objectives . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3
Scope of Thesis
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4
Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Literature Review
2.1
Material Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.1
Fracture Mechanics . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.2
Mixed-Mode Critical Energy Release Rate
. . . . . . . . . . . . . . .
2.1.3
Fatigue Crack Growth Models . . . . . . . . . . . . . . . . . . . . . .
2.1.4
Onset of Crack Growth . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.5
Crack Extension Direction . . . . . . . . . . . . . . . . . . . . . . . .
2.2
2.3
Methods of Fracture Analysis
. . . . . . . . . . . . . . . . . . . . . . . . . .
10
2.2.1
Analytical Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
2.2.2
Finite Element Methods
11
. . . . . . . . . . . . . . . . . . . . . . . . .
Stochastic Analyses in Fracture Mechanics
. . . . . . . . . . . . . . . . . . .
18
Crack Growth Rates for 7075-T6 Aluminum Alloy
19
3.1
Experimental Data
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
19
3.2
Fatigue Crack Growth Models . . . . . . . . . . . . . . . . . . . . . . . . . .
21
3.2.1
Paris Crack Growth Model . . . . . . . . . . . . . . . . . . . . . . . .
21
3.2.2
Walker Crack Growth Model . . . . . . . . . . . . . . . . . . . . . . .
22
iv
3.2.3
3.3
. . . . . . . . . . . . . . . . . . . . . .
23
Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
25
Estimating Fatigue Fracture Life
26
4.1
Center Cracked Plate Model . . . . . . . . . . . . . . . . . . . . . . . . . . .
27
4.2
Analytical Method
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
27
4.2.1
Model Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
28
4.2.2
Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
28
AFGROW . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
31
4.3.1
Model Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
31
4.3.2
Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
32
4.3
4.4
4.5
4.6
4.7
Forman Crack Growth Model
Finite Elements with Quarter Points (FEQP)
. . . . . . . . . . . . . . . . .
35
4.4.1
Model Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
35
4.4.2
Convergence Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . .
37
4.4.3
Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
39
Virtual Crack Closure Technique (VCCT)
. . . . . . . . . . . . . . . . . . .
42
4.5.1
Model Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
42
4.5.2
47
4.5.3
Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
49
eXtended Finite Element Method with Phantom Nodes (XFEM-PN)
. . . .
52
4.6.1
Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
52
4.6.2
55
4.6.3
Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
56
Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
59
Monte Carlo Analysis
61
5.1
Random Numbers Selection
. . . . . . . . . . . . . . . . . . . . . . . . . . .
61
5.2
Analytical Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
65
5.3
XFEM-PN Simulation
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
68
5.4
Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
70
Conclusions and Future Work
71
6.1
Research Goals and Objectives . . . . . . . . . . . . . . . . . . . . . . . . . .
71
6.2
Conclusions
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
71
6.3
Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
73
Bibliography
79
Appendix A Material Properties
85
Appendix B Center Cracked Plate Analysis Results
89
Appendix C Convergence Analysis for Finite Element Models
103
Appendix D Monte Carlo Analysis Results
107
Appendix E Compact Specimen Analysis Results
110
Appendix F MATLAB Code
116
Appendix G Python Code
130
vi
List of Tables
3.1
Paris Equation Coecients . . . . . . . . . . . . . . . . . . . . . . . . . . . .
22
3.2
Walker Equation Coecients . . . . . . . . . . . . . . . . . . . . . . . . . . .
23
3.3
Forman Equation Coecients
. . . . . . . . . . . . . . . . . . . . . . . . . .
24
4.1
Analyzed Loading Conditions
. . . . . . . . . . . . . . . . . . . . . . . . . .
26
4.2
Results from the Analytical Fatigue Fracture Analysis . . . . . . . . . . . . .
29
4.3
AFGROW Material Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
32
4.4
AFGROW Load Spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . .
33
4.5
Results from the AFGROW Fatigue Fracture Analysis
. . . . . . . . . . . .
33
4.6
Convergence analysis for the FEQP model
. . . . . . . . . . . . . . . . . . .
38
4.7
RMSE of
. . . . .
38
4.8
Results from the FEQP Fatigue Fracture Analysis . . . . . . . . . . . . . . .
39
4.9
Abaqus Crack Growth Constants. . . . . . . . . . . . . . . . . . . . . . . . .
46
KI
Contour Values at
a = 1.0 in.
(25.4 mm) in Fig. 4.7.
4.10 Stress Intensity Range Convergence Analysis.
. . . . . . . . . . . . . . . . .
4.11 Estimated Fatigue Fracture Life Convergence Analysis.
48
. . . . . . . . . . . .
49
4.12 Results from the FEQP Fatigue Fracture Analysis . . . . . . . . . . . . . . .
50
4.13 XFEM-PN Stress Intensity Range Convergence Analysis. . . . . . . . . . . .
55
4.14 XFEM-PN Estimated Fatigue Fracture Life Convergence Analysis. . . . . . .
56
4.15 Results from the XFEM-PN Fatigue Fracture Analysis
. . . . . . . . . . . .
57
4.16 Comparison of Critical Crack Lengths . . . . . . . . . . . . . . . . . . . . . .
59
4.17 Comparison of Life Estimates
60
. . . . . . . . . . . . . . . . . . . . . . . . . .
6.1
Convergence Analysis for the Compact Specimen
A.1
Average Tensile Properties for 7075-T6 Aluminum Alloy
vii
. . . . . . . . . . . . . . .
. . . . . . . . . . .
76
85
A.2
Fatigue Crack Growth Test Results for 7075-T6 Aluminum Alloy . . . . . . .
86
A.3
Residual Static Strength Test Results for 7075-T6 Aluminum Alloy
. . . . .
87
B.1
. . . . . . . . . . . . . . . . . . . . . . . . . .
90
B.2
RMSE Compared to the Analytical Method for the Normalized

in Figs. B.1-B.4.
B.3
KI
Curves
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
90
RMSE Compared to the Analytical Method for the Crack Growth Curves in
Figs. B.5-B.8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
KI
C.1
RMSE of
Contour Values in Convergence Analysis.
E.1
Rail Steel Material Properties
E.2
Compact Specimen Loading and Boundary Conditions
90
. . . . . . . . . . .
104
. . . . . . . . . . . . . . . . . . . . . . . . . .
111
viii
. . . . . . . . . . . .
111
List of Figures
2.1
Three modes of displacement at a crack tip.
. . . . . . . . . . . . . . . . . .
2.2
Typical crack growth curve on log-log axes. . . . . . . . . . . . . . . . . . . .
2.3
Schematic mesh of elements around a crack with mid-side nodes moved to the
quarter point location.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4
Schematic nite element representation of the VCCT method.
2.5
XFEM nodes with enrichment and jump functions.
2.6
Level Set Method used for XFEM.
13
. . . . . . . .
14
. . . . . . . . . . . . . .
15
. . . . . . . . . . . . . . . . . . . . . . .
15
2.7
Schematic mesh with phantom nodes. . . . . . . . . . . . . . . . . . . . . . .
17
3.1
Center Cracked Plate Specimen
20
3.2
Paris equation t to Hudson's experimental data.
3.3
Walker equation t to Hudson's experimental data.
. . . . . . . . . . . . . .
23
3.4
Forman equation t to Hudson's experimental data. . . . . . . . . . . . . . .
25
4.1
Simplied Center Crack Plate Model
27
4.2
Comparison of cycle count during crack propagation for the analytical method
and experimental data (R
4.3
4.4
= 0.2).
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . .
KI during
(R = 0.2). . .
Comparison of the normalized
crack propagation for the analyt-
ical method and AFGROW
. . . . . . . . . . . . . . . . . . .
= 0.2).
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
33
34
36
Coarse meshing scheme for the FEQP model with an inset zoomed in at the
crack tip.
4.7
31
Partitioning, loading, and boundary conditions used in the FEQP model with
an inset zoomed in at the crack tip. . . . . . . . . . . . . . . . . . . . . . . .
4.6
22
and AFGROW (R
4.5
. . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
FEQP normalized
KI
variation with contour intervals:
ix
a = 1.0 in.
(25.4 mm)
37
39
4.8
ical and FEQP methods

4.9
KI
(R = 0.2).
during crack propagation for the analyt. . . . . . . . . . . . . . . . . . . . . . .
40
Comparison of cycle count during crack propagation for the analytical and
FEQP methods (R
= 0.2).
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.10 Comparison of crack propagation direction for the MTS, MERR, and
criterion in the FEQP analysis (R
= 0.2).
41
KII = 0
. . . . . . . . . . . . . . . . . . . .
41
4.11 Partitioning, loading, and boundary conditions used in the VCCT model with
an inset zoomed in at the initial crack.
. . . . . . . . . . . . . . . . . . . . .
43
4.12 Coarse meshing scheme for the VCCT model with an inset zoomed in at the
initial crack. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.13 Example load amplitude curve for the static and cyclic steps with
4.14 Error in stress intensity range,
KI , during crack propagation in VCCT analysis.
4.15 VCCT Convergence plot of estimated fatigue fracture life.
KI
= 0.2).
4.16 Comparison of the normalized

ical and VCCT methods (R
R = 0.2.
. . . . . . . . . .
44
45
48
50
during crack propagation for the analyt. . . . . . . . . . . . . . . . . . . . . . .
51
4.17 Comparison of cycle count during crack propagation for the analytical and
VCCT methods (R
= 0.2).
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
51
4.18 Partitioning, loading, and boundary conditions used in the XFEM model with
. . . . . . . . . . . . . . . . . . . . .
53
4.19 Coarse meshing scheme for the XFEM model with an inset zoomed in at the
initial crack. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.20 Error in stress intensity range,
KI ,
during crack propagation in XFEM-PN
analysis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.21 XFEM-PN Convergence plot of estimated fatigue fracture life.
KI during
(R = 0.2). . .
4.22 Comparison of the normalized

ical and XFEM-PN methods
54
. . . . . . . .
55
56
crack propagation for the analyt. . . . . . . . . . . . . . . . . . .
57
4.23 Comparison of cycle count during crack propagation for the analytical and
XFEM-PN methods (R
= 0.2).
. . . . . . . . . . . . . . . . . . . . . . . . .
58
4.24 Comparison of crack propagation direction from FEQP and XFEM-PN analyses (R
= 0.2).
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
log10 (da/dN )
5.1
Histogram of
5.2
Q-Q plot of
5.3
Monte Carlo simulation of fatigue fracture life with 10,000 samples (R
log10 (da/dN )
from experimental data. . . . . . . . . . . . . . .
from experimental data with a normal PDF.
. . .
= 0.2).
58
64
64
66
5.4
Q-Q plot of the distribution of Monte Carlo simulation compared to a normal

PDF. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.5
Monte Carlo simulation of fatigue fracture life with 10,000 samples with no
correlation between Walker constants (R
= 0.2).
. . . . . . . . . . . . . . . .
5.6
Monte Carlo simulation of fatigue crack growth with 10,000 samples (R
5.7
Monte Carlo simulation of fatigue fracture life estimated with XFEM-PN and
the integrated polynomial method with 500 samples (R
5.8
= 0.2).
= 0.2)
. . . . . . . .
67
67
69
Q-Q plot of the distribution of the XFEM-PN Monte Carlo simulation compared to a normal PDF.
5.9
66
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
69
Q-Q plot of the distribution of the XFEM-PN Monte Carlo simulation compared to the distribution of the Analytical Method simulation. . . . . . . . .
70
6.1
Compact specimen partitioning scheme and coarse mesh. . . . . . . . . . . .
75
6.2
Variation in the normalized
6.3
Comparison of the median normalized
across a compact specimen half model. . . .
KI
during crack propagation in a 3-D
compact specimen for the analytical and XFEM-PN methods (R

6.4
75
= 0.2). .
. .
77
Comparison of cycle count during crack propagation in a 3-D compact specimen for the analytical and XFEM-PN methods (R
= 0.2).
. . . . . . . . . .
77
. . . . . . . . . . . . . . . . . . . . . .
88
A.1
Detailed center crack plate specimen
B.1
Normalized
KI
during crack propagation (R
= 0).
B.2
Normalized
KI
= 0.2).
. . . . . . . . . . . . .
91
B.3
Normalized
KI
= 0.5).
. . . . . . . . . . . . .
92
B.4
Normalized
KI
= 0.67).
. . . . . . . . . . . .
92
B.5
Cycle count during crack propagation (R
= 0).
. . . . . . . . . . . . . . . . .
93
B.6
= 0.2).
. . . . . . . . . . . . . . . .
93
B.7
= 0.5).
. . . . . . . . . . . . . . . .
94
B.8
= 0.67).
B.9
Crack extension direction during propagation (R
. . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
= 0).
91
94
. . . . . . . . . . . .
95
B.10 Crack extension direction during propagation (R
= 0.2).
. . . . . . . . . . .
95
= 0.5).
. . . . . . . . . . .
96
= 0.67).
. . . . . . . . . . .
96
. . . . . . . . . . .
97
B.13 FEQP, deformation mode
ux
with
a = 2.24 in. (R = 0.2).

xi
B.14 FEQP, deformation mode
uy
a = 2.24 in. (R = 0.2).
with
. . . . . . . . . . .
97
B.15 FEQP, principal strain
max
with
a = 2.24 in. (R = 0.2).
. . . . . . . . . . . .
98
B.16 FEQP, principal strain
min
with
a = 2.24 in. (R = 0.2).
. . . . . . . . . . . .
98
B.17 VCCT, deformation mode
ux
with
a = 2.25 in. (R = 0.2).
. . . . . . . . . . .
99
B.18 VCCT, deformation mode
uy
with
a = 2.25 in. (R = 0.2).
. . . . . . . . . . .
99
B.19 VCCT, principal strain
max
with
a = 2.25 in. (R = 0.2).
. . . . . . . . . . .
100
B.20 VCCT, principal strain
min
with
a = 2.25 in. (R = 0.2).
. . . . . . . . . . .
100
B.21 XFEM-PN, deformation mode
ux
with
a = 2.38 in. (R = 0.2).
. . . . . . . .
101
B.22 XFEM-PN, deformation mode
uy
with
a = 2.38 in. (R = 0.2).
. . . . . . . .
101
B.23 XFEM-PN, principal strain
mAX
with
a = 2.38 in. (R = 0.2).
. . . . . . . .
102
B.24 XFEM-PN, principal strain
min
. . . . . . . . .
102
a = 2.38 in. (R = 0.2).
with
C.1
FEQP normalized
KI
a = 0.1 in.
(2.54 mm) 104
C.2
FEQP normalized
KI
a = 1.0 in.
(25.4 mm) 104
C.3
FEQP normalized
mm).
C.4
KI
. . . . . . . . . . . . . . . . . . . . . . . .
VCCT constant residual stabilization ratio,
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
R0
R0 ,
C0 ,
106
during four cycle iterations
labeled once cycle 2253 has been stabilized. . . . . . . . . . . . . . .
Distribution of Walker constant,
105
R0 , during fatigue crack propaga-

with
D.1
105
Typical stabilization ratio for VCCT and XFEM-PN direct cyclic analyses
tion.
C.6
(76.2
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
during fatigue crack propagation.

C.5
a = 3.0 in.
106
with 10,000 samples and a log-normal
PDF. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
107
D.2
m,
D.3
D.4
c3 ,
with 500 samples and a log-normal PDF.
109
D.5
c4 ,
with 500 samples and a normal PDF. . .
109
E.1
Compact specimen dimensions.
. . . . . . . . . . . . . . . . . . . . . . . . .
112
E.2
Compact specimen deformation mode
ux
with
a = 1.46 in. (R = 0.2).
. . . .
112
E.3
uy
with
a = 1.46 in. (R = 0.2).
. . . .
113
with 10,000 samples and a normal PDF.

xii
108
108
uz
E.4
E.5
Compact specimen principal strain
E.6
Crack tip of compact specimen principal strain
E.7
E.8
max
min
xiii
with
with
with
a = 1.46 in. (R = 0.2).
a = 1.46 in. (R = 0.2).

max
with
with
113
. . . . .
114
a = 1.46 in. (R = 0.2).114
a = 1.46 in. (R = 0.2).

min
. . . .
. . . . .
115
a = 1.46 in. (R = 0.2).115
Nomenclature
Ratio of half crack length to half width of plate
Empirical exponents in the Power law
Mean coecient vector
Correlation coecient matrix
bji
Coecient vector
Displacement vector
Stress intensity range
Element length at crack tip
ij
Kronecker delta
max , min
Maximum and minimum strain
Empirical exponent in the BK law
Empirical exponent in the Reeder law
Contour surrounding the crack tip
Linear elastic energy release rate
Gc
Critical linear elastic energy release rate
GIc , GIIc , GIIIc
Mode dependent critical energy release rates
GI , GII , GIII
Mode dependent energy release rates
Gpl
Energy release rate at Paris limit
Gth
Threshold energy release rate
Poisson's ratio
Angular frequency
xiv
Equivalent stress intensity range at
R=0
Nodal residual in Abaqus' direct cyclic solver
Nodal displacements in Abaqus' direct cyclic solver
Level set functions
Far eld normal and shear stress
2,
Variance and covariance matrix
max , min
Maximum and minimum stress
m , a
Mean and alternating stress
x , y , z
Normal stress in the x-, y-, and z-direction
yld , ult
Yield and ultimate stress
xy , yz , xz
Shear stress in the x-y, y-z, and x-z plane
Random error
Half crack length
ac
Critical half crack length at fracture
af
Final measured half crack length at fracture
ai
Initial half crack length
C, m
Coecients for the Paris fatigue crack growth model
C0 , m,
Coecients for the Walker fatigue crack growth model
c1 , c 2
Coecients for Abaqus' onset of crack growth model
C2 , m 2
Coecients for the Forman fatigue crack growth model
c3 , c 4
Coecients for Abaqus' fatigue crack growth model
C3 , m 3
Coecients for Broek and Rice's fatigue crack growth model
CR0 , CRn
Ratio of the largest residual coecient to the time average force for
the constant and periodic terms, respectively
CU0 , CUn
Ratio of the largest correction to the displacement coecients over the

largest displacement coecient for the constant and periodic terms,
respectively
da/dN
Fatigue crack growth rate
xv
Young's modulus
Force component; Dimensionless geometry factor
Fi
Enrichment terms for the nite element shape functions
Shear modulus
Jump function; Heaviside function
I0
Frequency at which residuals of successive iterations are checked
IA
Maximum number of cutbacks allowed in an increment
IR
Starting point for the logarithmic rate of convergence checks
Non-linear elastic energy release rate
Jc
Critical non-linear elastic energy release rate
Stress intensity factor
Kc
Critical stress intensity factor
KIc , KIIc , KIIIc
Mode dependent critical stress intensity factor
KI , KII , KIII
Mode dependent stress intensity factors
Kmax , Kmin
Maximum and minimum stress intensity factor, respectively
Kth
Threshold stress intensity factor
Fatigue cycles
Nf
Fatigue fracture life
Ni (x)
Finite element shape functions
nj
Outward normal vector to contour
Pf
Final load at fracture
Forman plotting variable
Radial distance from crack tip; Ratio for modied Simpson's rule
r,
Axes in a polar coordinate system
R0
Unknown constant residual coecient
re
Enrichment zone radius
Rks Rkc
Periodic residual coecient
xvi
rp
Plane stress plastic zone radius
Set of nodes
u0
Unknown constant displacement coecient
usk uck
Periodic displacement coecient
Stress work density
x, y, z
Axes in a rectangular coordinate system
Half width of plate specimen
Half height of plate specimen
Thickness of plate specimen; time
AFGROW
Air Force Growth fracture mechanics and fatigue crack growth analysis
software tool
CCP
Center Crack Plate
CPD
Crack Propagation Direction
EPFM
Elastic Plastic Fracture Mechanics
FE
Finite Element
FEM
Finite Element Method
FEQP
Finite Elements with Quarter Points
LEFM
Linear Elastic Fracture Mechanics
LSM
Level Set Method
MERR
Maximum Energy Release Rate
MTS
Maximum Tangential Stress
SIF
Stress Intensity Factor
VCCT
Virtual Crack Closure Technique
XFEM
eXtended Finite Element Method
XFEM-PN
eXtended Finite Element Method with Phantom Nodes
CAE
Abaqus' Complete Abaqus Environment le extension
INP
Abaqus' analysis input le extension
ODB
Abaqus' output database le extension
STA
Abaqus' analysis status le extension
xvii
Chapter 1
Introduction
1.1
Overview and Statement of Need
Engineered structures can fail in various ways including yielding, buckling, and brittle fracture. In the presence of a aw or stress raiser, cracks can form and grow to a critical length
where the structure's strength has been reduced to a point where fracture occurs. When this
occurs with the application of a cyclic load, the process is known as fatigue crack growth [1].
The cyclic loads are typically subcritical service loads that would not normally cause failure.
Structures can be designed to minimize stress raisers, but aws are inherent in all materials.
Therefore, design and analysis of structures should take into account cracks that may form
from these aws.
Frequent inspections are used to look for cracks, and if a crack is found the component can
be xed or replaced as necessary. Inspection intervals could then be set at some fraction of
the fatigue fracture life of the structure. Fatigue fracture life is calculated as the number
of cycles required to grow a crack from a minimum detectable size to a critical size when
the structure fails. Analytical fatigue crack growth models are available in the literature for
many generalized geometries including an edge crack, center crack, and embedded penny
crack just to name a few [2]. As geometry, loading and boundary conditions become more
complex, the analytical equations quickly become more complicated and it may be dicult
to account for all eects on crack propagation.
The Finite Element Method (FEM) has been used for decades to assist engineers in analyzing complex, cracked structures. Until recently, cracks had to be modeled as part of the
structure's geometry. As the crack grew, the model would be rebuilt and remeshed, requiring
signicant user interaction or specialized programs. The eXtended Finite Element Method
(XFEM) was developed in 1999 by Belytschko and Black [3], where cracks could be dened
arbitrarily, independent of the mesh. A second iteration of the XFEM was introduced in 2006
by Song, Areias, and Belytschko [4], where Phantom Nodes (XFEM-PN) and a Level Set
Method (LSM) were used to locate the crack. This method naturally lends itself to fatigue
crack growth, where a crack propagates along a solution dependent path, independent of
the mesh. The commercial nite element code, Abaqus, employs XFEM-PN for automated
crack growth in a low-cycle fatigue analysis [5]. A fundamental aspect of this thesis is to
take a critical look at the code's ability to estimate fatigue fracture life.
Fatigue crack growth is highly stochastic, and the actual fatigue fracture life may be dicult
to determine. For this reason, distributions on the life estimate are typically used. Monte
Carlo simulations can be performed to model the complex interactions between the input
distributions such as the variations in geometry, loading, and material properties. This thesis
uses XFEM-PN as a tool for predicting the distribution of fatigue fracture life estimates based
on the variability in the fatigue crack growth models..
1.2
Research Goals and Objectives
The primary goals of this thesis are to investigate the capabilities of the eXtended Finite
Element Method with Phantom Nodes (XFEM-PN) for modeling crack propagation and
estimating fatigue fracture life. Then, use XFEM-PN to quantify the distribution of fatigue
fracture life. The following objectives were met to realize these goals:
1. Fit fatigue crack growth rate models to experimental data.
2. Analyze crack growth using methods established in the literature.
3. Analyze crack growth using XFEM-PN and compare results to the experimental data
and established methods.
4. Combine the variability in the tted crack growth rate models with the XFEM-PN to
determine the distribution on fatigue fracture life estimates.
1.3
Scope of Thesis
Before the XFEM-PN can be used to model fatigue crack growth and estimate fatigue fracture life in everyday structures with complex geometry, complex loading and boundary conditions, material variability, and numerous other external factors, we must rst focus on a
simplied model with a known analytical solution. This thesis focuses on pure Mode I crack
propagation in a Center Cracked Plate (CCP). The CCP oers simplied geometry, simplied boundary and loading conditions, and has been studied since the beginning of fracture
mechanics.
The scope is further limited to Linear Elastic Fracture Mechanics (LEFM) in two dimensions.
CCP crack propagation data was previously collected by Hudson [6] for 7075-T6 aluminum
alloy, thus the analysis will be limited to the same homogeneous, isotropic material and
specimen geometry.
contact.
Monotonic oscillatory loads are analyzed with no consideration for
All nite element analyses were performed in Abaqus 6-13.3, unless otherwise
stated, with no user dened subroutines.
A large number of stochastic variables are available for analyzing the variability of fatigue
fracture life; however, the scope of this thesis is limited to the variability in the constants of
the tted fatigue crack growth models.
1.4
Thesis Outline
This thesis begins in Chapter 2 with a review of fracture mechanics and fatigue crack growth
analysis procedures currently available in the literature. Fatigue crack growth models are t
to experimental data in Chapter 3. Resulting growth models are input into multiple fatigue
crack growth analysis tools in Chapter 4 and compared to XFEM-PN. Chapter 5 uses the
Monte Carlo method to analyze the variability in the life estimates. Finally, conclusions are
made in Chapter 6 with a brief look at future work in three-dimensional fatigue crack growth.
Additional data, results, and MATLAB and Python code can be found in the Appendices.
Chapter 2
Literature Review
The following literature review has been provided to give a brief background into fracture
mechanics and fatigue crack growth laws that will be used in the remainder of this thesis.
Established techniques in fatigue fracture analysis are explored, and examples of their use
as they apply to this thesis are provided. Finally, a look at the Monte Carlo method, its
assumptions, and applications to engineering analyses are given.
2.1
Material Models
The following sections give a brief background on fracture mechanics, fatigue crack growth
models, methods used by Abaqus for fatigue crack growth analyses. These methods include
approximating the critical mixed-mode energy release rates and rules for determining crack
propagation directions which are used in the Abaqus analysis presented in Chapter 4.
2.1.1
Fracture Mechanics
Fracture mechanics is broadly divided into two types: Linear Elastic Fracture Mechanics
(LEFM), and Elastic Plastic Fracture Mechanics (EPFM). LEFM assumes small deformations and minimal yielding at the crack tip, while EPFM can account for large deformations
and plastic eects [1].
In LEFM, the Stress Intensity Factor (SIF),
and is calculated with Eq. (2.1) where
nominal stress, and
K,
is a measure of the stress eld at a crack tip
is a dimensionless geometry factor,
is the crack length.
is a remote
This equation was developed from a theory of
elasticity solution of the stress eld around a sharp notch [7]. From the elastic solution it
was concluded that the stress eld is proportional to
1/ r
where
is the radial distance
from the crack tip. This stress proportionality results in a stress singularity at the crack tip
where
r 0.

K = F a
(2.1)
SIFs are divided into three modes based on the displacement at the crack tip, as shown in
Fig. 2.1. Mode I, the opening mode, is caused by a displacement perpendicular to the crack
plane, which is typically a result of tensile stresses. Mode I is primarily responsible for crack
growth. Displacements perpendicular to the crack tip edge from in-plane shear stresses cause
the sliding mode, Mode II. Out-of-plane shear stresses result in displacements parallel to the
crack tip edge and a tearing mode, Mode III. To dierentiate between modes, subscripts
II ,
and
III
are appended to
K.
I,
[8]
An energy approach to fracture mechanics was rst explored by Grith [9] in 1921 in a
paper analyzing the brittle fracture of glass. His approach states that the energy required
t (da),
to grow a crack by an area,
is equal to the change in potential energy,
dU ,
over the
change in crack area as shown in Eq. (2.2) [7]. In LEFM, the potential energy is comprised
of the internal strain energy.
G =
1 dU
t da
(2.2)
His work was extended to metals by Irwin [10] in 1957, who showed that the strain energy
release rate,
G,
is related to
in LEFM using the Young's modulus,
E,
and Poisson's ratio,
in Eq. (2.3).
G=
E=
K2
E
(2.3)
for plane stress

(2.4)
E
1 2
for plane strain
As described by Anderson [11], Rice [12] proposed the use of a line integral around the crack
J -integral, to characterize the strain energy release rate for EPFM when large
J , as calculated in Eq. (2.5), assumes an idealized non-linear
material response. In this equation, W is the strain energy density, is a contour
tip, called the
deformations are present.

elastic
Fig. 2.1:
Three modes of displacement at a crack tip.
surrounding the crack tip,

and
is the outward normal to this contour,
is the displacement,
is the Kronecker delta.

ui
J=
W 1j ij
nj d
x1
Shih and Asaro [13] related the
(2.5)
J -integral to K for homogenous, isotropic materials

E is calculated with Eq. (2.4) and G is the
LEFM with the following equation, where
using
shear
modulus.

1
1 2
2
K
KI2 + KII
+
2G III
E
J =G=
Fracture is assumed to occur when
(2.6)
reaches the critical stress intensity factor,
Kc , which is
dependent upon the material and specimen geometry. Should the specimen be large enough
that a state of plane strain exists,
Kc
will approach a minimum value called the plain strain
critical stress intensity factor, or fracture toughness.

extended to the critical energy release rate,
2.1.2
Gc .
The critical value concept can be
[7]
Mixed-Mode Critical Energy Release Rate
Fracture often occurs under a combination of modes, thus an equivalent critical value is
needed to include the contributions from each mode. The total energy release rate is calculated as the sum of the energy released from each mode,
G = GI + GII + GIII ,
but this is not
the case with the critical values. A universally accepted method for determining a mixed
Gc does not exist, thus an appropriate method must be selected on a case-by-case basis.
following mixed-mode Gc calculations are available in Abaqus 6-13.3.
mode
The
In 1984, Whitcomb [14] used the superposition of stresses to develop a power law relationship
describing the mixed mode delamination of composites during buckling assuming that Mode
I and II are primarily responsible for delamination. One of the results from his report was
Eq. (2.7) where
and
are empirically determined constants.

Gc =
GI
GIc

+
GII
GIIc

(2.7)
While testing glass/epoxy composites in 1996, Benzeggagh and Kaine [15] observed that
Gc
increased as fracture transitioned from pure Mode I to pure Mode II. They developed the
mixed Mode I-II, semi-empirical criterion in Eq. (2.8) as a combination of the critical values
based on the fraction of
GII
present in the analysis. The exponent,
is t to experimental
data.

Gc = GIc + (GIIc GIc )
6
GII
GI + GII

(2.8)
This was extended to three dimensions in 2006 by Reeder [16] at NASA. Eq. (2.9) combines
all three fracture modes and
is an empirical constant.
Reeder also concluded that Eq.
(2.9) has no physical signicance in two dimensions and thus should not be used in those
cases.

GII + GIII
GI + GII + GIII
1

+ (GIIIc GIIc )
2.1.3
GIII
GII + GIII

1
GII + GIII
GI + GII + GIII
(2.9)
Fatigue Crack Growth Models
Besides static loading, crack growth can occur when a subcritical load is repetitively applied.
In 1961 Paris, Gomez, and Anderson hypothesized that the fatigue crack growth rate,
da/dN ,
K = Kmax Kmin [17]. On a log-log plot, it is

now known that da/dN has a sigmoidal relation with K as seen Fig. 2.2 and any textbook
was related to the stress intensity range,
on fracture mechanics [8, 11].
The crack growth curve can be divided into three regions of generalized behavior. In the
threshold region, crack growth is slow as
Kth ,
asymptotically approaches the threshold value,
where crack growth may not occur. The slope of the crack growth curve is approxi-
mately linear in the intermediate, or Paris, region. And the unstable region is characterized
by rapid, unstable crack growth where
Kmax
asymptotically approaches
Kc
and fracture is
imminent. Numerous crack growth models have been developed to describe the relationship
between
da/dN
and
within these regions. A select few are presented in this thesis.
In 1963, Paris and Erdogan [18] approximated the intermediate crack growth region with a
power law relationship known today as the Paris equation, Eq. (2.10), where
and
are
empirically determined material constants. One shortcoming of this relationship is that it

applies to a single stress ratio,
R = min /max ,
where stresses are dened far from the crack
tip.
da
= C (K)m
dN
Walker [19] reported on the eects of the
in 1970, concluding that increasing
developed to address the eect of
R-ratio
(2.10)
on crack propagation for aluminum alloys
resulted in an increased growth rate. Eq. (2.11) was
R-ratio
on crack growth rate. Constants
C0 , m,
and
are
empirically determined from experimental data.
da
C0
=
(K)m
dN
(1 R)m(1)
7
(2.11)
da/dN (length/cycle)
Threshold
Region
Paris
Region
Kmax=Kc
Unstable
Region
Kth
K (Stresslength1/2)
Fig. 2.2:
Typical crack growth curve on log-log axes.
Paris and Walker equations work well in the Paris region of crack growth, but do not address
the asymptotic behavior in the unstable region. To address this behavior, Forman, Kearney,
and Engle proposed a model in 1967 [20], known as the Forman equation, Eq. (2.12). As
Kmax
Kc , the denominator approaches zero and da/dN tends toward

C2 and m2 are empirically determined through experimental testing.
approaches
Constants
C2 (K)m2
C2 (K)m2
da
=
=
dN
(1 R) Kc K
(1 R) (Kc Kmax )
innity.
(2.12)
To incorporate the curvature in the threshold region, Broek and Rice [21] developed Eq.
(2.13) for rail steels, as presented by Farris, Keer, and Steele [22].
adjusts the growth towards zero as
equation to innity as
Kmax
approaches
The bracketed term
Kth , and the denominator forces the

C3 and m3 are empirically determined
approaches
Kc .
Again,
constants.

(K)m3 1
da
2
2
= C3 (K) Kth
dN
Kc Kmax
Eqs.
(2.12) and (2.13) require data over a wide range of
K 's.
(2.13)
Since this is sometimes
dicult to obtain, Dowling [8] recommends using either Eq. (2.10) or (2.11) for growth in
the Paris region and bounding it by
range,
Kth
and
Kmax = Kc , or in terms of the stress intensity
Kc = Kc (1 R).
The previously presented equations do not address the many factors that inuence crack
growth including load frequency, load sequencing, temperature, and environmental factors
just to name a few. The following examples are given to show the narrow window in which
Eqs.
(2.10) through (2.13) apply.
Hartman and Schijve [23] showed that an increase in
loading frequency produced a decrease in growth rates in aluminum alloys. Yokobori and
Kiyoshi [24] followed this research by proposing the addition of a frequency term to Eq. (2.10)
to address these aects. Hartman and Schijve's report [23] also found signicant dierences
in the growth rates depending on the environments in which the tests were performed. Their
report focused on wet and dry oxygen, air, and argon environments. In 1973, Yokobori and
Aizawa [25] found that increasing the temperature of aluminum alloys increased the crack
growth rate.
Finally, sequencing of the loads can play a large role in crack growth rates.
is when overloading occurs prior to normal service loads.
An example
Overloading the crack can lead
to large plastic deformation at the crack tip. The plasticly deformed region is surrounded
by undeformed material that, upon removal of the applied overload, elasticly returns to its
original conguration. The undeformed region places a compressive stress on the deformed
region around the crack tip. The compressive stress around the crack tip results in crack
growth retardation [8].
2.1.4
Onset of Crack Growth
Abaqus 6-13.3 requires a form of Eq.
(2.14) to dene the fatigue crack growth criterion.
This section briey address it's applicability to the current analysis.

Murri, Salpekar, and O'Brien [26] proposed a method to determine the maximum
GI
that
can be applied for the onset of delamination of tapered, unidirectional laminate. Eq. (2.14)
proposes that the
N,
and constants
GImax required for the onset of delamination is a function of the cycles,

c and d. Analyses using this onset model are covered by Murri, O'Brien,
and Rousseau [27] and a mixed-mode analysis is presented by Giannis [28].

According to ASTM Standard D6115-97 [29], this method is only known to apply to unidirectional composites. Section 4.5.1 of this thesis discusses setting the constants
and
to
zero to remove Eq. (2.14) from the analysis.
GImax = cN d
2.1.5
(2.14)
Crack Extension Direction
Many models are available for predicting the direction of crack extension based on either
stress, strain, energy, or any combination of these. Appropriate models should be selected
based on material and loading conditions. The three models presented here were selected
since they are popular in the literature and are available in Abaqus 6-13.3 without the need
for a user dened subroutine.
Erdogan and Sih [30] proposed the Maximum Tangential Stress (MTS) which states that
crack extension will occur radially from the crack tip and perpendicular to the maximum
applied tensile load.
These two criterion are met when the tangential stress,
is maxi-
mized and the shear stress,
r ,
is zero. Erdogan and Sih developed parametric equations
to determine the propagation angle,
in
KI -KII
space. Setting
r = 0
and
/ = 0,
the parametric equations were solved in the Abaqus Theory Manual [31] and are provided
in Eq. (2.15). Validation experiments were performed on Plexiglas sheets and Erdogan and
Sih reported good agreement with the parametric equations..
= cos1
!
p
2
2
3KII
+ KI4 + 8KI2 KII
2
KI2 + 9KII
(2.15)
In 1974, Hussain, Pu, and Underwood [32] developed the Maximum Energy Release Rate
(MERR) criterion which was based on Grith and Irwin's work. MERR states that crack
extension will occur at an angle,
that maximizes
in Eq. (2.6) for a mixed Mode I-II
condition. Equations are provided in [32] that compute

assume
KIII = 0.
KI = KI (), KII = KII (),
and
Hussain et al. state that the MERR model compared well to experimental
data collected for crack growth in steel foil.

Cotterell and Rice [33] suggested in 1980 that crack extension occurs in a direction where
KII = 0 for isotropic, homogeneous materials. Cotterell and Rice also proposed that the MTS
and MERR models' solutions meet the KII = 0 criterion once the crack has extended. Curved
and kinked crack models were analyzed in their article, but no experimental verication was
reported.
Each crack extension direction criterion will give slightly dierent results. Analysis of crack
propagation in Abaqus should be run with each criterion described in this section and compared to experimental data for selection of the most appropriate one.
2.2
Methods of Fracture Analysis
A fracture analysis typically starts with an initial crack size, or crack initiation criteria based
on stress or strain, and propagates the crack until a critical value is reached such as
Gc .
Kc
or
Growth rates are calculated using a fatigue crack growth model such as the Walker
equation, Eq. (2.11). Several methods have been developed in the literature and a select
few are presented in the following sections.
2.2.1
Analytical Methods
Integrating a fatigue crack growth law such as Eq. (2.10) is the most basic approach explored
in this thesis for solving for the fatigue fracture life of a specimen. The integral approach
assumes self-similar Mode I crack propagation and is used a validation tool in this thesis.
Programs have been developed to aid in these calculations such as AFGROW, created by
The Air Force Research Laboratory and LexTech, Inc. [34], and NASGRO, created by NASA
Johnson Space Center and Southwest Research Institute
10
[35]. Programs like these access
a database of equations for SIFs, and apply a load spectrum and fatigue crack growth law to
analyze crack growth and estimate fatigue fracture life. AFGROW is used in this thesis for
the crack growth analysis. While these programs work well for the catalog of design situations, a nite element analysis will likely be required for more complex geometry, boundary
conditions, and loading conditions.
2.2.2
Finite Element Methods
Finite Element Methods (FEMs) are used to provide approximate solutions to partial differential equations in the form of either boundary-value or initial-value problems. The FEM
was rst developed as early as 1851 and gained wide-spread use with the advent of the
digital computer in the 1960s [36]. The FEM has been applied to problems in stress and
strain analysis, wave propagation, heat transfer, electrical and magnetic elds, and uid
ow. Currently, nite element analysis software is readily available in commercial codes like
Abaqus and ANSYS. These codes are continually updated with the most recent advances in
the literature such as the eXtended Finite Element Method for fracture mechanics.
Early Applications of the FEM.
Some of the rst applications of FEM to fracture
mechanics were reported in 1969 by Watwood [37] for a center cracked plate, and Anderson,
Garron, and Ruggles [38] for fracture of a solid propellant rocket motor cartridge.
These
early applications of the FEM to problems in fracture mechanics were followed in 1970 by
Chan, Tuba, and Wilson [39] who analyzed compact and rotating test specimens. All authors
noted that the stress singularity at the crack tip was unattainable with linear elements, and
thus their analyses quickly led to Grith's energy release rate or to the contour integral for
and
K.
All reports showed good agreement between energy release rates, contour integrals,
and analytical methods.
Specialized Crack Tip Elements.
In 1970, in an eort to address the issues with mod-
eling the stress singularity, Byskov [40] developed a triangular crack tip element containing
necessary shape functions to accurately model the singularity. This approach was expanded
to circular elements by Wilson [41] in 1973 and rectangular elements in 1974 by Hardy [42].
The circular element has had many versions including a mixed-mode element by Holston [43]
in 1976 and a bending specic version developed by Jiang and Cheung [44] in 1995.
Finite Elements with Quarter Points.
To remove the need for specialized elements,
Henshell and Shaw [45] published a paper in 1975 which concluded that moving the mid-side
nodes of quadratic elements adjacent to the crack tip to their quarter point location adequately models the
1/ r
stress singularity. Barsoum [46] simultaneously published similar
ndings without knowledge of Henshell and Shaw's work.
A schematic of this concept is
shown in Fig. 2.3. The following year, Barsoum [47] showed that all nodes of the elements
at the crack tip should be collapsed to a single node for linear elastic analyses, and that this
11
method could be extended to model the
1/r
singularity for a perfectly plastic analysis when
the collapsed nodes remain independent.

To analyze the fatigue crack growth using Finite Elements with Quarter Points (FEQP),
nite element models must be built at incremental crack lengths since explicit crack propagation is not possible. Typically a concentric meshing scheme is used around the crack, as
seen in Fig. 2.3. Hjfeldt and stervig [48] illustrated this procedure in 1986 for cracks in
shoulder llets of a shaft. In general this process proceeds as follows:
1. Build and mesh a nite element model including the crack in the geometry.
2. Apply quasi-static load from minimum to maximum value.
3. Determine
K 's
or
G 's.
4. Determine crack extension direction.

5. Calculate incremental crack growth length.
6. Determine the cycles required to grow the crack this incremental length based on a
fatigue crack growth model.
7. Add cycles to the previous cycle count.
8. Rebuild the model and mesh with new crack geometry.
9. Repeat steps 2 through 8 until critical values are reached.
Recent publications show that this process has remained basically unchanged, including Alegre and Cuesta's [49] 2010 paper where crack propagation paths where compared for various
initial crack angles. Both [48] and [49] observed good comparison between experimental and
nite element results.
Virtual Crack Closure Technique.
Rybicki and Kanninen [50] developed the modied
crack closure integral in 1977 which was based on Irwin's 1957 assertion that the energy
needed to grow a crack is equal to the work required to close a crack of the same length.
Eq. (2.16) was developed in the paper for Mode I and Mode II energy release rates for a
unit thickness. This was later extended to Mode III by Shivakumar, Tan, and Newman [51].
Referring to Fig. 2.4 to dene variables;
and
from separating,
u_,ab
F_,cd
is the incremental crack length over which

in front of the crack tip.
are force components required to keep nodes
are the displacement components between nodes
is being calculated, and
and
is the element length
These equations apply only to uniform meshes where
equations are provided in [50] for
a 6= l.
12
b, a
a = l;
Fig. 2.3:
Schematic mesh of elements around a crack with mid-side nodes moved to the
quarter point location.
GI =
1
Fy,cd uy,ab
a0 2a
(2.16a)
GII =
1
Fx,cd ux,ab
a0 2a
(2.16b)
1
Fz,cd uz,ab
2a
(2.16c)
GIII =
lim
lim
lim
a0
Today, Rybicki and Kanninen's method is known as the Virtual Crack Closure Technique
(VCCT) and the technique has been used in a myriad of publications for static and cyclic
analyses and three dimensional congurations that may be too dicult to mesh using FEQP.
In general, the process used in the literature to model VCCT cyclic crack growth is similar
to that for FEQP as shown in Mabson et al. [52] and Hosseini-Toudeshky, Ghaari, and
Mohammadi [53].
Mabson et al.
analyzed the fatigue fracture life of double cantilevered
beam specimens while Hosseini-Toudeshky et al. looked at the crack growth path in curved
aluminum panels with stieners and composite repair patches. Both studies observed good
agreement between experimental data and VCCT results.
eXtended Finite Element Method.
The nite element methods discussed so far re-
quire that the crack be incorporated into the geometry so that elements align with the crack
boundary. The eXtended Finite Element Method (XFEM) was developed in 1999 by Belytschko and Black [3] as a method where cracks could be dened arbitrarily within the mesh.
Their method was based on the partition of unity presented by Melenk and Babuka [54] in
1996 which broadly states that the sum of all shape functions must be one, thus allowing
additional terms to be added to the shape functions with
a priori
knowledge of the solution.
Belytschko and Black suggested using enriched shape functions shown in Eq. (2.17) for the
displacement eld,
u. N (x) are standard shape functions, b and c are nodal coecient, (r, )
13
Fig. 2.4:
Schematic nite element representation of the VCCT method.
S1 is a set containing all nodes

F , are only activated in the element containing
is a polar coordinate system originating at the crack tip, and

in the FE model. Crack tip enrichment terms,
the crack tip, represented by node set
function enrichment term,
H (x),
S3 .
Mos, Dolbow, and Belytschko [55] added the jump
which is used to represent the discontinuous displacement
eld in elements that have been completely bisected by a crack.
S2
is the node set where the
jump function is active. Fig. 2.5 illustrates where crack tip and jump function enrichment
terms would be used. SIFs calculated with contour integrals around the crack tip in [3, 55]
were in excellent agreement with analytical values.
"
u (x) =
Ni (x)ui +
iS1
{Fi (r, )}4i=1
Nj (x)bj H (x) +
jS2
Nk (x)
4
X
cnk Fn (r, )
(2.17a)
n=1
kS3
r cos , r sin , r sin sin , r cos sin

2
2
2
2
(
1
H (x) =
1
y>0
for y < 0
for
(2.17b)
(2.17c)
A Level Set Method (LSM) was developed by Osher and Sethian [56] in 1988 to locate moving
boundaries such as grain boundaries and inclusions in material models. Stolarska, Chopp,
Mos, and Belytschko [57] incorporated the LSM concept into XFEM in 2001.
The LSM
implementation in XFEM uses two nearly orthogonal signed distance functions, (x, t) and
(x, t), to locate the crack front. (x, t) = 0 locates the crack face and i (x, t) = 0 denes
th
the i
crack tip. Nodal values for and are incorporated into the mesh with the standard
shape functions with Eq. (2.18). Level set regions used to dene a crack are illustrated in
Fig. 2.6. As the crack propagates through the mesh, values of
14
are stored in the output le
Fig. 2.5:
XFEM nodes with enrichment and jump functions.
since the previously formed crack face is assumed not to change. Stolarska et al. modied
Eq. (2.17c) to include the LSM values resulting in Eq. (2.19).
i (x, t) =
ij (t) Nj (x)
(2.18a)
j (t) Nj (x)
(2.18b)
jS2
(x, t) =
X
jS2
(
1
H (x) = H ( (x, t)) =
1
(x, t) > 0
for (x, t) < 0
for
(2.19)
XFEM has become extremely popular since its inception in 1999. Most crack propagation
analyses follow a similar processes used by FEQP and VCCT; however, remeshing is typically
not an issue. The basic procedure is as follows:
1. Build and mesh a nite element model.

2. Dene the crack location by nodal level set values.
Fig. 2.6:
Level Set Method used for XFEM.
15
3. Apply quasi-static load from minimum to maximum value.

4. Determine
K 's
or
G 's.
5. Determine crack extension direction.

6. Calculate incremental crack growth length.
7. Determine the cycles required to grow the crack this incremental length based on a
fatigue crack growth model.
8. Add cycles to the previous cycle count.
9. Redene crack location by adjusting nodal level set values without remeshing.
10. Repeat steps 3 through 9 until critical values are reached.
XFEM was incorporated into an Abaqus user dened subroutine by Shi, Chopp, Lua, Sukumar, and Belytschko [58] in 2010 to analyze cyclic crack growth in standard and modied
compact specimens, as well as a complex helicopter part. Shi et al.'s subroutine used a
ied
mod-
VCCT method for calculating the energy release rates at the crack tip. A signicant
amount of published XFEM fatigue crack growth articles deal with complex geometry that
would be time consuming or impossible to mesh using FEQP or VCCT methods.
These
include Singh, Mishra, Bhattacharya, and Pati's [59] analysis of the eects on crack growth
by randomly adding inclusions, holes, and minor cracks to two dimensional center cracked
and edge cracked plates. Or Pathak, Singh, and Singh's [60] analysis of embedded and surface elliptical cracks, lens cracks, and arbitrary shaped cracks of various aspect ratios and
inclinations. Each of these articles noted good agreement with analytical and experimental
values when available.
16
eXtended Finite Element Method with Phantom Nodes.
Song, Areias, and Be-
lytschko [4] developed an adaptation of XFEM to model arbitrary crack locations within a
mesh using phantom nodes in 2006.
Drawing conclusions from previous work by Hansbo
and Hansbo [61] and Belytschko and Black [3], the eXtended Finite Element Method with
Phantom Nodes (XFEM-PN) represents the discontinuous displacement eld with superposed elements each containing phantom nodes and phantom degrees of freedom. Prior to
fracture, superposed elements move together, acting as a single element. Once bisected by a
crack, real nodes and phantom nodes move independently of each other, schematically shown
in Fig. 2.7. Eq. (2.20) represents the displacement eld,
as the sum of the displacements
S1 and S2
are the sets of nodes corresponding to each element and
is the Heaviside step function.
XFEM-PN has the advantage over standard XFEM in
from elements 1 and 2.
u,
that no crack tip enrichment terms are needed to model the stress singularity, thus it can be
integrated directly into most FE codes. Abaqus 6-13.3 combines phantom nodes and LSM
for crack propagation analyses [5].
u (x, t) =
(1)
ui (t) Ni (x)H (f (x)) +
(2)
ui (t) Ni (x)H (f (x))
(2.20a)
iS2
iS1
(
1
H(x) =
0
x>0
for x 0
for
(2.20b)
Steps followed to model fatigue crack growth with XFEM-PN would be the same as those
with the standard XFEM. Published literature on XFEM-PN for fatigue crack propagation
was not as abundant compared to XFEM when this thesis was written. Kucharczyk, Sharaf,
and Mnstermann [62] successfully modeled the reduction in load capacity due to crack
growth through steel grains; however this paper used traction separation laws with cohesive
elements, and not the fatigue crack growth models such as the Paris or Walker equations.
Fig. 2.7:
Schematic mesh with phantom nodes.
17
2.3
Stochastic Analyses in Fracture Mechanics
Fatigue crack growth is inherently stochastic and dicult to predict.
Virkler, Hillberry,
and Goel [63] published a paper in 1979 which analyzed the statistical distribution of crack
propagation in 2024-T3 aluminum center cracked plates.
Tests were conducted on sixty-
eight identical specimens with constant amplitude loading.
Experimental results showed
large amounts of variability in crack growth rates not only between samples but also within
each sample. The variability within the each sample was attributed to the non-homogeneity
of the material at the microscopic level, thus growth rate Eqs. (2.10)-(2.13) only apply in
an average sense macroscopically.
One way to model a stochastic process is to use a Monte Carlo method where constants are
randomly selected from a Probability Density Function (PDF). A normal PDF is based on
2
the constant's mean value, , and variance, . The analysis is performed using samples
from the distribution of the random variables. Repeated analyses are performed continuing
to draw samples from the distributions of the random variables until a mean outcome from
the analysis is realized. The fundamental principles of the Monte Carlo method are the law
of large numbers, and the central limit theorem [64]. The law of large numbers states that
the mean value of a sample approaches the mean value of the population as the number
of samples approach innity. The central limit theorem is used to estimate the uncertainty
in the mean value, thus relating the standard deviation of the sample,
deviation of the population,
s,
to the standard
Annis [65] performed a Monte Carlo analysis and compared results to experimental data
presented in Virkler's report [63]. If random numbers were selected for
and
m in the Paris
equation, the Monte Carlo simulation resulted in a variation more than seven times greater
than the variation in the experimental data. To x the error, he proposed that
picked independently of
cannot be
thus a bivariate normal distribution must be used. Multivariate
normal distributions could be used if more than two constants are correlated. Applying this
concept; Annis reduced the error to approximately 1%.
18
Chapter 3
Crack Growth Rates for 7075-T6
Aluminum Alloy
In 1969, C. Michael Hudson published a report for NASA on the eects of the
R-ratio
on
fatigue crack growth rates for two types of aluminum alloys, 7075-T6 and 2024-T3 [6]. The
data collected for the 7075-T6 aluminum alloy is the focus of this thesis.
The following
sections discuss the experimental data that was collected and the tests performed, three
fatigue crack growth models t to the collected data, and a selection of a model to be used
in the remainder of the thesis.
3.1
Experimental Data
Hudson [6] performed three material characterization tests on the 7075-T6 aluminum alloy
specimen to determine the properties needed for a fatigue crack growth analysis: (1) tensile
test for elastic properties, (2) fatigue crack growth test for a crack growth rate model, and
(3) residual static strength test for the critical stress intensity factor,
Kc .
The aluminum alloy specimen measures 12 in. (305 mm) wide by 35 in. (889 mm) tall and
0.09 in. (2.29 mm) thick. It was loaded by ve pins on both the top and bottom as seen in
Fig. 3.1. An initial stress raiser 0.1 in. (2.54 mm) long was made in the center of the plate
with an electrical discharge machine (EDM). A detailed drawing can be found in Fig. A.1
of Appendix A.
19
Fig. 3.1:
Center Cracked Plate Specimen
Elastic properties were determined from a standard tensile test and are listed in Table A.1.
Experimental data from the tensile tests was not provided in Hudson's report.
Fatigue crack growth tests were run on multiple axial-load fatigue-testing machines which
applied monotonic oscillatory loads with a mean stress,
m , and alternating stress, a , levels
listed in Table A.2 of Appendix A. As the crack propagated through the specimens, the cycle
count was recorded at incremental crack lengths. Two tests were performed at each load level
and the reported results are the average results from both tests.
It should be noted that
the incremental crack lengths provided in Hudson's paper are coarse compared to current
a 0.05 in.
with an R 0
ASTM Standard E647. Section 8.8 of ASTM Standard E647 [66] suggests a
while Hudson's measurements ranged from
a = 0.1
were analyzed in this thesis since Hudson found that

specimens loaded with
0.2 in. Only

R < 0 had the
to
tests
same growth rate as
R = 0.
Approximately one-third of the tests were stopped prior to failure in the fatigue crack growth
test, and these specimens were used in the residual static strength test. It was not specied in
Hudson's report which specimens were stopped prior to failure in fatigue crack growth test.
For the static strength test, an increasing axial load was applied to the specimens until rapid
unstable crack growth was visually observed with a camera. When failure was observed, the
nal load,
Pf ,
and nal crack length,
in Table A.3 were used to calculate
af ,
Kc
were recorded. Resulting
Pf
and
af
values listed
for the specimen. Visually determining when failure
occurs does not conform to current ASTM Standard E399. Section 9 of ASTM Standard
E399 [67] suggests the use of a load versus displacement, or
P -v ,
curve to determine when
failure occured. The visual method may result in an overestimated nal load and nal crack
length, which would lead to an overestimated
Kc .
The eects of an overestimated
Kc
will
be mentioned again in Section 4.2 when analyzing the estimated critical crack lengths,
from analytical equations.
20
ac ,
3.2
Fatigue Crack Growth Models
Many crack growth models have been developed to describe the relationship between
and
da/dN
within the three regions of generalized crack growth behavior seen in Fig. 2.2. Three
models explored in this thesis are the Paris, Walker, and Forman crack growth models.
3.2.1
Paris Crack Growth Model
The rst crack growth model investigated in this thesis is the Paris equation, Eq.
(3.1).
Crack length and cycles provided in Table A.2 were used to determine the crack growth
rate,
da/dN .
A three-point incremental polynomial method, as described in Appendix XI of
ASTM Standard E647 [66], was selected to calculate
da/dN
since it gave smoother results
than the secant method. Using more than three points to dene the polynomial removed
too much data without any noticeable gains. The stress intensity range,
with Eq. (3.2) where
is the stress range calculated as
length. For the dimensionless geometry factor,
F,
= 2a
Eq. (3.3) from Tada [2] was used since
this form is accurate to within 1% of experimental data for any

where
is the half-height and
K , was calculated
a is the half-crack
and
= a/b
and
h/b 1.5
is the half-width of the specimen. For this specimen, the
half-height was measured arbitrarily
1.5 in.
below the centerline of loading holes as it was
assumed that the load would be uniformly distributed at this point. Measuring down to the
horizontal centerline of the plate resulted in an
h/b = 2.0.
da
= C (K)m
dN
(3.1)
K = F () a
F () =
Constants
and
(3.2)
1 0.5 + 0.3262
(3.3)
were found using a least squares linear regression in MATLAB after
transforming the data to log-log space. Constants for models t to each
R-ratio are listed in
Table 3.1 and a plot is provided in Fig. 3.2 showing the data points and tted models.
The Paris equation ts well for each
are relatively small.
R-ratios
R-ratio
since all Root Mean Square Errors (RMSE)
of 0.2 and 0.33 had the smallest error since all data appeared
be recorded within the intermediate region of crack growth where growth rates are linear in
log-log space. Models with
R = 0,
0.5, and 0.8 have a RMSE
> 0.1
which is a result of the
Kmax approaches Kc . The largest

R = 0.8 which may also be attributed to the smaller data set. Error in the
6
model with R = 0.67 is attributed to the slower growth rates where da/dN < 10
in./cycle.
increased growth rates near the end of each data set where
RMSE was with
The slower growth rates could be a result of the variability in growth rates or growth rates
could be approaching the threshold region of crack growth, but no other data points were
available in this region for comparison. The Walker equation was used next, combining all
Paris models into one, using the entire data set, and addressing the
21
R-ratio
eect.
Table 3.1:
Paris Equation Coecients
R-ratio
RMSE
(MPa mm)
9.83 1014
3.64
0.114
4.28 1010
4.72 1015
4.13
0.083
3.79 10
13
2.86 10
3.59
0.095
5.47 109
3.40 1013
3.64
0.128
14
4.14
0.108
15
4.68
0.184
in./cycle
mm/cycle
(ksi in.)
1.56 109
0.2
0.33
0.5
3.36 10
0.67
2.50 10
0.8
3.55 10
3.86 10
1/2
10
10
K (MPamm
10
10
3
10
10
Ratio
_____
5
10
0
0.2
0.33
0.5
0.7
0.8
10
10
10
10
10
10
da/dN (mm/cycle)
da/dN (in./cycle)
10
10
2
10
K (ksiin1/2)
Fig. 3.2:
3.2.2
Paris equation t to Hudson's experimental data.
Walker Crack Growth Model
R-ratios to a single line when

K is dened as the equivalent stress intensity
K values calculated previously were used for the
The Walker equation, Eq. (3.4), collapses data from dierent
da/dN versus K in Eq.

at R = 0. The same da/dN
plotting
(3.5).
range
and
Walker equation. A least squares multiple linear regression was used in MATLAB to t the
model after transforming the data to log-log space. The coecients and a plot of the results
are shown in Table 3.2 and Fig. 3.3, respectively.
22
da
C0
=
(K)m
m(1)
dN
(1 R)
K =
(3.4)
K
(1 R)1
(3.5)
The Walker model ts well to all data with an RMSE value comparable to the Paris models.
R = 0, 0.5, and 0.8

R = 0.67 experiment. The
Error was attributed to the increased growth rates in experiments with

as they approached fracture and decreased growth rates in the
Forman equation was used next to address the increased growth rates in the unstable region.
Table 3.2:
Walker Equation Coecients
C0
in./cycle
mm/cycle
(ksi in.)
8.63 1010
(MPa mm)
2.66 1014
Fig. 3.3:
3.2.3
RMSE
3.84
0.564
0.130
Walker equation t to Hudson's experimental data.
Forman Crack Growth Model
The Forman equation, Eq. (3.6), was used to model the asymptotic increases in the crack
growth rate as
Kmax
approached
Kc .
When
Q
23
in Eq. (3.7) was plotted against
K ,
data
from dierent
R-ratios
collapse to a single line in log-log space and the curvature in the
unstable region was removed.
da
C2 (K)m2
=
dN
(1 R) Kc K
Q=
(3.6)
da
[(1 R) Kc K]
dN
(3.7)
Data from the residual static strength test was used to determine the critical stress intensity
factor for the material and CCP geometry. Eq. (2.1) was rewritten in the form of Eq. (3.8)
where
is the plate thickness and all other variables were previously dened.
this equation to the data in Table A.3 and averaging the results gives
(2500
MPa
Applying
Kc = 72.0 ksi
in.
mm).
Kc = F
a P
f
f
af
b 2bt
Similar to the Walker equation, this model addresses the
(3.8)
R-ratio eects, thus only one model
was needed to describe the entire data set. Following the same procedure as before, and using
a least square linear regression on the transformed data, the coecients in Table 3.3 and the
plot Fig. 3.4 were developed.
Table 3.3:
Forman Equation Coecients
C2
in./cycle
m2 1
(ksi in.)
3.39 107
mm/cycle
m2
RMSE
3.14
0.119
m2 1
(MPa mm)
1.24 1010
24
10
10
K (MPamm1/2)
10
Q (in./cycle ksiin.1/2)
10
RRatio
_______
10
10
0
0.2
0.33
0.5
0.7
0.8
10
10
10
Q (mm/cycle MPamm1/2)
10
10
10
10
10
10
1/2
K (ksiin
Fig. 3.4:
Forman equation t to Hudson's experimental data.
The Forman equation ts well to all data with an RMSE similar to the Walker and Paris
models. Error was attributed to the decreased growth rates from the experiment with
0.67.
R=
Comparing Figs. 3.3 and 3.4, it appears that only a few data points where located in
the unstable region due to the coarse crack measurements. It is recommended that additional
data be collected in the unstable region to adequately t this model.
3.3
Conclusion
The Paris equation is widely used in fracture mechanics and ts well to the data for each
R-ratio.
The six Paris equations were combined into one with the Walker equation which
accounts for the stress ratio eects and uses the entire data set to t the model. The Forman equation has the added advantage of modeling the asymptotic curvature near fracture.
However, Hudson's experimental data only has a few data points in this region due to the
coarse measurements. The Walker and Forman models had similar RMSE values, suggesting
that the models equally represent data set. As will be discussed in Chapter 4, Abaqus 6-13.3
is not able to incorporate the Forman equation without the use of a user dened subroutine
which is outside the scope of this thesis.
Since no data appeared to be collected in the
threshold region, the Broek and Rice's model could not be used.
For these reasons, the
Walker equation has been selected for the analysis in the remainder of this thesis.
25
Chapter 4
Estimating Fatigue Fracture Life
R-ratio
It has been well documented in the literature that the
[6, 19, 68], mean stress level
[69, 70], and loading frequency [23, 24] can aect fatigue fracture life. To limit the aect of
these fatigue fracture parameters, four loading conditions have been selected from Hudson's
report. Each selection has the same mean stress,
m ,
while increasing the
not possible to select a single loading frequency while maintaining
m .
R-ratio.
It was
Selected loading
conditions are listed in Table 4.1.
Table 4.1:
Loading
Frequency
Condition
(ksi)
(MPa)
(ksi)
(MPa)
(Hz)
15.0
103
15.0
103
0.5/13.7*
15.0
103
10.0
68.9
0.2
0.5
15.0
103
5.00
34.5
0.5
13.7
15.0
103
3.00
20.7
0.67
13.7
* Indicates multiple tests run a dierent frequencies.
The following sections explore ve dierent methods for evaluating the fatigue fracture life of
a Center Cracked Plate (CCP) model: (1) analytical methods based on published equations,
(2) AFGROW, a program common in the aerospace industry used for crack analysis, (3)
Finite Elements with Quarter Points (FEQP), an analysis procedure where the crack is
manually grown which requires remeshing, (4) Virtual Crack Closure Technique (VCCT),
where a crack propagates along a predened path, and (5) the eXtended Finite Element
Method with Phantom Nodes (XFEM-PN) which allows a crack to propagate along a solution
dependent path [5]. First, the specimen model and general assumptions are established.
26
4.1
Center Cracked Plate Model
The tests described in Hudson's 1969 NASA report were performed on a thin CCP specimen
which was assumed to be in a two-dimensional state of plane stress since
where
is the thickness,
is the half-height, and
th
and
t b,
is the half-width. The plate was loaded
by ve pins near the top and bottom edges. For simplication, the pins have been removed
and an evenly distributed stress was assumed across the top and bottom edges of the plate.
According to a study by Isda in 1971 [71], the eects of the boundary condition on the stress
intensity factor are minimal when
h/b 1.5 for a center cracked plate.
As shown in Fig. 4.1,
the simplied model arbitrarily begins 1.5 inches below the centerline of the inner loading
holes. This results in
h/b = 2.0
Fig. 4.1:
4.2
which is within the acceptable range of values.
Simplied Center Crack Plate Model
Analytical Method
The analytical method was selected as a baseline for the comparison to all other techniques.
The concept for the analytical method was to integrate the Walker equation, Eq. (4.1), from
an initial crack length,
fatigue fracture life,
ai , to a critical crack length at fracture, ac , resulting in the estimated
Nf ,
of the specimen. The following sections briey describe the process
and present the results with comparisons to the experimental data.
da
C0
=
(K)m
dN
(1 R)m(1)
27
(4.1)
4.2.1
Model Setup
ac , was calculated by rearranging Eq. (2.1) and using the critical

Kc , resulting in Eq. (4.2). This equation requires an iterative solution
The critical crack length,

stress intensity factor,
since
F = F (a/b).
This equation was solved in MATLAB using a system of non-linear

6
equations and the function fsolve() with the default termination tolerance of 10 .
1
ac =
Kc
F (ac /b) max
2
(4.2)
Eq. (4.3) was substituted into Eq. (4.1) and rearranged for integration. The resulting Eq.
(4.4) does not have a closed form solution in general, because of the dimensionless geometry
factor, thus it was numerically integrated in MATLAB with the
integral()
function. The
code used for this analysis can be found in Appendix F.
K =
Nc
Nf =
Ni
K
(1 R)1
1
dN =
C0
ac
(4.3)
da

ai
K (a)
m
(4.4)
This method assumes self-similar crack growth, thus only the Mode I Stress Intensity Factor
(SIF),
KI , is present.
Applying the Maximum Tangential Stress (MTS) criterion, restated in

= cos1 (1) = 0,
Eq. (4.5), for the crack extension direction conrms this assumption since
with
KII = 0.
= cos1
4.2.2
!
p
2
2
+ KI4 + 8KI2 KII
3KII
2
KI2 + 9KII
(4.5)
Results
Table 4.2 list the results for critical crack length and estimated fatigue fracture life. For a
comparison, the last recorded crack length and fatigue cycles from Hudson's report are also
listed.
The analytical critical crack lengths are approximately twice as large as the nal
measured crack lengths, while the life estimates are within 10-20% of the experiments. One
reason for the discrepancy in crack length is that
Kc
may have been overestimated based on
the visual testing method used in Hudson's report and mentioned in Section 3.1. Also, the
coarse crack measurements of
fracture occurred.
0.1-0.2 in.
are not able capture the exact crack length when
In the unstable region, high growth rates only require a few cycles to
propagate a crack over a signicant length, thus
Nf
is not sensitive to the nal crack length.
This was conrmed by a sensitivity analysis presented in AFGROW's Handbook for Damage
Tolerance Design [72]. Adding to the discrepancy of critical crack length, Hudson's report
28
did not specify whether specimens were fractured in the fatigue crack growth test or were
removed prior to failure for the residual strength test.
Table 4.2:
Loading
Results from the Analytical Fatigue Fracture Analysis
Critical Crack Length,
ac ,
in. (mm)
Fatigue Fracture Life,
Nf ,
cycles
Condition
Experimental
0.8
(20.3)
1.68
(42.7)
3 050
2 280
0.2
1.4
(35.6)
2.24
(56.9)
8 420
7 540
0.5
1.4
(35.6)
2.99
(75.9)
42 500
49 600
0.67
1.8
(45.7)
3.34
(84.8)
154 000
179 200
Analytical
Experimental
Analytical
R = 0 and 0.2, the analytical method gave a conservative life estimate, while it overestimated Nf for R-ratios of 0.5 and 0.67. The inconsistency in fatigue cycle estimates may
For
be attributed to several factors including the loading frequency and plastic zone size at the
crack tip. The loading frequency was not consistent across all experiments, and analysis of
its eects is outside the scope of this thesis.
Fatigue cycles may also be aected by the plastic zone ahead of the crack tip, which may require the use a
J -integral analysis using Elastic Plastic Fracture Mechanics (EPFM). Plastic
zone size was analyzed by Rice [73] in 1967 who developed the formula for the plane stress
plastic zone radius,
rp ,
presented in Eq.
(4.6) for plane stress conditions.
A ow chart
provided by Dowling in Chapter 8 of [8] provides a logical progression to determine if LEFM

is applicable in the presence of the plastic zone. The path followed in the ow chart in [8] is
described as follows. First, for the analysis to be plane strain Eq. (4.7) must be satised for
t is the thickness, h is the half-height, b is the half-width,

K = K (a). Eq. (4.7) was not satised for any of the loading
hand side of the equation was greater than t = 0.09 in. (2.29 mm)
all plate dimensions listed, where
is the crack half-length, and
conditions since the right

for all
a's,
therefore a state of plane stress exists.
rp
1
=
2
K
yld
2

t, (b a) , h 2.5
The next criterion in Eq.
(4.6)
K
yld
2
(4.7)
(4.8) was used to determine if the region of
aected by the geometry of the plate.
In the region of
K -dominance,
K -dominance
at the crack tip can be neglected, and SIFs adequately characterize the stress eld.
K -dominance
K -dominance
be warranted.
The
region must be small relative to the distance to any of the listed geometric
boundaries. A distance of
on the
was
limited plasticity
Eq.
8rp
is typically acceptable [8]. If a geometric boundary intrudes
region, then gross yielding may occur and an EPFM analysis may
(4.8) was satised for all loading conditions.
29
Therefore,
continues
to describe the stress eld surrounding the crack tip and LEFM is still applicable for the
analysis in this thesis.
a, (b a) , h 8rp
(4.8)
Probably the most likely reason for the discrepancy in fatigue fracture life is the fact that
only two tests were run at each loading condition, and only the average results were provided
in Hudson's report.
As Virkler et al. [63] showed from sixty-eight fatigue crack growth
tests, there is large variability in growth rates for the same material, geometry, and loading
condition. As a result it is not likely to draw any denitive conclusions from a two specimen
test.
Fig. 4.2 shows a plot of crack length versus fatigue cycle count for the analytical method and
R = 0.2. To build the analytical curve, Eq. (4.4) was incrementally

ac until a smooth curve was formed. Experimental data was plotted
experimental data with

solved between
ai
and
directly from Table A.2. The experimental data has a Root Mean Square Error (RMSE) of
13.2% calculated with Eq. (4.9), where
calculated with Eq.
is the number of samples. The maximum error
(4.10) was 15.6% compared to the analytical curve.
Eqs.
(4.9) and
(4.10) are used throughout the remainder of this thesis, with the appropriate substitutions
for
N.
Results for all loading conditions can be found in Appendix D.
v
u n
2
u 1 X NCurrent Analysis
t
1 100%
RMSE =
n k=1 NAnalytical Result
k

% Error
= max

NCurrent Analysis
1 100%
NAnalytical Result
30
(4.9)
(4.10)
0.4
Analysis
________
0.35
R = 0.2
= 20.0 ksi
(695 MPa)
Hudson
Analytical
0.3
a/b
0.25
0.2
0.15
0.1
0.05
0
Cycles, N (x10 cycle)
Fig. 4.2:
and experimental data (R
4.3
= 0.2).
AFGROW
AFGROW was used as a comparison tool to a well-established commercial code in the

aerospace, mechanical, and civil industries. Conceptually, all that was needed for this analysis was specimen geometry, a loading spectrum, and material properties. A description of
the model setup and a comparison to the analytical method are provided in the following
sections.
4.3.1
Model Setup
The center cracked plate is a classic, predened model in AFGROW. Required inputs include
the geometric dimensions shown in Fig.
4.1, the elastic material properties and Walker
equation constants listed in Table 4.3, and a constant amplitude loading spectrum. All elastic
material properties are from Hudson's report unless otherwise noted.
Material properties
that did not apply to the crack analysis in this thesis were left at default values. The loading
spectrum consisted of the maximum stress,
max , and R-ratio listed in Table 4.4.
AFGROW
processed crack growth incrementally by increasing the crack length by a length of

until
Kc
was reached.
0.01 in.
Cycles were summed over each increment for the nal estimate of
fatigue fracture life.
31
Table 4.3:
AFGROW Material Data
AFGROW Material Data
(1)
8.63 1010
n(2)
in./cycle
ksi
in.
m/2
3.84
(3)
0.564
10100 ksi
(4)
0.33
Coe. Thermal Expan.
1.25E-05 in./in./ F
yld
75.9 ksi
(default)
Kc
72.012 ksi
KIc (4)
26.0 ksi
Kth
2.0 ksi
in.
in.
in.
(default)
Lower Limit on R shift
-0.3
(default)
Upper limit on R shift
0.8
(default)
C0 in Walker equation Eq. (4.1)

m in Walker equation Eq. (4.1)
(3)
in Walker equation Eq. (4.1)
(1)
(2)
(4)
4.3.2
Found at www.matweb.com [74]
Results
Critical crack lengths and fatigue fracture life estimates for the AFGROW analysis are listed
in Table 4.5.
Compared to the analytical method, critical crack lengths are within 0.6%
and life estimates are within 2.5%.
R = 0.2
A plot of the normalized stress intensity range for
in Fig. 4.3 shows that the stress intensity factors compare well with those from the
analytical method, with an RMSE of 0.4% and a maximum error 0.8% at the nal crack
length. Discrepancies in
KI
values were a result of dierent equations for the dimensionless
F . Additionally, the error in KI 's are within an acceptable range since

F , Eq. (3.3), was only expected to be accurate within 1% of experimental
shows that the cycle count during crack propagation for R = 0.2 compare
geometry factor,
Tada's equation for

data. Fig. 4.4
well, with an RMSE of 2.1% and a maximum error of -2.2% at the nal crack length. Similar
results for
K , ac ,
and
Nf
were found for all loading conditions in Appendix D.
32
Table 4.4:
AFGROW Load Spectrum
max
Loading
Table 4.5:
Loading
Condition
(ksi)
(MPa)
30.0
207
25.0
172
0.2
20.0
138
0.5
18.0
124
0.67
Results from the AFGROW Fatigue Fracture Analysis
Critical Crack Length, in. (mm)
Condition
Experimental
0.8
(20.3)
1.68
(42.9)
1.69
0.2
1.4
(35.6)
2.24
(56.6)
0.5
1.4
(35.6)
2.99
0.67
1.8
(45.7)
3.34
Analytical
AFGROW
Fatigue Fracture Life, cycles

Experimental
Analytical
AFGROW
(42.7)
3 050
2 280
2 220
2.23
(56.9)
8 420
7 540
7 370
(76.2)
3.00
(76.3)
42 500
49 600
48 510
(84.7)
3.33
(85.1)
154 000
179 200
175 300
1.1
Analysis
________
K/(( a)1/2)
Analytical
AFGROW
R = 0.2
1.05
0.95
Fig. 4.3:
0.05
0.1
0.15
method and AFGROW (R
0.2
a/b
KI
= 0.2).
33
0.25
0.3
0.35
0.4
during crack propagation for the analytical
0.4
Analysis
________
0.35
R = 0.2
= 20.0 ksi
(695 MPa)
Hudson
Analytical
AFGROW
0.3
a/b
0.25
0.2
0.15
0.1
0.05
0
Fig. 4.4:
and AFGROW (R
= 0.2).
34
4.4
Finite Elements with Quarter Points (FEQP)
The rst FEM explored uses the method of moving the mid-side nodes of elements adjacent to
the crack tip to their quarter point location. This method was selected because it represents
a classical FEM solution to a fracture mechanics problem, as mentioned in Section 2.2.2.
Only static loads are possible, and the crack must be manually grown and remeshed at each
crack increment. Finite element analyses were performed in Abaqus 6-13.3. Results from
Abaqus were post processed in MATLAB to estimate the cycles to failure,
Nf .
The following
sections discuss the model setup, mesh convergence analysis, and results compared to the
analytical equations..
4.4.1
Model Setup
A two-dimensional half model was developed due to the symmetry condition about the
vertical mid-plane of the plate. Quarter models could not be used due to the limitations of
the contour integral used by Abaqus to calculate
K.
Zero vertical displacement was specied
along the bottom plate edge and zero horizontal displacement along the symmetry mid-plane.
A pressure load was applied along the top surface of the plate to simulate a uniform nominal
axial stress. Material properties included
and
were assigned with a section thickness of
0.09 inches (2.29 mm).

A focused meshing strategy was employed to increase the number of elements in the vicinity
of the crack tip, requiring the partitioning scheme in Fig.
4.5 where the crack has been
highlighted in red. The model used a structured meshing algorithm everywhere, except in
the focused region where a sweeped meshing algorithm was used with element density biased
toward the crack tip.
Two-dimensional plane stress, quadratic quadrilateral elements (CPS8) were selected. Quadratic
elements were required since they have the mid-side nodes needed for the FEQP analysis.
Edges of elements adjacent to the crack tip were collapsed to a single point, and mid-side
nodes were relocated to the element's quarter point location to better represent the stress
singularity at the crack tip. The resulting course mesh is shown in Fig. 4.6.
Circular partitioning used in the meshing strategy serves a dual purpose. Abaqus evaluates
the
J -integral
using contour integrals. A ring of elements around the crack tip are selected
for this calculation. Building the model so the elements fall within a regular circular region
helps avoid uctuations in the results. SIFs are extracted from
for a linear elastic material
using Eq. (4.11) [31].
J =G=

1 2
1
2
KI2 + KII
+
K
2G III
E
(4.11)
An FEQP analysis does not allow for crack propagation, thus the crack was manually grown
from
ai
to
ac .
With this limitation, the Crack Propagation Direction (CPD) was assumed
to be perpendicular to the vertical mid-plane of the plate and uniform nominal axial stress.
35
Fig. 4.5:
Partitioning, loading, and boundary conditions used in the FEQP model with an
inset zoomed in at the crack tip.
This assumes that the crack would propagate by pure Mode I fracture. Each crack increment
required a new FEQP model which was remeshed to the new geometry. A modied Simpson's
rule was used to select the minimum number of crack lengths needed to adequately estimate
life. In Chapter 11 of [8], Dowling suggests using the modied Simpson's rule presented in
Eqs. (4.12) and (4.13) to numerically integrate Eq. (4.4). This set of equations places more
crack increments,
an , close to ai
where the slope of the
y = dN/da is steepest, and gradually

r was found by selecting an
increase the spacing as the slope attens at larger crack lengths.

even number of crack lengths,
(4.12) with
an = ac .
n,
starting at zero for the initial crack length and solving Eq.
The value of
was selected by increasing its value in even increments
until the result of Eq. (4.13) did not change signicantly.
an = r n ai
aj+2
y da =
aj
n2
X
aj (r2 1)
j=0
6r
yj r (2 r) + yj+1 (r + 1)2 + yj+2 (2r 1)
(4.12)
(4.13)
Abaqus calculates the CPD based on the Maximum Tangential Stress (MTS), Maximum
Energy Release Rate (MERR), and
KII = 0
crack extension criteria mentioned in Section
2.1.5. CPD results are used here to conrm the previous crack growth direction assumption.
Abaqus' q-vector, which describe the current crack extension direction, was set parallel to
the crack plane for each model which is essential to extract meaningful SIFs.
36
Fig. 4.6:
Coarse meshing scheme for the FEQP model with an inset zoomed in at the crack
tip.
Considering the time it takes to build and remesh the models for each of the loading conditions and crack lengths, it was decided to parametrize the FEQP model. A script was developed in MATLAB and Python to access Abaqus' Complete Abaqus Environment (CAE)
le to adjust the partitions, remesh the model, and run the analysis. Stress intensity factors
were extracted from the output database (ODB) les, and the values were averaged excluding the rst and second contours as these typically uctuated compared to other contours
which were approximately constant. The uctuation in SIFs for the rst and second contour
was a result of the stress singularity at the crack tip and the selection of the modied shape
functions for the quarter point nodes.
The rst contour was dened by the nodes at the
crack tip, and the second contour included all nodes within elements adjacent to the crack
tip. Subsequent contours include additional rings of elements. After the second contour, the
eects of the stress singularity are reduced.
It is recommended in Abaqus' Analysis User
Guide to discard the rst few contours that exhibit this uctuation [5].
4.4.2
Convergence Analysis
Due to the number of models at dierent crack lengths and loading conditions used in this
analysis, three models were selected to be converged with a single loading condition. Convergence of the three models bounds the accuracy of the solution, assuring that the solutions of
the intermediate crack lengths will also be converged. Additionally, in a Linear Elastic Fracture Mechanics (LEFM) analysis the stress intensity is proportional to the loading. Crack
lengths of
0.1 in. (2.54 mm), 1.0 in. (25.4 mm),

R = 0.
and
3.0 in. (76.2 mm)
were selected to be the
convergence models with
All three models converged immediately with average
KI
values less than 1% error relative
to the analytical results for the rst mesh which is within an acceptable error since Tada's
37
equation for
F , Eq.
(3.3), was only expected to be accurate within 1% of experimental data.
Successive mesh renements showed less than a 0.1% dierence in
KI .
Percent dierences
were calculated with Eq. (4.14), which will be used throughout the remainder of this thesis
with appropriate values substituted for
K .
Table 4.6 lists the percent error in average
KI
values relative to the analytical method and percent dierence between mesh renements.
Mesh renements were based on the number of elements along the radial distance of the
focused region,
re .

% Di
= max
Table 4.6:

KCurrent M esh
1 100%
KP revious M esh
(4.14)
Convergence analysis for the FEQP model

Avg.KI Percent Error Relative to Analytical Results
and Percent Dierence Between Mesh Renement
a,
Crack Length
Analytical
KI ksi
Mesh
in. (mm)
in.
(MPa
mm)
Num. Radial
0.1 (2.54)
1.0 (25.4)
3.0 (76.2)
5.61 (195)
18.0 (624)
36.1 (1250)
% Error
% Di.
% Error
% Di.
% Error
% Di.
Elements
1
-0.22
0.29
0.94
-0.12
0.09
0.31
0.02
0.96
0.02
16
-0.10
0.02
0.31
0.00
0.96
0.00
Fig. 4.7 shows the normalized
KI
variation at dierent contour intervals, where
is the
radius from the crack tip. Contour locations are represented by the markers on each line.
The number of contours in each mesh was increased so that all elements in the circular
focused region were included in each calculation.
Fig.
4.7 shows the uctuation in
KI
values for the rst and second contour as a result of the stress singularity and modied
shape functions. These two contour values were not included when averaging
KI
for the
crack length. The RMSE of all contour intervals and each mesh are listed in Table 4.7. Mesh
#2 had an RMSE of 0.31% relative to the analytical result. Proceeding to the third mesh
does not oer signicant improvement in the RMSE. Results for convergence models with
a = 0.1
and 3.0 in. had similar
Table 4.7:
RMSE of
KI
KI
Contour Values at
Mesh
The average
KI
variations, and the gures are found in Appendix C.
KI
a = 1.0 in.
(25.4 mm) in Fig. 4.7.
RMSE Relative
to Analytical Results (%)
0.33
0.31
0.31
values for rst mesh were in excellent agreement with the analytical
method, but the second mesh was selected for added condence that the mesh density would
38
1.02
a/b = 0.17
K/(( a)1/2)
1.019
1.018
1.017
Mesh
____
#1
#2
#3
Analytical
1.016
1.015
1.014
1.013
Fig. 4.7:
0.1
0.2
FEQP normalized
0.3
KI
0.4
0.5
r/re
0.6
0.7
0.8
be adequate for the intermediate crack lengths. Additionally, more
0.9
a = 1.0 in.
KI
(25.4 mm)
contour values were
averaged with this selection. All models solved in under a minute so the added processing
time was insignicant.
4.4.3
Results
Based on the modied Simpson's rule analysis, it was found that ten crack lengths are
adequate to estimate life. However, twenty crack lengths were used for added delity in the
crack growth curves.
Critical crack lengths for this analysis were based on the analytical
method since these were needed to setup the models. Therefore, critical crack lengths for
the FEQP and analytical analyses are the same, as seen in Table 4.8. Fatigue fracture life
estimates compare well with the analytical method with errors less than 0.5%.
Table 4.8:
Loading
Results from the FEQP Fatigue Fracture Analysis
ac ,
in. (mm)
Condition
Experimental
0.8
(20.3)
1.68
(42.7)
1.68
0.2
1.4
(35.6)
2.24
(56.9)
0.5
1.4
(35.6)
2.99
0.67
1.8
(45.7)
3.34
Analytical
FEQP
Nf ,
cycles
Experimental
Analytical
FEQP
(42.7)
3 050
2 280
2 270
2.24
(56.9)
8 420
7 540
7 530
(75.9)
2.99
(75.9)
42 500
49 600
49 490
(84.8)
3.34
(84.8)
154 000
179 200
178 780
Plots of the normalized stress intensity variation and cycle count during crack growth for
R = 0.2 are shown in Figs.
4.8 and 4.9, respectively. SIFs compare well as expected from the
convergence analysis with an RMSE of 0.3% and maximum error of 0.8% at the nal crack
length, both within an acceptable range compared to Tada's equation for
39
F.
As a result of
the minimal error in the stress intensity factors, the crack length versus cycle count was also
in excellent agreement with the analytical method with an RMSE of 0.1% and a maximum
error of -0.2% at the nal crack length. Results in Appendix B for the remaining loading
conditions show similar behavior for both SIFs and fatigue cycles versus crack length.
1.1
K/(( a)1/2)
Analysis
________
Analytical
AFGROW
FEQP
1.05
R = 0.2
0.95
Fig. 4.8:
0.05
0.1
0.15
and FEQP methods (R
0.2
a/b
KI
0.25
0.3
0.35
0.4
= 0.2).
A CPD analysis shows that the MTS, MERR, and
KII = 0
methods resulted in crack
extension angles of zero degrees for all a/b < 0.3. Crack extension angles for a/b > 0.3 do
not exceed 0.05 , which is a result of numerical noise and was ignored. This analysis conrms
the prior assumption that pure Mode I crack growth occurs perpendicular to the mid-plane
of the plate. CPD angles were plotted in Fig. 4.10 for
are found in Appendix B.
40
R = 0.2
and other loading conditions
0.4
Analysis
________
0.35
0.3
0.25
a/b
R = 0.2
= 20.0 ksi
(695 MPa)
Hudson
Analytical
AFGROW
FEQP
0.2
0.15
0.1
0.05
0
Fig. 4.9:
Comparison of cycle count during crack propagation for the analytical and FEQP
methods (R
= 0.2).
Crack Propagation Direction, (degrees)
0.15
0.1
R = 0.2
0.05
0
0.05
CPD
Method
__________
MTS
MERR
KII=0
0.1
0.15
Fig. 4.10:
0.05
0.1
0.15
0.2
a/b
0.25
0.3
0.35
0.4
Comparison of crack propagation direction for the MTS, MERR, and
criterion in the FEQP analysis (R
= 0.2).
41
KII = 0
4.5
Virtual Crack Closure Technique (VCCT)
A goal of this thesis was to evaluate XFEM-PN's ability to model crack propagation. The
VCCT analysis serves as a logical step to this goal since XFEM-PN uses a
method which will be discussed in Section 4.6.
modied
VCCT
Abaqus' VCCT low-cycle fatigue analysis
allows for automated crack propagation along a predened path using a single model with no
remeshing. Principles of LEFM were used to propagate the crack which provided a parallel
comparison to the previous methods.
Analyses were performed on an IBM iDataPlex system due to the increased memory requirements and solution time. Abaqus 6-13.1 was available on the system; however, the results of
the rst mesh were compared to those from Abaqus 6-13.3 and no dierences were observed.
The following sections discuss the model setup, mesh convergence, and results.
4.5.1
Model Setup
A two-dimensional half model was developed in this analysis using two parts; similar top
and bottom halves of the plate.
Each part had the same material properties and section
as dened in the FEQP analysis. The halves of the plate were placed in surface-to-surface
contact; specifying the propagation plane, labeled in Fig.
4.11, as the initially bonded
region. The unbonded portion represents the initial crack. The assembled model had the
same boundary conditions as the FEQP model.
The goal of the focused meshing strategy in the VCCT model was to maximize the number of
elements along the propagation plane while reducing the density elsewhere, as it was expected
to require many elements to converge
Nf .
The partitioning scheme shown in Fig. 4.11 allows
the number of elements to gradually and uniformly increase while using a structured meshing
algorithm. Several of these square patterns were stacked together to build a 12:1 element
transition, meaning for every one partition at the top and bottom plate edges, there are
twelve partitions at the crack propagation plane.
Two-dimensional plane stress, linear quadrilateral elements (CPS4) were selected since these
elements are the only quadrilateral elements supported by VCCT in a low-cycle fatigue
analysis [5]. Tetrahedron elements are supported; however, these elements tend to be too
sti. The resulting course mesh is shown in Fig. 4.12.
Two steps were required in this analysis: a static general, and direct cyclic step. The static
general step serves two purposes: rst, it initializes the crack prior to the direct cyclic step,
and second, it loads the structure up to
max
and
min .
max
so that the cyclic load can oscillate between
An example of the loading amplitude curve is shown in Fig. 4.13 for
R = 0.2.
Dening a sinusoidal curve tended to reduce the stabilization ratios, discussed next, by an
order or two of magnitude compared to a triangular shaped curve. Assisting stabilization of
each cycle reduced the computation times.
42
Fig. 4.11:
Partitioning, loading, and boundary conditions used in the VCCT model with
Application of a cyclic load requires Abaqus' direct cyclic solver, which uses a truncated
Fourier series to solve for the displacement at each node at all time increments. Dierences
between the applied load and the internal forces, known as the residuals, are also calculated at
each node and time increment. Abaqus' representations for nodal displacement and residuals
based on a truncated Fourier Series are presented in Eqs.
(4.15) and (4.16), respectively.
Displacement coecients are corrected after each iteration until a stabilized solution is found.
The number of Fourier terms,
n,
can increase between iterations. Stabilized solutions are
dened by four ratios which must be below a set tolerance. Stabilization ratios are dened
as: (1) the ratios of the largest correction to the displacement coecients over the largest
CU0 , and periodic, CUn , terms and (2) the ratios of

the largest residual coecient to the time average force for the constant, CR0 , and periodic,
CRn , terms. [5, 31]
displacement coecient for the constant,
u(t) = u0 +
n
X
(usk sin kt + uck cos kt)
(4.15)
(Rks sin kt + Rkc cos kt)
(4.16)
k=1
R(t) = R0 +
n
X
k=1
43
Fig. 4.12:
Coarse meshing scheme for the VCCT model with an inset zoomed in at the
initial crack.
where:
u is
u0
the nodal displacement
is an unknown constant displacement coecient
usk , uck
R is
R0
the nodal residual
is an unknown constant residual coecient
Rks ,
are unknown periodic displacement coecients
Rkc
are unknown periodic residual coecients
is angular frequency
t is
time.
Several steps were taken to aid in the stabilization of the VCCT low-cycle fatigue analysis.
For a propagating crack analysis, Abaqus recommends increasing the parameter
I0 ,
which
is the frequency at which residuals of successive iterations are checked to determine if they
are increasing.
analysis option.
I0
was increased from a default value of 4 to 8 using the discontinuous
The discontinuous analysis option also increases
IR ,
the starting point
for the logarithmic rate of convergence checks, from a default of 8 to 10.
However, this
convergence check was not performed in the direct cyclic step since a xed time increment
was used.
IA
was increased from the default value of 5 to a recommended value of 10 in
the static step and 20 in the direct cyclic step. The
IA
parameter determines the maximum
number of cutbacks allowed in an increment. [5, 31, 75]

Processing time and solution accuracy were used as metrics in a tradeo study to determine
the number of Fourier terms, values for stabilization ratios, and time increment size needed
for a stabilized cycle. The number of Fourier terms was limited to 12 since the default value
of 25 resulted in an increase in processing time while not converging to the stabilization
44
max
Load Amplitude
R = 0.2
0.75max
0.50max
Step
____
0.25max
Static
Cyclic
0
Fig. 4.13:
0.5
1
Step Time, t
1.5
Example load amplitude curve for the static and cyclic steps with
tolerances. Ratios
R = 0.2.
CU0 and CUn were left at default values of 0.005 and CR0 and CRn where
both increased to values between 0.15 to 0.30 depending on mesh density. Additional results
for
CR0
can be found in Appendix C. Since no temperature eects were analyzed, Abaqus
used a xed time increment set to 0.02.
Larger time increments resulted in destabilized
cycles.
Fracture and fatigue criterion were needed for the static and direct cyclic steps, respectively.
The Benzeggagh-Kenane, or BK, model [15] was selected to describe the equivalent critical
energy release rate,
Gc ,
for a mixed mode fracture analysis. Eq. (4.17) presents the BK law
as interpreted by Abaqus [5]. This form of
Gc
was selected because failure was assumed to
occur from pure Mode I fracture as seen in the FEQP analysis and the following order of
magnitude study showed that
Gc GIc .
In a two dimensional problem
assumed that the CCP in this thesis failed by pure Mode I fracture, so
was no material data to empirically t the exponent,
GIII = 0.
GII GI .
Gc GIc , it was assumed that GIIc G Ic .

results in Gc GIc as needed for this analysis.
(4.17)

There
so a value of unity was assumed.
No data was available to determine the pure Mode II critical energy release rate,
ensure that
It was
GIIc .
To
Applying these assumptions to Eq.
GII + GIII
GI + GII + GIII

(4.17)
The fatigue crack propagation portion of the analysis also needed a denition for crack
growth. Abaqus has dened the onset of fracture with Eq. (4.18) and growth dened by the
Paris equation, modied for the energy release rate range,
c2 , c3 ,
and
c4
G ,
in Eq. (4.19). Constants
c1 ,
are empirically determined based on material and geometry; and the energy
release rate range is dened as
G = Gmax Gmin
[5]. Constants
c1
and
c2
were both set to
zero in this study since no data was available to t them and this would assume growth to
45
start immediately. To solve for constants c3 and c4 , Eq. (4.19) was set equal to Eq. (4.1).
2
Using G = Gmax (1 R ) and the relation between G and K , the relationship between
Abaqus' constants and Walker equation constants was formed in Eq. (4.20). Values of
and
c4
c3
for each loading condition are listed in Table 4.9.
N = c1 (G)c2
(4.18)
da
= c3 (G)c4
dN
(4.19)
c3 = C0 (1 R)
c4 =
Table 4.9:
Loading
E
1 R2
m/2
(4.20a)
m
2
(4.20b)
Abaqus Crack Growth Constants.
c3
c4
in./cycle
(ksi in.)c4
mm/cycle
(MPa mm)c4
4.19 102
5.26 105
1.92
0.2
2.79 10
3.51 10
1.92
0.5
1.62 102
2.04 105
1.92
1.92
Condition
1
0.67
1.20 10
Abaqus requires bounds on the range of
1.51 10
applied to the Paris equation since the equation
is only valid in the Paris region of crack growth.
The energy release rate threshold,
Gth ,
denes the lower bound where no crack growth occurs and the Paris limit energy release
rate,
Gpl ,
denes the upper bound where the crack growth rate enters the unstable growth
region [5].
Dowling [8] recommends using the critical stress intensity as the upper bound
to this region; thus,
Gpl GIc .
Gpl /Gc
was set to 0.999, allowing growth to occur up until fracture at
Since data provided in Hudson's report did not exhibit signicantly slower growth
rates approaching a threshold value, the ratio of
Gth /Gc
was set to 0.001, generally allowing
crack growth to occur at all applied loads.

During Abaqus' VCCT fatigue crack growth analysis, values of
quantities at the crack tip with Eq.
Inserting
are calculated as nodal
(4.21), restated here and dened in Section 2.2.2.
into Eq. (4.19) and using the characteristic element length of the model, the
number of cycles required to fracture the element ahead of the crack tip are calculated. In a
three-dimensional model, the element, or elements, with the least number of cycles remaining
within a given tolerance are fractured [5]. For a VCCT analysis, this means the contacting
nodes on the propagation plane are released. In a two-dimensional model, there is only one
46
element ahead of the crack tip, thus only one element can be fractured during each cycle.
Once the nodes have been released, the calculated number of cycles is added to the total
cycle count and the analysis repeats, redistributing the load over the reduced cross section.
This process continues until the fracture ratio,
Gpl /Gc , is reached, at which point propagation
occurs in one cycle increments until reaching the user dened number of cycles to analyze [5].
GI =
1
Fy,cd uy,ab
2a
(4.21a)
GII =
1
Fx,cd ux,ab
2a
(4.21b)
GIII =
1
Fz,cd uz,ab
2a
(4.21c)
One of the objectives of this thesis is to compare the estimated life of dierent analysis
procedures. The nal cycle count, as calculated by Abaqus, was unknown, and these analyses
were computationally expensive with some of the nal mesh renements taking several days
to complete. Having an analysis run after it had reached
Gpl /Gc
was undesirable. Therefore,
a Python code was written to run in parallel with Abaqus to monitor the analysis's status
(STA) le, looking for consecutive cycles to increase by one. If this occurred, it was assumed
that
GIc
had been reached, fracture had occurred, and the analysis could be terminated.
This code is provided in Appendix G.

Fatigue crack propagation was achieved by coupling VCCT, the direct cyclic solver, and the
fatigue and fracture criterion; however, this operation is not directly supported by Abaqus
CAE. All denitions made in the direct cyclic step were manually added to the input (INP)
le.
4.5.2
Two convergence criteria were used for the VCCT analysis: the stress intensity range,
KI ,
Nf . Abaqus reports the energy release

rate range, G , and Eq. (4.22) was used to convert to K based on Abaqus' denition of
G = Gmax Gmin . Loading condition #1 with R = 0 was used, similar to the FEQP
as the crack propagated, and the nal cycles at failure,
convergence analysis.
r
K =
1R
EG
1+R
(4.22)
KI , calculated with Eqs. (4.9) and

ai through ac where Kmax = Kc . Abaqus does not
= Kc so results were truncated to meet this condition.
Table 4.10 lists the RMSE and maximum error of

(4.10), as the crack propagates from
terminate the analysis when
Kmax
Error calculations are based on the analytical method. As shown in Fig. 4.14, the maximum
47
errors occur with small
a/b 0.03,
a/b
ratios.
By the third mesh, nodal errors are less than 2% for
and the RMSE including all nodes is 1.0%.
acceptable range compared to Tada's Eq.

within 1% of experimental data.
Fig.
The
KI
RMSE is within an
(3.3) which was only expected to be accurate
4.14 indicates the importance of a ne mesh with
short crack length ratios.

Percent dierences in Table 4.10 were calculated between consecutive meshes at nodal locations consistent with Mesh #1. The Root Mean Square Dierence (RMSD) calculation
is presented in Eq. (4.23), where
is the total number of nodes in Mesh #1. Maximum
dierences were calculated with Eq. (4.14). RMSD were 0.5% between Mesh #3 and #4.
RMSE and RMSD for
KI
in Mesh #3 indicate that this mesh would be sucient for the
rst metric of this analysis.
v
u n
2
u 1 X KCurrent M esh
1 100%
RMSD = t
n k=1 KP revious M esh
k
Table 4.10:
(4.23)
Stress Intensity Range Convergence Analysis.
Elements Along
KI
% Error
Propagation Plane
Max.
RMSE
Max.
RMSD
120
-17
4.5
240
-7.3
2.0
11
2.7
480
-3.5
1.0
4.1
1.0
960
-1.8
0.9
1.8
0.5
Mesh
KI
% Dierence
K Percent Error (%)
10
15
20
Fig. 4.14:
Mesh
____
#1
#2
#3
#4
0
0.05
0.1
Error in stress intensity range,
0.15
a/b
0.2
0.25
0.3
KI , during crack propagation in VCCT analysis.

48
Fatigue fracture life estimates were the second metric for convergence, and results of the
analysis are listed in Table 4.11.
Percent error was based on the analytical method, and
percent dierence was between consecutive meshes. Two methods of computing life estimates
were used in this analysis. The rst estimate is the fatigue cycle as reported by Abaqus when
Kmax = Kc .
The second method, the integrated polynomial method, integrates the
da/dN
curve. This was achieved by tting a quadratic polynomial to three consecutive points with
the nodal location on the abscissa and
da/dN
calculated from Eq. (4.19) on the ordinate. The
incremental curves were integrated and the results summed over the crack length, resulting
in
Nf .
This method is similar to the previously described modied Simpson's rule used
in the FEQP analysis; however, in general this method does not restrict the location of
middle point since it's location is mesh dependent. This was performed in MATLAB with
the
polyfit(), polyint(),
and
polyval()
functions. Code can be found in Appendix F.
It is obvious from Table 4.11 that the estimated fatigue fracture life is highly mesh dependent,
with errors in the Abaqus estimate relative to the analytical results greater than 97% for the
rst mesh and reducing to 9% by the fourth mesh. Using the integrated polynomial approach
shows a signicant reduction in error.
Fig.
4.15 shows that the integrated polynomial
method gives life estimates approximately one mesh renement before Abaqus' calculation,
thus requiring fewer elements and less computational time for a similar answer. Mesh #3's
estimate of
Nf
is within 2.8% of the analytical results. Further rening the mesh shows a
3% change in the answer, while there is a signicant increase in processing time.
Table 4.11:
Mesh
Estimated Fatigue Fracture Life Convergence Analysis.
Elements Along
Abaqus,
Nf
Integrated Polynomial,
Nf
Processing
Propagation Plane
% Error
% Di.
% Error
% Di.
Time (hrs)
120
97
38
0.5
240
34
-32
11
-19
3.0
480
15
-14
2.8
-7.4
21
960
9.2
-5.3
-0.2
-3.0
184
Mesh #3 was considered converged since RMSE was within 1% for
KI
and 3% for
Nf
using the integrated polynomial method. RMSD's for both metrics do not show a signicant
improvement by proceeding to Mesh #4. Processing time was also taken into consideration
with Mesh #4 taking 184 hours to run, which was not feasible to run multiple models.
4.5.3
Results
Results from the VCCT analysis are shown in Table 4.12. Critical crack lengths for the VCCT
analysis compare well with the analytical method with errors less than 0.5% attributed to the
discretization of the mesh. Fatigue fracture life estimates, using the integrated polynomial
method, are approximately 8% higher than the analytical results, which is a result of a slower
growth rate from under estimating
KI
at the shorter crack lengths where the most cycles
49
Life Estimate Percent Error (%)
100
Method
______
Abaqus
Int. Poly.
Analytical
80
60
40
20
0
100
200
Fig. 4.15:
300
400
500
600
700
800
Number of Elements Along Propagation Plane
900
1000
VCCT Convergence plot of estimated fatigue fracture life.
are accumulated. Normalized
KI
values are plotted in Fig. 4.16, where the error in shorter
crack lengths is obvious. This is a result of the converged mesh selection, where a similar
curved shape is seen for Mesh #3 in Fig. 4.14. The RMSE for the
R = 0.2 load was 0.9%

KI RMSE is within
with a maximum error of -3.6% found at the initial crack length. The
an acceptable range compared to Tada's dimensionless geometry factor error.
Table 4.12:
Loading
Results from the FEQP Fatigue Fracture Analysis
Critical Crack Length, in. (mm)
Condition
Experimental
0.8
(20.3)
1.68
(42.7)
1.69
0.2
1.4
(35.6)
2.24
(56.9)
0.5
1.4
(35.6)
2.99
0.67
1.8
(45.7)
3.34
Analytical
Crack growth versus cycle count for
Fatigue Fracture Life, cycles
VCCT
Experimental
Analytical
(42.9)
3 050
2 280
2 460
2.25
(57.2)
8 420
7 540
8 150
(75.9)
3.00
(76.3)
42 500
49 600
53 630
(84.8)
3.35
(85.1)
154 000
179 200
193 900
R = 0.2
is shown in Fig. 4.17 with an RMSE of 8.6%
and a maximum error of 13% found during the rst growth cycle. Similar
cycle results were found for all
R-ratios
VCCT
KI
and fatigue
with gures in Appendix D. No CPD analysis was
performed since the crack is restricted to growth along the plane perpendicular to plate's
mid-plane. As a check on the solution, all values of
with
=0
was expected.
50
KII KI ,
thus pure Mode I fracture
1.1
K/(( a)1/2)
Analysis
________
Analytical
AFGROW
FEQP
VCCT
1.05
R = 0.2
0.95
Fig. 4.16:
0.05
0.1
0.15
and VCCT methods (R
0.2
a/b
KI
0.25
0.3
0.35
0.4
= 0.2).
0.4
Analysis
________
0.35
0.3
0.25
a/b
R = 0.2
= 20.0 ksi
(695 MPa)
Hudson
Analytical
AFGROW
FEQP
VCCT
0.2
0.15
0.1
0.05
0
Fig. 4.17:
Comparison of cycle count during crack propagation for the analytical and VCCT
methods (R
= 0.2).
51
4.6
eXtended Finite Element Method with Phantom Nodes

(XFEM-PN)
XFEM-PN was the nal FEM explored which allows a crack to propagate along a solution
dependent path without remeshing. This method is advantageous when analyzing complex
loading and geometry where no analytical solution or experimental data are available. In
the present analysis, analytical methods and experimental data are available so a direct
comparison was made.
XFEM with crack tip enrichment functions are only available in Abaqus 6-13.3 with a nonpropagating crack and only in three dimensions. Thus, a comparison to this method is not
possible within the scope of this thesis.
Analyses were performed with Abaqus 6-13.3 on a 64-bit personal computer with 20.0 Gb of
RAM. Abaqus 6-13.1 on the IBM iDataPlex system was initially used; however, a problem
with the Abaqus 6.-13.1 code was found that generated incorrect values for
G .
The model
setup for XFEM-PN low-cycle fatigue analysis used a similar process to the VCCT analysis
in Section 4.5. The model setup, convergence analysis, and a comparison of results to the
previous analyses are presented below.
4.6.1
Methods
Two parts were needed to dene a model with an initial crack: (1) a two-dimensional half
model of the plate, and (2) a two-dimensional wire representing the initial crack. The plate
had the same section and boundary conditions as dened in the FEQP and VCCT analyses.
Traction separation and damage evolution parameters, discussed later, were added to the
material denition along with the elastic properties used in the previous models.
A focused meshing strategy similar to that used in the VCCT model was developed here with
the same 12:1 element transition. The model was assembled and meshed such that the initial
crack was centered on several elements within the enrichment region. Placing initial cracks
on element boundaries was avoided as it tended to give varying and incorrect results which
was attributed to a misrepresentation of the initial crack.
Two-dimensional plane stress,
linear quadrilateral elements (CPS4) were selected since it is the only quadrilateral element
supported in an XFEM-PN low-cycle fatigue analysis [5]. Linear and quadratic triangular
elements were avoided since they produced uctuation inG values and the triangular elements made the structure too sti. Partitioning and meshing results are show in Figs. 4.18
and 4.19, respectively.
The main dierence between VCCT and XFEM-PN analyses was the representation of the
crack. XFEM-PN cracks are placed in an enrichment region where elements have coincident
phantom nodes and phantom degrees of freedom, similar to the process development by Song,
Areias, and Belytschko [4] discussed in Section 2.2.2. Once fatigue or fracture requirements
were met in an element, Abaqus fractures the element, releasing the real node from its
coincident phantom node. The two nodes are then free move independently of each other.
52
Fig. 4.18:
Partitioning, loading, and boundary conditions used in the XFEM model with
XFEM-PN assumes that the stress singularity at the crack tip is not needed in the analysis
and denes failure on
Gc
criterion [5].
A Level Set Method (LSM), similar to the one proposed by Osher and Sethian [56], is
used to locate the discontinuity in the displacement eld of fractured elements.
representation uses two nearly orthogonal signed distance functions,
crack. The intersection of the distance to the crack surface,
crack tip,
and
Abaqus'
to locate the
, and the distance to the initial
denes the crack surface [5]. Only one crack tip is allowed in an element, and
the crack cannot bifurcate or coalesce. Once propagation begins, Abaqus assumes that the
crack will completely fracture the element; therefore, only
and
values need to be stored. Both
are returned as nodal values in the ODB le for the initial cycle associated with
cracking the element, and only
for all subsequent cycles.
Two steps were used for the XFEM-PN low-cycle fatigue analysis in this study: (1) the static
general step loaded the plate up to
max ,
and (2) the cyclic load was applied in the direct
cyclic step to propagate the crack. Similar loading proles to Fig. 4.13 where used. The
direct cyclic solver uses the truncated Fourier series discussed in Section 4.5.1 using Eqs.
(4.15) and (4.16) with control parameters
I0 , IR ,
and
IA
adjusted as in the VCCT analysis.
Processing time and solution accuracy were used as metrics in a similar tradeo study to
determine the number of Fourier terms, values for stabilization ratios, and time increment
size needed for a stabilized cycle. Abaqus' XFEM-PN low-cycle fatigue analysis allows for
parallel computations, thus dramatically reducing processing time compared to VCCT which
solves on a single processor.
Each iteration solved quickly with an increment size of 0.05
and within the default ratio tolerance of 0.005 for
CR0 , CRn , CU0 ,
and
CUn .
Fourier terms
were limited to 12; however, never more than 2 were needed in all XFEM-PN analyses.
53
Fig. 4.19:
Coarse meshing scheme for the XFEM model with an inset zoomed in at the
initial crack.
Similar fracture and fatigue criterion were used for the constants of Eqs. (4.17), (4.18), and
(4.19), and bounding ratios
Gth /Gc
and
Gpl /Gc .
Traction separation and damage evolution
are required by Abaqus for an XFEM-PN analysis; however, these parameters are primarily
used for static load applications, so they did not aect the fatigue crack growth analysis.
For completeness, the maximum principal stress criterion was selected using the ultimate
stress,
ult = 83.2 ksi
(574 MPa). Damage evolution was selected to be a mode-independent
fracture energy equal to
GIc .
Energy release rate ranges were calculated based on a
modied
VCCT method similar to
the method described above; however, nodal forces and displacements are not available at
the crack tip since the crack tip is internal to the element. Abaqus' methods for calculating
the values needed for Eq. (4.21) are proprietary information, not available to the public. A
similar method was mentioned in the paper by Shi et al. [58] where the authors developed
a user dened subroutine for Abaqus, but no insight was provided on the
Once
modied
VCCT.
was known, cycles were calculated from Eq. (4.19), and the crack was propagated
entirely through the element ahead of the crack tip. This was continued until
Gpl /Gc
was
reached. Python code was run in parallel to terminate the analysis should cycles increment
by one, the limit case for critical crack growth.
Abaqus calculates CPD based one on three criteria: MTS, normal to the element 1-direction,
or normal to the element 2-direction.
The latter two methods grow the crack in a mesh
dependent fashion, which removes the need for crack extension models. For comparison to
the FEQP CPD analysis, only the MTS criterion was selected.
All parameters were entered manually into the INP le since Abaqus CAE does not directly
support an XFEM-PN low-cycle fatigue analysis.
54
4.6.2
The stress intensity ranges and cycles to failure were used as convergence criteria with
R = 0,
similar to the VCCT analysis. Table 4.13 shows that the RMSE, computed relative to the
analytical method with Eq. (4.9), holds relatively constant at approximately 4.5% for all
meshes.
Fig.
4.20 shows that the percent error in
KI
at the nodes remains relatively
constant throughout crack growth in Meshes #2 through #4. The uctuations in error seen
in Mesh #1 are attributed to the coarse meshing at the borders of the enriched region.
RMSD in
KI ,
calculated with Eq. (4.23), is 0.2% between Meshes #2 and #3. While the
RMSE is larger than the error reported in Tada's Eq. (3.3) for
F,
proceeding to Mesh #3
does not provided signicantly better results. For these reasons, Mesh #2 was considered
converged for the rst convergence criteria.
Table 4.13:
XFEM-PN Stress Intensity Range Convergence Analysis.

Elements Along
KI
% Error
Propagation Plane
Max.
RMSE
Max.
RMSD
120
-6.4
4.4
240
-4.6
4.5
2.0
0.4
480
-4.8
4.7
-0.3
0.2
960
-5.0
4.9
-0.3
0.3
Mesh
KI
% Dierence
0
Mesh
____
#1
#2
#3
#4
K Percent Error (%)
1
2
3
4
5
6
0
Fig. 4.20:
0.05
0.1
0.15
a/b
Error in stress intensity range,
KI ,
0.2
0.25
0.3
during crack propagation in XFEM-PN
analysis.
The second convergence criterion was the estimated
Nf
calculated by Abaqus and the in-
tegrated polynomial method. Results are listed in Table 4.14. The integrated polynomial
approach has dierences less than 1% between the second and third meshes. Errors for the
55
integrated polynomial method are similar for each mesh, with Mesh #2 21% higher than the
analytical result. Fig. 4.21 shows that Abaqus' results are converging towards the results
from integrated polynomial method.
Table 4.14:
Mesh
XFEM-PN Estimated Fatigue Fracture Life Convergence Analysis.
Elements Along
Abaqus,
Nf
Integrated Polynomial,
Nf
Processing
% Error
% Di.
% Error
% Di.
Time (hrs)
120
60
23
0.25
240
37
-15
21
-2.0
1.0
480
29
-5.8
22
0.5
5.0
960
27
-1.6
22
0.7
37
Life Estimate Percent Error (%)
Propagation Plane
60
Method
______
Abaqus
Int. Poly.
Analytical
50
40
30
20
10
0
100
200
Fig. 4.21:
300
400
500
600
700
800
Number of Elements Along Propagation Plane
900
1000
XFEM-PN Convergence plot of estimated fatigue fracture life.
Mesh #2 was considered converged by both convergence criteria since progressing to a more
rened mesh showed less than a 1% change in results.
Processing time was signicantly
lower compared to the VCCT analysis due to parallel computing.
4.6.3
Results
Results from the XFEM-PN analysis of the four loading conditions are listed in Table 4.15.
Critical crack lengths are 4-7% longer than those calculated with the analytical method.
The longer crack length is a result of low SIFs values.
KI
for
R = 0.2
Fig.
4.22 shows the normalized
which has an RMSE of 4.4% and a maximum error of -4.6% at the initial
crack length relative to the analytical equations. Low SIFs increased the estimated cycles to
failure with errors of approximately 20% for each
56
R-ratio
using the integrated polynomial
method. Fig. 4.23 shows that XFEM-PN has the largest estimate of fatigue fracture life of
the analyses performed in this thesis with a RMSE of 20% and a maximum error of 22% for
the rst growth cycle.
Table 4.15:
Loading
Results from the XFEM-PN Fatigue Fracture Analysis
ac ,
in. (mm)
Experimental
0.8
(20.3)
1.68
(42.7)
1.80
0.2
1.4
(35.6)
2.24
(56.9)
2.38
0.5
1.4
(35.6)
2.99
(75.9)
0.67
1.8
(45.7)
3.34
(84.8)
Condition
Analytical
XFEM-PN
Nf ,
cycles
Experimental
Analytical
XFEM-PN
(45.7)
3 050
2 280
2 750
(60.3)
8 420
7 540
9 090
3.13
(79.4)
42 500
49 600
59 680
3.48
(88.3)
154 000
179 200
215 600
1.1
K/(( a)1/2)
Analysis
________
Analytical
AFGROW
FEQP
VCCT
XFEMPN
1.05
R = 0.2
0.95
Fig. 4.22:
0.05
0.1
0.15
and XFEM-PN methods (R
0.2
a/b
KI
0.25
0.3
0.35
0.4
= 0.2).
A CPD analysis showed that the MTS method predicted
=0
for all crack lengths, as seen
in Fig. 4.24 and as expected with pure Mode I fracture. Abaqus does not return a value for
CPD, so it was post-processed with Eq. (4.24), where
level set values,
is the current fractured element, and
x and y coordinates are based on

i + 1 is the fractured element in the
next cycle. Similar CPD results were found for the other loading conditions are provided in
Appendix D.
= tan1
yi yi+1
xi xi+1
57
(4.24)
0.4
Analysis
________
0.35
0.3
0.25
a/b
R = 0.2
= 20.0 ksi
(695 MPa)
Hudson
Analytical
AFGROW
FEQP
VCCT
XFEMPN
0.2
0.15
0.1
0.05
0
10
Fig. 4.23:
Comparison of cycle count during crack propagation for the analytical and
XFEM-PN methods (R
= 0.2).
0.15
0.1
0.05
0
CPD
Method
__________
0.05
MTS
MERR
K =0
0.1
II
XFEMPN
0.15
Fig. 4.24:
(R
R = 0.2
0.05
0.1
0.15
0.2
0.25
a/b
0.3
0.35
0.4
0.45
0.5
Comparison of crack propagation direction from FEQP and XFEM-PN analyses
= 0.2).
58
4.7
Conclusion
Five methods were used to estimate the fatigue crack growth in a center crack plate: (1)
analytical equations, (2) AFGROW, (3) FEQP, (4) VCCT, and (5) XFEM-PN. Stress intensity ranges from
ai
to
ac
for all methods compared well to each other with the largest
discrepancy seen at short crack lengths for the VCCT, and a constant oset of -5% throughout crack growth for XFEM-PN. Critical crack lengths for all of the analyses are summarized
in Table 4.16. XFEM-PN had the longest critical crack lengths that were 4-7% larger than
the analytical results, which was attributed to under estimating the SIFs.
Table 4.16:
Comparison of Critical Crack Lengths
Loading
Critical Crack length,
ac ,
in. (mm)
Condition
Experimental
0.8
(20.3)
1.68
(42.7)
1.69
(42.7)
1.68
(42.7)
1.69
(42.9)
1.80
(45.7)
0.2
1.4
(35.6)
2.24
(56.9)
2.23
(56.9)
2.24
(56.9)
2.25
(57.2)
2.38
(60.3)
0.5
1.4
(35.6)
2.99
(75.9)
3.00
(76.3)
2.99
(75.9)
3.00
(76.3)
3.13
(79.4)
0.67
1.8
(45.7)
3.34
(84.8)
3.33
(85.1)
3.34
(84.8)
3.35
(85.1)
3.48
(88.3)
Analytical
AFGROW
FEQP
VCCT
XFEM-PN
A crack propagation direction analysis showed good agreement between MTS, MERR, and
KII = 0 in the FEQP analysis and MTS in the XFEM-PN analysis. Each analysis indicated
pure Mode I growth with = 0 as seen in Fig. 4.24. 's that did deviate from zero where
considered numerical noise and were ignored.
One of the main goals of this thesis is to estimate fatigue fracture life. Table 4.17 summarizes the results from each of the analyses. Resulting cycles to failure from the analytical,
AFGROW, and FEQP analyses compared well to each other.
AFGROW and the FEQP
analyses had errors of 0.5% and 2.5%, respectively, relative to the analytical results.
An
integrated polynomial method for calculating life based on the results from Abaqus was
introduced to aid convergence of the VCCT and XFEM-PN models. Compared to the analytical method, VCCT and XFEM-PN overestimated life by 8% and 20%, respectively, which
is a result of the under estimating the SIFs.
None of the methods compare well with the experimental results from Hudson's report.
Several factors are not accounted for with the Walker equation; however, discrepancies in
cycles to failure are most likely a result of only having two tests run at each loading condition.
Large variability in growth rates is common as described in [63], and many more tests would
need to be run to completely quantify the variation associated with the fracture variables
studied.
Methods (1)-(4) oer a good baseline for comparison as they have been fully developed in
the literature. Analytical equations do not oer much insight into how the crack interacts
with the geometry and boundary conditions as the crack propagates through the material.
AFGROW allows for a quick analysis and check of classical problems as well as some custom
problems. Numerical solutions such as the nite element method are required for complex
structures, load transfer and other complex loading and boundary conditions.
59
Table 4.17:
Loading
Condition
Comparison of Life Estimates

Nf ,
cycles
Experimental
Analytical
AFGROW
FEQP
VCCT
XFEM
3 050
2 280
2 220
2 270
2 460
2 750
0.2
8 420
7 540
7 370
7 530
8 150
9 090
0.5
42 500
49 600
48 510
49 490
53 630
59 680
0.67
154 000
179 200
175 300
178 800
193 900
215 600
FEQP is the most accurate of the nite element methods explored here relative to the
analytical solution studied. The FEQP method can be extended to the generalized 2-D and
3-D problems with good expected results.
However, cracks must be grown manually and
remeshed at each increment. This requires either signicant user interaction and time, or an
external program which may have to be customized for each problem.
The signicant user interaction required with manually or scripted remeshing for each increment associated with crack propagation led to a VCCT low-cycle fatigue analysis which
propagated a crack along a predened path without remeshing. Predened paths are acceptable, but in complicated geometry and loading the path is often unknown.
Finally, the XFEM-PN low-cycle fatigue analysis was examined. Using the XFEM-PN allows cracks to grow along a solution dependent path, independent of the mesh given it is
suciently rened. The XFEM-PN gives engineers the resource to evaluate dierent crack
locations and orientations with minimal remodeling.
XFEM-PN will be used in the next
analysis which looks at a Monte Carlo simulation to bound the life estimate.
60
Chapter 5
Monte Carlo Analysis
Analyzing the variability in fatigue fracture life of engineering structures can be dicult
when the structures include complex geometry and complex loading and boundary conditions
combined with fatigue crack growth models t to experimental data. Each of these aspects,
including many others, has a statistical distribution, and investigating how these aect the
structure's life in the presence of a crack may lead to safer, more robust designs.
If no
known analytical solution exists for combining these components, one method to analyze
the distribution on life is the Monte Carlo method. In a Monte Carlo simulation, constants
such as loading magnitude and direction, geometric dimensions, or material properties are
replaced by Probability Distribution Functions (PDF). If the PDF is developed from a normal
distribution, the PDF is dened by two parameters, the mean,
, and standard deviation, ,
of the sample. Random values are selected from each PDF and the analysis is run as many
times as needed to satisfy the law of large numbers and central limit theorem [64] mentioned
in Section 2.3.
The analysis of a Center Cracked Plate (CCP) has a known solution so it stands to be a good
validation tool between the analytical equations and the eXtended Finite Element Method
with Phantom Nodes (XFEM-PN) coupled with a low-cycle fatigue analysis in Abaqus 6-13.3.
To limit the scope of the Monte Carlo analysis, distributions were only used for the Walker
equation coecients,
C0 , m,
and
t in Section 3.2.2.
Additionally, due to the number
of runs required to get a representative distribution on the fatigue fracture life estimate,
Nf ,
only the
R = 0.2
loading condition was analyzed. The following sections describe the
random numbers selection process, the analytical and XFEM-PN Monte Carlo simulations,
and results.
5.1
Random Numbers Selection
Annis concluded in [65] that growth constants cannot be independently selected at random
since they are jointly distributed based on the linear regression t.
Three constants are
present in the Walker equation and each is assumed to be normally distributed about their
61
mean in log-log space where the coecients were t.
n-dimensional,
jointly distributed
PDFs are described by a multivariate normal distribution presented in Eq. (5.1), where
is the mean vector and
is the non-singular symmetric covariance matrix describing the
interaction between constants. Several approaches exist for selecting multivariate normally
distributed values.
x1 based on a
and 2 , and so on
One way is to select
x2 conditionally on x1 , 2 ,
x1 , x2 , . . ., xn1 , n , and n [76]. MATLAB
then nd
normal distribution of
until
xn
and
1 ,
is found conditionally on
has a standard function for generating random
numbers based on a multivariate normal distributions called
mvnrnd()
which was used in
this analysis.
h
f (x) =
1
n/2
(2)
||
(x )T 1 (x )
exp
(5.1)
Eq. (5.2) was transformed into log-log space as shown in Eq. (5.3), which was t with a
least-squares multiple linear regression in MATLAB using the form of Eq. (5.4).
da
C0
=
(K)m
dN
(1 R)m(1)

log10
da
dN
(5.2)

= log10 (C0 ) + m log10 (K) m (1 ) log10 (1 R) +
y = b + m1 x1 + m2 x2 +
where:
da
dN
(5.3)
(5.4)
= log10
= log10 (C0 )
m1 = m
m2 = m (1 )
x1
= log10 (K)
x2
= log10 (1 R)
random error
log10 units. If
were zero, it would imply no interrelation between the variances
The resulting mean vector and covariance matrix are shown in Eq. (5.5) in
all o diagonal values of
on the diagonal; however, the converse is not necessarily true [76]. To determine interdependency, the correlation coecient matrix,
was extracted, where a value of
1 indicates
a direct correlation and zero no correlation. Positive values imply both variances change in
a similar direction, while negative values imply variances change in inverse directions. It is
obvious from Eq. (5.5c) that all coecients in the correlation matrix are correlated since the
o-diagonal terms are signicant relative to 1.0 and must be modeled as such, conrming
Annis' previous conclusion.
62
9.064
= m1 = 3.839
m2
1.674
b2
12 b m1
= 12 b m1
13 b m2
2
m
1
13 b m2 23 m1 m2
(5.5a)
5.680 3.777
3.916
23 m1 m2 = 3.777
2.625 2.365 103
2
m
3.916 2.365
3.812
2
11 12 13
1 0.978
0.842
(5.5b)
= 12 22 23 = 0.978
1 0.748
13 23 33
0.842 0.748
1
(5.5c)
Least-squares linear regression assumes a normally distributed random error term.

describe the distribution of the random error term,
log10 (da/dN )
in Eq.
To
(5.4), a histogram of the
is shown in Fig. 5.1 with a normal PDF superimposed. The histogram shows
a decent comparison with the normal PDF. To conrm that the
log10 (da/dN )
and the nor-
mal PDF have similar distributions, a quantile-quantile, or Q-Q, plot is presented in Fig.
5.2 comparing the two distributions. The reference line in Fig. 5.2 represents both samples
originating from the same distribution.
The bulk of the
log10 (da/dN )
data lie near the
reference line. Data at the left and right ends of the plot deviate from the reference line.
Points that deviate from the reference line on the left end are attributed to the non-linearity
of the crack growth rate approaching
Kth .
Points that deviate from the reference line on
the right end are explained by the non-linearity of the crack growth rate when approaching
Kmax = KIc .
From the Q-Q plot, it can be concluded that the
log10 (da/dN )
data, and
consequently the random error is normally distributed in log-log space.

In the Monte Carlo analyses, random variables used in this study were selected in log-log
space using
and in Eqs.
(5.5a) and (5.5), then converted back to
inverse of Eq. (5.4) and used to estimate
Nf
C0 , m,
and
with the
in the analytical and XFEM-PN simulations
described next.
63
40
35
Observations
30
25
20
15
10
5
0
7
Fig. 5.1:
5
4
3
log10(da/dN) from Experiments
Histogram of
log10 (da/dN )
from experimental data.
10
Quantiles of log (da/dN)
4
5
Exp. Data
Reference
6
4
Fig. 5.2:
Q-Q plot of
1
0
1
2
Standard Normal Quantiles
log10 (da/dN )
from experimental data with a normal PDF.
64
5.2
Analytical Simulation
The mean vector and covariance matrix were used to generate random numbers for the
variables in the Walker equation, which were inserted into Eq.
ai
to
ac .
(5.7) for integration from
A convergence study was performed to determine an adequate sample size for the
Monte Carlo simulation where the number of samples was increased in increments of 1000
until the mean and standard deviation on
Nf
did not change by more than 0.5% between
iterations. This process was repeated ve times resulting in sample sizes ranging from 20007000 samples.
A total number of 10,000 samples were selected since the increase in the
computational time was minimal.
K =
Nc
Nf =
Ni
K
(1 R)1
1
dN =
C0
ac
ai
(5.6)
da

K (a)
(5.7)
m
Results from the analytical Monte Carlo simulation are shown in Fig. 5.3 with a mean value
of
= 7550 cycles
and a standard deviation of
= 240 cycles.
The coecient of variation
(CV) was used to normalize the standard deviation such that CV
= / 100%.
A 3.2%
variation in estimated fatigue fracture life is excellent. Results from Hudson's report [6] are
also shown by the vertical line at the recorded average data point, emphasizing the fact that
only two experimental runs at
R = 0.2
m = 15.0 ksi is not an adequate sample size to

and ts very well to the analysis which was
Distributions of the 10,000 samples of C0 , m, and
and
describe the system. A normal PDF using

conrmed by the Q-Q plot in Fig. 5.4.
are found in Appendix D.
For a comparison on the correlation eects of the covariance matrix, Fig.
5.5 illustrates
what happens to the standard deviation when no correlation is assumed between the Walker
constants. The CV for the uncorrelated analysis increased to 23%, signicantly larger than
the correlated simulation above, thus highlighting the importance of correctly selecting the
random variables in a Monte Carlo simulation.
Additionally, the distribution no longer
appears to be described by a normal PDF, and may more closely be compared to a lognormal PDF.
Crack length and cycle count are plotted in Fig. 5.6 where cycles were calculated from the
10,000 sampled Walker coecients. Life was calculated at several crack increments by integrating Eq. (5.7) from
ai
to the incremental crack length. Mean life and standard deviation
were determined at each crack incremented based on a normal distribution. Recorded data
from Hudson's report and XFEM-PN results from Section 4.6 were included. Neither the
experimental data nor XFEM-PN results fall within three standard deviations of the mean
analytical values.
This gure further demonstrates that the two experiments from Hud-
son's report are not adequate to realize the variability in crack growth. The overestimated
XFEM-PN fatigue cycles are a result of the underestimated energy release rates as described
in Section 4.6.
65
800
700
= 7550
= 240
C.V. = 3.2%
normal
Hudson
Observations
600
500
400
300
200
100
0
6.5
7.5
8.5
Life Estimate, Nf (x10 cycles)
Fig. 5.3:
Monte Carlo simulation of fatigue fracture life with 10,000 samples (R
= 0.2).
9000
Quantiles of Nf
8500
8000
7500
7000
Data
Reference
6500
4
Fig. 5.4:
1
0
1
Q-Q plot of the distribution of Monte Carlo simulation compared to a normal
PDF.
66
1000
= 7730
= 1780
C.V. = 23.1%
Observations
800
600
400
lognormal
normal
Hudson
200
0
10
12
14
16
18
20
22
Life Estimate, N (x10 cycles)
Fig. 5.5:
Monte Carlo simulation of fatigue fracture life with 10,000 samples with no
correlation between Walker constants (R
Fig. 5.6:
= 0.2).
Monte Carlo simulation of fatigue crack growth with 10,000 samples (R
67
= 0.2)
5.3
XFEM-PN Simulation
A similar method was followed for the XFEM-PN Monte Carlo simulation, except running
10,000 samples was not feasible with each nite element analysis requiring approximately
1.0-1.5 hours to complete. A total of 500 runs were completed, each one with a randomly
selected, correlated
converted to
c3
and
C0 , m,
c4 with
and
with the methods described in Section 5.1. Values were
Eq. (5.8) and inserted into a standard Abaqus INP le. Python
code was developed to handle the multiple runs, minimizing user interaction.
c3 = C0 (1 R)
c4 =
m/2
E
1 R2
(5.8a)
m
2
(5.8b)
Results from the XFEM-PN Monte Carlo simulation are shown in Fig.
5.7 with fatigue
fracture life estimated with the integrated polynomial method. A mean fatigue fracture life
of
= 9090 cycles
was found with CV
= 3.3%.
Mean life values for XFEM-PN are 20%
higher than the analytical estimation as a result of Abaqus' predicted

being about 5% too low, as explained in Section 4.6.
KI 's
for the CCP
A Q-Q plot comparing the results
from the XFEM-PN simulation with a normal PDF is provided in Fig. 5.8 showing good
correlation between the two distributions. More importantly, the normal distribution of the
XFEM-PN estimate was similar to the distribution from the analytical method as seen with
the Q-Q plot in Fig. 5.9, verifying that XFEM-PN is an acceptable method for analyzing
the variability in fatigue fracture life. Probability distributions for the Abaqus coecients
c3
and
c4
used in this analysis are found in Appendix D.
68
100
= 9100
= 310
C.V. = 3.4%
Observations
80
60
normal
Hudson
40
20
0
8.2
8.4
8.6
8.8
9.2
9.4
9.6
9.8
10
10.2
Life Estimate, Nf (x10 cycles)
Fig. 5.7:
Monte Carlo simulation of fatigue fracture life estimated with XFEM-PN and the
integrated polynomial method with 500 samples (R
= 0.2).
10500
Quantiles of N
10000
9500
9000
8500
Data
Reference
8000
4
Fig. 5.8:
1
0
1
Q-Q plot of the distribution of the XFEM-PN Monte Carlo simulation compared
to a normal PDF.
69
Quantiles of XFEMPN Nf (103 cycles)
10.5
10
9.5
9
8.5
8
7.5
6.5
Data
Reference
7
7.5
8.5
Quantiles of Analytical Nf (10 cycles)
Fig. 5.9:
Q-Q plot of the distribution of the XFEM-PN Monte Carlo simulation compared
to the distribution of the Analytical Method simulation.
5.4
Conclusion
This chapter looked at two Monte Carlo simulations:
the rst following the traditional
analytical approach allowing for a large number of samples, and a second approach based on
crack growth estimates using the Abaqus XFEM-PN low-cycle fatigue solver, the proposed
method in this thesis of estimating the variability in fatigue fracture life. It was shown that
a multivariate normal distribution was needed to appropriately select random numbers from
the distributions for the random variables in the Walker equation. Crack growth constants
C0 , m,
and
which were all found to be strongly correlated. Both simulation methods gave
similar results, with comparable coecients of variation of approximately 3%. Distributions

of the results appeared to be normal by analyzing Q-Q plots.
An XFEM-PN Monte Carlo simulation has proven to be a viable method for analyzing the
variability in fatigue fracture life estimates for a center cracked plate. It is expected that the
Abaqus XFEM-PN low-cycle fatigue solver could be used as the basis for extending the life
predictions to more complex geometries, loading and boundary conditions.
70
Chapter 6
Conclusions and Future Work
The nal chapter of this thesis summarizes the work performed and conclusions found
throughout the analyses. Recommendations are provided for future work in three-dimensional
fatigue crack propagation analyses. First, the goals of the thesis are reviewed.
6.1
Research Goals and Objectives
Research goals were stated in Chapter 1 with the primary goal to investigate the capabilities
of the eXtended Finite Element Method with Phantom Nodes (XFEM-PN) for modeling
crack propagation and estimating fatigue fracture life.
Then, use XFEM-PN to quantify
the distribution of fatigue fracture life. If this is possible with a simple model, condence
will be gained when using this method with complex models that have not been studied in
the literature, and may not have an analytical description. The objectives that were met to
realize these goals include:
1. Fit fatigue crack growth rate models to experimental data.

2. Analyze crack growth using methods established in the literature.
3. Analyze crack growth using XFEM-PN and compare results to the experimental data
and established methods.
4. Combine the variability in the tted crack growth rate models with the XFEM-PN to
determine the distribution on fatigue fracture life estimates.
6.2
Conclusions
The Walker crack growth model was selected to estimate the fatigue crack growth rates.
Unlike the Paris equation, the Walker equation addresses the eects of the
71
R-ratio,
thus
the entire data set was used to t the crack growth model. The Walker equation was also
selected because Hudson's data set did not contain an adequate number of samples in the
unstable and threshold regions of crack growth to model the Forman equation and Broek and
Rice equation. The Walker equation was easily modied for implementation into Abaqus'
low-cycle fatigue analysis tool used for fatigue crack growth.
Four models were selected to estimate the fatigue fracture life of the Center Cracked Plate
(CCP) for comparison to XFEM-PN: (1) analytical equations, (2) AFGROW, (3) FEQP,
and (4) VCCT. First, integration of the Walker equation was performed from an initial crack
ac , computed from the critical stress intensity factor, Kc .

factor, F , was found in Tada's handbook [2]. This method
length to the critical crack length,

The dimensionless geometry
provided a solid reference for comparison as these equations have been thoroughly developed
in the literature.
AFGROW provided another baseline for comparison. The program, developed by The Air
Force Research Laboratory and LexTech, Inc. [34], uses a database of equations to calculate Stress Intensity Factors (SIFs) for numerous crack and specimen geometry. Cracks are
incrementally grown until reaching
Kc ,
and life is summed between crack increments. The
CCP is a predened model in AFGROW and required only dimensions, material data, and
a load spectrum.
Stress intensity factors (SIFs) from the AFGROW analysis had a Root
Mean Square Error (RMSE) less than 0.5% relative to the analytical analysis.
in the SIFs was attributed to dierent equations for
F.
Dierence
AFGROW's estimate of cycles to
failure for the four loading conditions were approximate 2.5% low.
A Linear Elastic Fracture Mechanics (LEFM) Finite Element formulation using Quarter
Points (FEQP) was developed in the mid-1970's by Henshell and Shaw [45] and Barsoum [46]
as one of the rst nite element methods that did not require the use of specialized elements
and could adequately model the
1/ r
stress singularity at the crack tip.
The method
requires the mid-side nodes of quadratic elements adjacent to the crack tip to be moved
to the quarter point location, and the edges adjacent to the crack tip to be collapsed to a
single point. The Abaqus nite element code only allows quasi-static analyses with quarter
point method, thus the crack must be manually grown and the model remeshed between
increments.
Dowling's [8] modied Simpson's method was used to reduce the number of
models needed to calculate the fatigue fracture life.
Twenty models were selected and a
Python code was developed to grow the crack, remesh the model, and extract pertinent
results from the Abaqus Output Database (ODB) le. The FEQP analyses provided excellent
agreement with the analytical analyses.
K 's throughout propagation had an RMSE of less
than 0.5% and fracture life was within 0.5% of the analytical result. A Crack Propagation
Direction (CPD) analysis with the three criteria available in Abaqus showed that crack
growth would occur perpendicular to the mid-plane of plate in a pure Mode I fashion, as
expected.
Rybicki and Kanninen [50] developed the Virtual Crack Closure Technique (VCCT) in the
late 1970's based on Irwin's assertion that the energy required to grow a crack is equal to
the energy required to close a crack of the same length. VCCT is naturally accommodated
in a nite element framework as the force required to hold two nodes together at a crack
tip and the displacement behind the crack tip are readily available. This method is easier
72
to use than FEQP since it does not require a focused meshed around the crack tip and can
be easily applied to crack growth models if the crack propagation plane is known.
singularities are not modeled in this analysis, and only the energy release rate,
for crack propagation.
SIFs, converted from
a/b < 0.03.
is needed
in the ODB le, compared well with the
analytical analysis with nodal errors were less than 2% for

when
G,
Stress
a/b 0.03,
but were too low
The RMSE for all SIFs was within 1% of the analytical values.
The
low SIFs at the beginning of crack growth led to overestimating the fatigue fracture life
by approximately 8%. The integrated polynomial approach was developed to calculate life
based on Abaqus' results extracted from the ODB le. As a result, a less rened mesh was
required to estimate life, saving signicant computational time.
XFEM-PN was developed by Song, Areias, and Belytschko [4] in 2006 as a method to arbitrarily dene a crack, independent of the mesh, given the mesh is suciently rened. A Level
Set Method (LSM) is used to dene the crack location using signed distance functions stored
as nodal values.
Abaqus combines XFEM-PN, LSM, and a direct cyclic solver propagate
fatigue cracks along a solution dependent path [5]. A

calculate
G ,
modied
VCCT method was used to
and a modied Paris equation determines the cycles between crack lengths.
SIFs had a RMSE 5% lower than the analytical analysis which led to a 20% overestimation
of expected life.
The integrated polynomial method was used in this analysis as well.
CPD analysis showed crack propagation perpendicular to the mid-plane of the plate under
pure Mode I growth, the same as the FEQP solution.
A Monte Carlo simulation was selected as the method to model variability in the fatigue
fracture life estimate. This study focused on random numbers selected from the distributions
of the random variables representing the constants
C0 , m,
and
of the Walker equation.
As suggested by Annis [65] and conrmed in this thesis, the growth constants are highly
correlated and cannot be selected independent of each other, thus a multivariate normal
PDF was used. Two simulations were performed: (1) an analytical simulation allowed the
selection of 10,000 samples to ensure that the law of large numbers and central limit theorem
were satised, and (2) an XFEM-PN simulation. Both the analytical simulation and XFEMPN simulations resulted in normal distributions on life and Coecients of Variation (CV) of
approximately 3%.
These results show that for the analysis of a CCP the XFEM-PN provides a viable alternative
to estimating the fatigue fracture life while providing an understanding of the variability
in that solution.
It is expected that the methods used in the study can be extended to
other crack propagation problems for estimating life under more complex geometry, loading
scenarios and boundary conditions.
Applying the methods developed in this thesis has
great potential, leading to safer, more robust design decisions if fatigue fracture life can be
accounted for earlier in the design phase.
6.3
Future Work
Work performed in this thesis was limited to a center cracked plate made from a 7075T6 Aluminum alloy in two dimensions because experimental data on crack growth and
73
an analytic solution for the LEFM analysis existed in various forms to provide a baseline
for evaluating nite-element-based crack propagation methods. An obvious rst extension
would be modeling crack propagation using XFEM-PN in three dimensions. Some XFEMPN three-dimensional crack growth analyses have already been performed by the author for
other research projects, including a compact specimen used in a crack growth test. Material
properties used in the model are listed in Table E.1 in Appendix E for standard rail steel.
Compact specimen dimensions are provided in Fig.
E.1 in Appendix E. The partitioning
scheme and coarse structured mesh of a half-model with symmetry about the mid-plane of
the specimen are shown in Fig. 6.1.
Compact specimens are loaded by two pins through the top and bottom holes of the specimen.
Loading and boundary conditions at the pins were modeled as reference points at the center
of the pins. The reference points used kinematic coupling constraints for all active degrees-offreedom to distribute the load and boundary conditions to surfaces on the compact specimen.
The surfaces were dened by the top and bottom 45 of the top and bottom pin holes,
respectively. Loading and boundary conditions of each reference point are listed in Table E.2.
Contact between the pins and the compact specimen was explored; however, Abaqus 6-13.3
currently contains a programming error where the direct cyclic solver would not stabilize if
surface-to-surface contact was present. The programming error resulted in incorrect results,
thus analyzing contact was abandoned.
Dassault Systmes Simulia Corporation has been
notied and veried the error. The error is to be corrected in a future release.
The nite element model was built with three dimensional, linear continuum elements (C3D8)
as these are the only brick element types allowed for an XFEM-PN analysis in Abaqus. The
incompatible mode option (C3D8I) was used for the elements outside the enrichment region,
since this option adds degrees of freedom for curvature which aids in convergence when
bending is present.
Before a convergence analysis was performed, a post-processing routine was needed to determine how to extract and analyze the output data from the model. Three dimensional fatigue
K , across the crack front.

Fig. 6.2 shows an example of the variation in the normalized K values for four randomly
selected cycles from the second mesh renement. The reduction in K on the right side
of the gure is a result of the free surface at z/b = 0.5, where z is the through-thickness
coordinate and b is the specimen thickness. The mid-plane of the specimen is at z/b = 0.
For the purposes of this analysis, the median value of the K prole was used to represent
the total K used to calculate the crack growth rate, da/dN .
crack growth analyses have varying stress intensity factor ranges,
The varying
K values also resulted in dierent crack growth rates across the crack front, and
thus non-uniform crack propagation was seen. The non-uniform crack growth was observed
as some fraction of the elements adjacent to the crack front releasing since the elements
met the criteria for fracture.
It often took another solution cycle for remaining elements
to release and to become a uniform crack front again. The primary issue with non-uniform
crack fronts, is that the sharp change in the crack geometry resulted in energy release rates
three to four times higher for elements not fractured compared to the other elements on the
crack front. Once the crack front was uniform, the energy release rates for all elements on
74
Viewport: 1
Model: PressModel
Step: StaticLoad
Viewport: 2
Pin_LD
Y
X
Pin_BC
Pin_BC
Step: Initial
Pin_LD
Model: PressModel
Y
X
Fig. 6.1:
Compact specimen partitioning scheme and coarse mesh.
10
1
17555
22947
26529
0.3
0.4
K tw1/2/ P
6
0
0.1
0.2
z/b
Fig. 6.2:
Variation in the normalized
75
across a compact specimen half model.
0.5
the crack front were similar. Throughout this analysis, results from non-uniform crack fronts
were discarded.
An XFEM-PN convergence analysis was performed and the results at listed in Table 6.1.
Mesh sizes are dened by the number of elements along the crack propagation plane times the
through-thickness elements.
values during propagation where compared to analytical
results for the convergence criteria.
The equation for
is provided in Eq.
(6.1) from
= a/w, a is the crack length, w is the distance from the

the specimen, P is the load range, and b is the specimen
ASTM Standard E647 [66], where

center line of the pins to back of
thickness. The Root Mean Square Dierence (RMSD) is 0.11% between Mesh #1 and #2
and both have a Root Mean Square Error (RMSE) of less than 1%. Models were not run to
complete fracture due to time constraints which resulted in dierent nal crack lengths for
each model. Since the nal crack lengths were not the same, the fatigue cycle counts of the
models were not comparable to each other, thus it was not used for measuring convergence.
Mesh #2 was considered converged for this analysis.

P (2 + )
0.886 + 4.64 13.322 + 14.723 5.64
K =
3/2
b w (1 )
Fig. 6.3 shows the median normalized
and the analytical results.
(6.1)
values at each crack increment during propagation
The XFEM-PN results compare well with an RMSE of 0.67%
and a maximum error of -1.2% occurring at the nal crack increment. Fig. 6.4 shows the
cycle count during crack growth for the analytical analysis, Abaqus' results from the 3-D
compact specimen, and the integrated polynomial method of calculating life from Abaqus'
results. The integrated polynomial method has an RMSE of 1.1% and a maximum error of
1.7% occurring at the nal crack length, while Abaqus' results have a larger error with an
RMSE of 11% and a maximum error of 27% at the nal crack increment. The lower RMSE
again shows that the integrated polynomial method can be applied to various crack growth
analyses. Deformation modes and principal strain elds for the compact specimen are found
in Appendix E.
Results from the compact specimen analysis show that XFEM-PN and the methods discussed
in this thesis can be applied to models other than the center cracked plate, including models
in three dimensions. Some of the issues for future work include optimizing Abaqus' damage
tolerance setting to possibly remove the non-uniform crack growth.
for analyzing the
Also, other methods
through-thickness prole should be explored, including looking at the
Table 6.1:
Convergence Analysis for the Compact Specimen
Elements On
2D Propagation Plane
Max.
RMSE
Max.
RMSD
96
-1.6
0.92
320
-1.2
0.67
0.21
0.11
768
-0.26
0.14
0.20
0.12
Mesh
76
% Error
% Dierence
10
R = 0.2
K tw1/2/ P
7
Analysis
________
XFEMPN
Analytical
6
0.3
0.34
0.38
0.42
0.46
0.5
a/w
Fig. 6.3:
Comparison of the median normalized
KI
during crack propagation in a 3-D
compact specimen for the analytical and XFEM-PN methods (R
= 0.2).
0.6
Analysis
________
0.5
a/w
R = 0.2
Abaqus
Int. Poly
Analytical
0.55
0.45
0.4
0.35
0.3
10
15
20
25
30
Fatigue Cyles, N (103cycles)
Fig. 6.4:
Comparison of cycle count during crack propagation in a 3-D compact specimen
for the analytical and XFEM-PN methods (R
= 0.2).
77
eects the prole has on crack growth projections along the crack front.
These analyses
should be coupled with experimental data to determine if the crack growth patterns are
similar to those seen in the test specimens.
This thesis used a Monte Carlo simulation to quantify the distribution on fatigue fracture
life for the CCP specimen.
Extension of the Monte Carlo method into three-dimensional
crack propagation may be dicult due to the computational time and le size requirements.
For instance, the second mesh renement had
and only propagated a crack
28 900
nodes and
86 700
degrees-of-freedom
0.45 in. (11.4 mm), taking 16 hours on 64-bit personal computer
with 20.0 Gb of RAM, resulting a ODB le size of 3 Gb. Other methods should be explored
for analyzing fatigue life distributions in 3-D including the advanced Uncertainty Quantication method proposed by Bigoni, Engsig-Karup, and True [77]. The advanced Uncertainty
Quantication method required fewer simulations than the traditional Monte Carlo method,
which would be advantageous for 3-D crack propagation analyses.
78
Bibliography
[1] Broek, D., 1988.
The Practical Use of Fracture Mechanics. Kluwer Academic Publishers,
Norwell, MA.
The Stress Analysis of Cracks Hand-
[2] Tada, H., Paris, P., and Irwin, G. R., eds., 1973.
rd
ed. Del Research Corp., Hellertown, Pa, ch. 2, pp. 101305.
, 3
book
See also URL
http://www.abc.edu.
[3] Belytschko, T., and Black, T., 1999.

Minimal Remeshing.
Elastic Crack Growth in Finite Elements with
International Journal for Numerical Methods in Engineering, 45,
pp. 601620.
[4] Song, J.-H., Areias, P. M. A., and Belytschko, T., 2006.
Crack and Shear Band Propagation with Phantom Nodes.
Numerical Methods in Engineering,
67(6),
International Journal for
pp. 868893.
[5] Dassault Systmes Simulia Corp., 2013.
Guide.
A Method for Dynamic
Abaqus 6.13: Abaqus Analysis User's
Providence, RI.
[6] Hudson, C. M., 1969. Eect of Stress Ratio on Fatigue-Crack Growth in 7075-T6 and
2024-T3 Aluminum-Alloy Specimens. NASA Technical Note NASA TN D-5390, Langley
Research Center, Langley Station, Hampton, VA, August.
[7] Broek, D., 1986.
Elementary Engineering Fracture Mechanics.
Dordrecht, Hingham,
MA.
Mechanical Behavior of Materials: Engineering Methods for

Deformation, Fracture, and Fatigue, 3rd ed. Pearson Education, Inc., Upper Saddle
[8] Dowling, N. E., 2007.

River, NJ.
[9] Grith, A., 1920.
The Phenomena of Rupture and Flow in Solids.
Transactions, Vol. 221 of
A, pp. 163198.
Philosophical
[10] Irwin, G., 1956. Onset of Fast Crack Propagation in High Strength Steel and Aluminum
Alloys.
Sagamore Research Conference Proceedings,
[11] Anderson, T., 1995.
2,
pp. 289305.
Fracture Mechanics: Fundamentals and Applications, 2nd
Press, Inc., Boca Raton, NY.
79
ed. CRC
[12] Rice, J., 1968. A Path Independent Integral and the Approximate Analysis of Strain
Concentration by Notches and Cracks.
Journal of Applied Mechanics,
35(2),
pp. 379
386.
[13] Shih, C., and Asaro, R., 1988. Elastic-Plastic Analysis of Cracks on Bimaterial Interfaces: Part I-Small Scale Yielding.
[14] Whitcomb, J. D., 1984.
Journal of Applied Mechanics, pp. 299316.
Analysis of Instability-Related Growth of a Through-Width
Delamination. Report NASA TM-86301, NASA Langley Research Center.

[15] Benzeggagh, M., and Kenane, M., 1996. Measurement of Mixed-Mode Delamination
Fracture Toughness of Unidirectional Glass/Epoxy Composites with Mixed-Mode Bending Apparatus.
Composites Science and Technology,
56(4),
pp. 439449.
[16] Reeder, J. R., 2006. 3D Mixed-Mode Delamination Fracture Criteria - An Experimentalist's Perspective. Report No. 2006004860, NASA Langley Research Center.
[17] Paris, P., Gomez, M., and Anderson, W., 1961. A Rational Analytic Theory of Fatigue.
The Trend in Engineering,
13,
pp. 914.
[18] Paris, P., and F.Erdogan, 1963.
Journal of Basic Engineering,
A Critical Analysis of Crack Propagation Laws.
85(4),
December, pp. 528533.
[19] Walker, K., 1970. The Eects of Stress Ratio During Crack Propagation and Fatigue for
2024-T3 and 7075-T6 Aluminum. Eects of Environment and Complex Load History
on Fatigue Life, ASTM STP 462, American Society for Testing and Materials, pp. 114.
[20] Forman, R., Kearney, V., and Engle, R., 1967. Numerical Analysis of Crack Propagation in Cyclic-Loaded Structures.
89(3),
September,
pp. 459463.
[21] Broek, D., and Rice, R., 2003. Fatigue Crack Growth Properties of Rail Steels. Report
No. DOT-TSC-1076, Battelle Columbus Laboratory for U.S. DOT, FRA.
[22] Farris, T., Keer, L., and Steele, R., 1990. Life Prediction for Unstable Shell Growth in
Rails.
Journal of Engineering for Industry,
112,
pp. 175180.
[23] Hartman, A., and Schijve, J., 1970. The Eects of Environment and Load Frequency
on the Crack Propagation Law for Macro Fatigue Crack Growth in Aluminium Alloys.
Engineering Fracture Mechanics,
1(4),
pp. 615631.
[24] Takeo, Y., and Kiyoshi, S., 1976. The Eect of Frequency on Fatigue Crack Propagation
Rate and Striation Spacing in 2024-T3 Aluminium Alloy and SM-50 Steel.
Fracture Mechanics,
8(1),
Engineering
pp. 8288.
[25] Yokobori, T., and Aizawa, T., 1973. The Inuence of Temperature and Stress Intensity
Factor Upon the Striation Spacing and Fatigue Crack Propagation Rate of Aluminum
Alloy.
International Journal of Fracture,
80
9,
pp. 489491.
[26] Murri, G., Salpekar, S., and O'Brien, T., 1991. Fatigue Delamination Onset Prediction
in Unidirectional Tapered Laminates. In
Composite Materials Fatigue and Fracture,
ASTM STP 1110, T. O'Brien, ed., Vol. 3. American Society for Testing and Materials,
Philadelphia, PA, pp. 312399.
[27] Murri, G., O'Brien, T., and Rousseau, C., 1997. Fatigue Life Methodology for Tapered
Composite Flexbeam Laminates.
53rd
,
American Helicopter Society
Annual Forum
May.
[28] Giannis, S., 2013. Utilising Fracture Mechanics Principles for Predicting the MixedMode Delamination Onset and Growth in Tapered Composite Laminates.
Structures,
102,
Composite
pp. 294 305.
Standard Test Method for Mode I Fatigue Delamination Growth Onset of Unidirectional Fiber-Reinforced Polymer Matrix Composites.
[29] ASTM Standard D6115-97, 2013.
ASTM International, West Conshohocken, PA.

[30] Erdogan, F., and Sih, G., 1963. On the Crack Extension in Plates Under Plane Loading
and Transverse Shear.
85(4),
pp. 519525.
Abaqus 6.13: Abaqus Theory Guide.
Providence, RI.
[32] Hussain, M., Pu, S., and Underwood, J., 1974. Strain Energy Release Rate for a Crack
Under Combined Mode I and Mode II. Fracture Analysis, ASTM STP 560, American
Society for Testing and Materials, pp. 228.
[33] Cotterell, B., and Rice, J., 1980. Slightly Curved or Kinked Cracks. pp. 155169.
[34] AFGROW, 2014. Fracture Mechanics and Fatigue Crack Growth Analysis software tool.
On the WWW, January. URLhttp://www.afgrow.net/.
[35] NASGRO, 2014. Fracture Mechanics and Fatigue Crack Growth Analysis Software. On
the WWW, January. URLhttp://www.swri.org/4org/d18/mateng/matint/nasgro/.
Concepts and
[36] Cook, R. D., Malkus, D. S., Plesha, M. E., and Witt, R. J., 2002.
th
, 4
ed. John Wiley and Sons, Inc., Hoboken,
Applications of Finite Element Analysis

NJ.
[37] Jr., V. W., 1970.
The Finite Element Method for Prediction of Crack Behavior.
Nuclear Engineering and Design,
11(2),
pp. 323332.
[38] Anderson, G. P., Ruggles, V. L., and Stibor, G. S., 1971. Use of Finite Element Computer Programs in Fracture Mechanics.
7(1),
International Journal of Fracture Mechanics,
pp. 6376.
[39] Chan, S., Tuba, I., and Wilson, W., 1970. On the Finite Element Method in Linear
Fracture Mechanics.
81
2(1),
pp. 117.
[40] Byskov, E., 1970. The Calculation of Stress Intensity Factors Using the Finite Element
International Journal of Fracture Mechanics,
Method with Cracked Elements.
6(2),
pp. 159167.
[41] Wilson, W., 1973. Finite Element Methods for Elastic Bodies Containing Cracks. In
Mechanics of Fracture, G. Sih, ed. Noordho International Publishing, Leyden.
[42] Hardy, R. H., 1974. A High-Order Finite Element for Two-Dimensional Crack Problems. PhD Thesis, Georgia Institute of Technology, Atlanta, GA, June.
[43] Holston, A., J., 1976. A Mixed Mode Crack Tip Finite Element.
of Fracture,
12(6),
International Journal
pp. 887899.
[44] Jiang, C., and Cheung, Y., 1995.
A Special Bending Crack Tip Finite Element.
71(1),
pp. 5769.
[45] Henshell, R., and Shaw, K., 1975. Crack Tip Finite Elements are Unnecessary.
national Journal for Numerical Methods,
[46] Barsoum, R., 1976.

Fracture Mechanics.
9(3),
Inter-
pp. 495507.
On the Use of Isoparametric Finite Elements in Linear Elastic
International Journal for Numerical Methods,
10(1),
pp. 2537.
[47] Barsoum, R., 1977. Triangular Quarter Point Elements as Elastic and Perfectly Plastic
Crack Tip Element.
International Journal for Numerical Methods,
[48] Hojfeldt, E., and Ostervig, C. B., 1986.

Shoulder Fillets.
11(1),
pp. 8598.
Fatigue Crack Propagation in Shafts with
25(4),
pp. 421 427.
[49] Alegre, J., and Cuesta, I., 2010. Some Aspects about the Crack Growth FEM Simulations Under Mixed-Mode Loading.
International Journal of Fatigue, 32(7), pp. 1090
1095.
[50] Rybicki, E., and Kanninen, M., 1977. A Finite Element Calculation of Stress Intensity
Factors by a Modied Crack Closure Integral.
9(4),
pp. 931 938.

[51] Shivakumar, K., Tan, P., and J.C. Newman, J., 1988. A Virtual Crack-Closure Technique for Calculating Stress Intensity Factors for Cracked Three Dimensional Bodies.
36,
pp. 43 50.
[52] Mabson, G., Deobald, L., Dopker, B., Hoyt, D., Baylor, J., and Graesser, D., 2007.
th
Fracture Interface Elements for Static and Fatigue Analysis. 16
International Con-
ference on Composite Materials.
[53] Hosseini-Toudeshky, H., Ghaari, M. A., and Mohammadi, B., 2013.
Mixed-Mode
Crack Propagation of Stiened Curved Panels Repaired by Composite Patch Under

Combined Tension and Shear Cyclic Loading.
pp. 344 363.
82
Aerospace Science and Technology, 28(1),
[54] Melenk, J., and Babuka, I.
The Partition of Unity Finite Element Method: Basic
Computer Methods in Applied Mechanics and Engineering,
Theory and Applications.
139.
[55] Mos, N., Dolbow, J., and Belytschko, T., 1999. A Finite Element Method for Crack
Growth without Remeshing.
neering,
46,
International Journal for Numerical Methods in Engi-
pp. 131150.
[56] Osher, S., and Sethian, J. A., 1988.
Fronts Propagating with Curvature-Dependent
Speed: Algorithms Based on Hamilton-Jacobi Formulations.
Physics,
79(1),
Journal of Computational
pp. 12 49.
[57] Stolarska, M., Chopp, D. L., Mos, N., and Belytschko, T., 2001.
Growth by Level Sets in the eXtended Finite Element Method.
for Numerical Methods in Engineering,
51(8),
Modelling Crack
pp. 943960.
[58] Shi, J., Chopp, D., Lua, J., Sukumar, N., and Belytschko, T., 2010. Abaqus Implementation of eXtended Finite Element Method Using a Level Set Representation for
Engineering Fracture
Three-Dimensional Fatigue Crack Growth and Life Predictions.
Mechanics,
77(14),
pp. 2840 2863.
[59] Singh, I., Mishra, B., Bhattacharya, S., and Patil, R., 2012. The Numerical Simulation of Fatigue Crack Growth Using eXtended Finite Element Method.
Journal of Fatigue,
36(1),
International
pp. 109 119.
[60] Pathak, H., Singh, A., and Singh, I. V., 2013. Fatigue Crack Growth Simulations of
3-D Problems Using XFEM.
International Journal of Mechanical Sciences, 76, pp. 112
131.
[61] Hansbo, A., and Hansbo, P. A Finite Element Method for the Simulation of Strong
and Weak Discontinuities in Solid Mechanics.
and Engineering,
193.
Computer Methods in Applied Mechanics
[62] Kucharczyk, P., Sharaf, M., and Mnstermann, S., 2012.
On the Inuence of Steel
Microstructure on Short Crack Growth Under Cyclic Loading.
of Fatigue,
41,
pp. 83 89.
[63] Virkler, D., Hillberry, B., and Goel, P., 1979. The Statistical Nature of Fatigue Crack
Propagation.
Journal of Engineering Materials and Technology,
[64] Dunn, W., and Shultis, J., 2012.
Exploring Monte Carlo Methods.
101,
pp. 148 153.
Academic Press, San
Diego, CA.
[65] Annis, C., 2004. Probabilistic Life Prediction Isn't as Easy as It Looks.
ASTM International,
1(2).
[66] ASTM Standard E647-13e1, 2013.
Crack Growth Rates.
Journal of
Standard Test Method for Measurement of Fatigue
ASTM International, West Conshohocken, PA.
83
Standard Test Method for Linear-Elastic PlaneStrain Fracture Toughness KIc of Metallic Materials. ASTM International, West Con-
[67] ASTM Standard E399-12, 2013.

shohocken, PA.
[68] Liaw, P., Leax, T., and Donald, J., 1987.

Steels.
Acta Metallurgica,
[69] Frost, N., 1962.
35(7),
Fatigue Crack Growth Behavior of 4340
pp. 14151432.
Eect of Mean Stress on the Rate of Growth of Fatigue Cracks in
Sheet Materials.
Journal of Mechanical Engineering Science,
4(1),
pp. 2235.
[70] Wang, C., and Miller, K., 1992. The Eects of Mean and Alternating Shear Stresses
Fatigue I& Fracture of Engineering Materials
on Short Fatigue Crack Growth Rates.
I& Structures,
15(12),
pp. 12231236.
[71] Isida, M., 1971. Eect of Width and Length on Stress Intensity Factors of Internally
International Journal of Fracture
Cracked Plates Under Various Boundary Conditions.
Mechanics,
[72] McGinty,
of
7(3),
R.,
Fatigue
Mercer
pp. 301316.
2011.
Crack
Handbook
Growth
Engineering
Rates
Research
for
to
Damage
Operating
Group,
Warner
Tolerance
Design:
Conditions.
Robins,
GA.
Afgrow
See
Sensitivity
example,
also
URL
http://www.afgrow.net/applications/DTDHandbook/.
[73] Rice, J., 1967. Mechancis of Crack Tip Deformation and Extension by Fatigue. Fatigue Crack Propagation, ASTM STP 415, American Society for Testing and Materials,
pp. 247309.
[74] MatWeb, 2014. Material Property data. See also URL
Guide.
http://www.matweb.com, April.
Abaqus 6.13: Abaqus Keywords Reference
Providence, RI.
Probability and Statistics for
[76] Walpole, R., Myers, R., Myers, S., and Ye, K., 2002.
th
, 7
ed. Prentice Hall, Upper Saddle River, NJ.
Engineers and Scientists
[77] Bigoni, D., Engsig-Karup, A., and True, H., 2013. Modern Uncertainty Quantication
Methods in Railroad Vehicle Dynamics.
Conference.
Rail Transportation Division Fall Technical
[78] Orringer, O., Tang, Y., Gordon, J., Jeong, D., Morris, J., and Perlman, A., 1988.
Crack Propagation Life of Detail Fractures in Rails. Final Report DOT-TSC-FRA-881, Federal Railroad Administration, Washington, DC, October.
[79] Feddersen, C., and Broek, D., 1978. Fatigue Crack Propagation in Rail Steels. Rail
Steels-Developments, Processing, and Use, ASTM STP 644, D. Stone and G. Knupp,
eds., American Society for Testing and Materials, pp. 414429.
[80] Jeong, D., Tang, Y., Orringer, O., and Perlman, A., 1998.
Propagation Analysis of
Transverse Defects Originating at the Lower Gage Corner of Rail. Final Report DOTVNTSC-FRA-98-14, Federal Railroad Administration, Washington, DC, December.
84
Appendix A
Material Properties
Table A.1 lists the elastic material properties of 7075-T6 aluminum alloy as determined in
Hudson's report [6] with a standard tensile test.
Poisson's ratio,
was not provided in
Hudson's report, so a value of 0.33 listed on www.matweb.com [74] was assumed. Table A.2
lists the results from fatigue crack growth tests, listing the mean stress,
stress,
a ,
m ,
and alternating
along with the cycles required to grow a crack from an initial half length of
0.1 in.
to the half length listed in each column. Table A.3 shows the results from the residual static
strength test.
Using the nal load,
Pf ,
and the nal crack length,
af ,
the critical stress
intensity factor can be calculated for the material and geometry. A detailed drawing of the
center crack plate specimen is shown in Fig. A.1.
Table A.1:
Average Tensile Properties for 7075-T6 Aluminum
Alloy
yld (1)
ult
(2)
Percent
(ksi)
(MPa)
(ksi)
(MPa)
(ksi)
(GPa)
Elongation
83.2
574
75.9
523
10 100
69.6
12
Source: [6]
(1)
(2)
0.2% Oset
2-inch (51-mm) gage length
85
86
172
138
103
69
34
207
172
138
103
207
172
138
103
69
207
172
138
103
69
207
172
138
103
69
207
172
25.0
20.0
15.0
10.0
5.0
30.0
25.0
20.0
15.0
30.0
25.0
20.0
15.0
10.0
30.0
25.0
20.0
15.0
10.0
30.0
25.0
20.0
15.0
10.0
30.0
25.0
Source: [6]
(MPa)
(ksi)
2.8
3.0
1.7
3.0
3.0
4.4
5.0
3.0
5.0
6.7
8.3
10.0
5.0
7.5
10.0
12.5
15.0
10.0
13.3
16.7
20.0
5.0
10.0
15.0
20.0
25.0
(ksi)
19
21
12
21
21
30
34
21
34
46
57
69
34
52
69
86
103
69
92
115
138
34
69
103
138
172
(MPa)
0.8
0.8
0.7
0.7
0.7
0.7
0.7
0.5
0.5
0.5
0.5
0.5
0.3
0.3
0.3
0.3
0.3
0.2
0.2
0.2
0.2
13.7
13.7
13.7
13.7
13.7
13.7
13.7
30.0
13.7
13.7
0.5
0.5,13.7
30.0
13.7
0.5,13.7
0.5
0.5,13.7
0.5
0.5
0.5
0.25
30.0
20.0
0.5,13.7
0.5
0.5
(Hz)
Freq
71 000
42 500
780 000
74 000
69 000
19 800
12 300
105 000
18 500
9 500
4 100
2 230
16 100
7 350
3 040
1 500
810
3 800
1 680
760
410
70 000
6 140
1 670
700
330
0.2 (5.08)
Table A.2:
94 000
58 750
1 260 000
98 000
97 000
27 000
16 300
140 000
28 000
12 400
5 600
3 040
23 700
10 100
4 260
2 060
1 050
5 650
2 340
990
490
88 300
8 800
2 360
930
413
0.3 (7.62)
106 000
68 750
1 540 000
115 000
112 000
30 900
18 300
160 000
32 800
13 820
6 380
3 430
28 900
11 700
4 900
2 320
1 130
6 550
2 640
1 100
98 300
10 370
2 660
1 030
0.4 (10.2)
111 000
75 000
1 690 000
125 000
121 000
33 000
19 400
175 000
35 600
14 730
6 800
3 660
32 600
12 800
5 280
2 470
1 170
7 150
2 830
1 150
105 000
11 300
2 820
1 075
0.5 (12.7)
116 000
77 500
1 780 000
132 000
127 000
34 200
20 200
185 000
37 600
15 400
7 100
3 810
35 500
5 550
2 570
7 500
2 960
1 180
110 000
11 800
2 930
1 100
0.6 (15.2)
120 000
80 000
1 825 000
137 000
130 000
35 200
20 700
192 500
39 000
7 260
3 920
37 800
5 720
2 630
7 750
3 060
113 000
12 300
3 000
1 110
0.7 (17.8)
122 000
82 000
1 850 000
141 000
134 000
36 000
21 000
197 500
40 000
7 380
3 990
39 500
5 850
2 680
7 950
3 120
117 000
12 570
3 050
0.8 (20.3)
123 000
83 750
1 870 000
144 000
136 000
36 500
202 500
40 800
7 480
40 800
5 940
8 080
120 000
12 800
0.9 (22.9)
Average Number of Cycles Required to Propagate a Crack from a Half-Length,
12 950
124 000
84 000
1 885 000
147 000
138 000
205 000
41 400
7 500
42 000
8 190
122 000
126 000
85 000
1 915 000
150 000
141 000
212 000
42 000
43 500
8 340
127 000
13 150
1.2 (30.5)
1 930 000
152 000
143 000
215 000
42 500
44 500
8 420
129 000
1.4 (35.6)
1 935 000
153 000
144 000
217 500
131 000
1.6 (40.6)
of 0.10 in. to a Half-Length of _ in. (mm)
1.0 (25.4)
a,
Fatigue Crack Growth Test Results for 7075-T6 Aluminum Alloy
154 000
132 000
1.8 (45.7)
Table A.3:
Residual Static Strength Test Results for 7075-T6 Aluminum Alloy
Observation
af
Pf
(in.)
(mm)
(kips)
(kN)
1.85
46.99
29.9
133
1.80
45.72
30.4
135
1.79
45.47
30.4
135
2.36
59.94
27.2
121
1.16
29.46
38.5
171
1.19
30.23
38.8
173
0.96
24.38
43.1
192
1.09
27.69
41.7
185
1.15
29.21
37.7
168
10
0.94
23.88
44.5
198
11
1.75
44.45
32.0
142
12
2.19
55.63
25.0
111
13
0.73
18.54
50.4
224
14
1.53
38.86
35.9
160
15
2.10
53.34
30.2
134
16
1.73
43.94
32.3
144
17
1.43
36.32
36.4
162
18
0.71
18.03
49.8
222
19
1.65
41.91
39.0
173
20
2.73
69.34
24.3
108
21
2.81
71.37
22.8
101
22
2.13
54.10
27.1
121
23
0.93
23.62
47.1
210
24
0.78
19.81
48.4
215
25
1.20
30.48
39.0
173
26
1.37
34.80
36.0
160
27
0.72
18.29
47.5
211
Source: [6]
87
Fig. A.1:
Detailed center crack plate specimen
88
Appendix B
Center Cracked Plate Analysis Results
The following gures show additional results from the analyses performed in Chapter 4.
Figs.
B.1-B.4 plot the normalized
half-width,
a/b,
KI
versus the ratio of crack half-length to specimen
for the loading conditions listed in Table B.1.
(RMSE) for the normalized
KI
Root Mean Square Error
curves compared to the analytical method are provided in
Table B.2. Cycle count during propagation is shown in Figs. B.5-B.8 with RMSE values for
each loading condition and analysis method in Table B.3. Crack propagation direction for
the three criteria in an FEQP analysis, as well as Maximum Tangential Stress (MTS) for
XFEM-PN are plotted in Figs. B.9-B.12.
Resulting nite element model deformations modes and principal strains are provided for each
of the nite element analyses. Deformations modes show the zero displacement boundary
conditions along the axis of symmetry in the
direction and along the bottom in the
direction. Maximum and minimum principal strains,
max
and
min ,
were provided as they
are independent of coordinate system. All gures depict the critical crack length,
ac , denoting
fracture is imminent. Overall specimens are shown on the left of each gure with insets at
the crack tip on the right.
contours are unaveraged.
All images have a deformation scale factor of 40 and color
Two-dimensional quadratic quadrilateral elements in the FEQP
Figs. B.13-B.16 result in smooth contours, as well as signicantly higher max at the crack tip
since mid-side nodes were moved to the quarter point location. Both VCCT in Figs. B.17B.20 and XFEM-PN in Figs. B.21-B.24 use two-dimensional linear quadrilateral elements.
Discretization of the elements is more noticeable in the XFEM-PN model since a courser
mesh was used compared to the VCCT model. Neither VCCT nor XFEM-PN are able to
adequately model the singularity at the crack tip as expected since these methods were not
developed for this purpose.
89
Table B.1:
Loading
max
Frequency
Condition
(ksi)
(MPa)
(ksi)
(MPa)
(ksi)
(MPa)
(ksi)
(MPa)
(Hz)
15.0
103
15.0
103
30.0
207
30.0
207
15.0
103
10.0
68.9
25.0
172
20.0
138
0.2
0.5
15.0
103
5.00
34.5
20.0
138
10.0
68.9
0.5
13.7
15.0
103
3.00
20.7
18.0
124
6.00
41.4
0.67
13.7
0.5/13.7*
* Indicates multiple tests run a dierent frequencies.
Table B.2:
RMSE Compared to the Analytical Method for the Normalized
KI
Curves in
Figs. B.1-B.4.
Loading
Condition
Table B.3:
KI
RMSE Compared to the Analytical Method (%)
AFGROW
FEQP
VCCT
XFEM
0.28
0.23
1.0
4.5
0.2
0.38
0.32
0.86
4.4
0.5
0.48
0.41
0.80
4.3
0.67
0.49
0.43
0.80
4.3
RMSE Compared to the Analytical Method for the Crack Growth Curves in
Figs. B.5-B.8.
Loading
Condition
RMSE Compared to the Analytical Method (%)
Experiments
AFGROW
FEQP
VCCT
XFEM
44
2.2
0.08
8.6
20.3
0.0
13
2.1
0.09
8.6
20.2
0.0
13
2.1
0.12
8.6
20.2
0.00
13
2.0
0.13
8.6
20.2
90
1.1
K/(( a)1/2)
Analysis
________
Analytical
AFGROW
FEQP
VCCT
XFEMPN
1.05
R=0
0.95
0.05
0.1
0.15
0.2
0.25
0.3
0.35
a/b
Fig. B.1:
Normalized
KI
= 0).
1.1
K/(( a)1/2)
Analysis
________
Analytical
AFGROW
FEQP
VCCT
XFEMPN
1.05
R = 0.2
0.95
0.05
Fig. B.2:
0.1
Normalized
0.15
KI
0.2
a/b
0.25
0.3
91
0.35
= 0.2).
0.4
Analysis
________
Analytical
AFGROW
FEQP
VCCT
XFEMPN
K/(( a)1/2)
1.15
1.1
R = 0.5
1.05
0.95
0.1
Fig. B.3:
0.2
Normalized
KI
0.3
a/b
0.4
0.5
0.6
= 0.5).
1.25
Analysis
________
K/(( a)1/2)
1.2
Analytical
AFGROW
FEQP
VCCT
XFEMPN
1.15
R = 0.67
1.1
1.05
1
0.95
0.1
Fig. B.4:
0.2
Normalized
KI
0.3
a/b
0.4
0.5
92
= 0.67).
0.6
0.35
Analysis
________
0.3
R=0
= 30.0 ksi
(1042 MPa)
Hudson
Analytical
AFGROW
FEQP
VCCT
XFEMPN
0.25
a/b
0.2
0.15
0.1
0.05
0
0.5
1.5
2.5
3.5
Fig. B.5:
= 0).
0.4
Analysis
________
0.35
0.3
0.25
a/b
R = 0.2
= 20.0 ksi
(695 MPa)
Hudson
Analytical
AFGROW
FEQP
VCCT
XFEMPN
0.2
0.15
0.1
0.05
0
Fig. B.6:
93
= 0.2).
10
0.6
Analysis
________
0.4
a/b
R = 0.5
= 10.0 ksi
(347 MPa)
Hudson
Analytical
AFGROW
FEQP
VCCT
XFEMPN
0.5
0.3
0.2
0.1
0
10
20
30
40
50
60
Fig. B.7:
= 0.5).
0.6
Analysis
________
Hudson
Analytical
AFGROW
FEQP
VCCT
XFEMPN
0.5
a/b
0.4
0.3
R = 0.67
= 6.0 ksi
(208 MPa)
0.2
0.1
0
50
100
150
200
Fig. B.8:
94
= 0.67).
250
0.15
0.1
R=0
0.05
0
CPD
Method
__________
0.05
MTS
MERR
K =0
0.1
II
XFEMPN
0.15
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
a/b
Fig. B.9:
= 0).
0.15
0.1
R = 0.2
0.05
0
CPD
Method
__________
0.05
MTS
MERR
KII=0
0.1
XFEMPN
0.15
0.05
Fig. B.10:
0.1
0.15
0.2
0.25
a/b
0.3
0.35
0.4
0.45
95
0.5
= 0.2).
0.15
0.1
R = 0.5
0.05
0
CPD
Method
__________
0.05
MTS
MERR
K =0
0.1
II
XFEMPN
0.15
0.1
0.2
0.3
0.4
0.5
0.6
0.7
a/b
Fig. B.11:
= 0.5).
0.15
0.1
R = 0.67
0.05
0
CPD
Method
__________
0.05
MTS
MERR
KII=0
0.1
XFEMPN
0.15
0.1
Fig. B.12:
0.2
0.3
0.4
a/b
0.5
0.6
0.7
96
0.8
= 0.67).
Fig. B.13:
FEQP, deformation mode
ux
with
a = 2.24 in. (R = 0.2).
Fig. B.14:
FEQP, deformation mode
uy
with
a = 2.24 in. (R = 0.2).
97
Fig. B.15:
FEQP, principal strain
max
with
a = 2.24 in. (R = 0.2).
Fig. B.16:
FEQP, principal strain
min
with
a = 2.24 in. (R = 0.2).
98
Fig. B.17:
VCCT, deformation mode
ux
with
a = 2.25 in. (R = 0.2).
Fig. B.18:
VCCT, deformation mode
uy
with
a = 2.25 in. (R = 0.2).
99
Fig. B.19:
VCCT, principal strain
max
with
a = 2.25 in. (R = 0.2).
Fig. B.20:
VCCT, principal strain
min
with
a = 2.25 in. (R = 0.2).
100
Fig. B.21:
XFEM-PN, deformation mode
ux
with
a = 2.38 in. (R = 0.2).
Fig. B.22:
XFEM-PN, deformation mode
uy
with
a = 2.38 in. (R = 0.2).
101
Fig. B.23:
Fig. B.24:
XFEM-PN, principal strain
XFEM-PN, principal strain
102
mAX
min
with
with
a = 2.38 in. (R = 0.2).
a = 2.38 in. (R = 0.2).
Appendix C
Convergence Analysis for Finite Element
Models
The following gures show additional convergence analysis results for each of the nite element models. For the Finite Element with Quarter Points (FEQP) analyses, Figs. C.1-C.3
depict the normalized
re .
KI
variation at contour integrals within the focused mesh region,
Radial distances to the contour,
r,
are measured from the crack tip. Root Mean Square
Error (RMSE) of all contours in each mesh are listed in Table C.1 with errors relative to the
analytical solution.
Typical results of the stabilization ratios
R0 , Rn , U0 , and Un
for both VCCT and XFEM-PN
analyses are shown in Fig. C.4. Stabilized cycle ratios were plotted versus the half crack
length,
a,
normalized with half width of the plate,
b,
for each of the mesh density. All ratios
R0 in the VCCT analysis.

R0 increases with crack length and mesh density, requiring a user dened
critical value for CR0 as noted in Section 4.5. Spikes seen in Meshes #3 and #4 were caused
were signicantly smaller than the default value of 0.005, except

Fig. C.5 shows that
when the analyses were terminated due to single cycle increments, also noted in Section 4.5.
After termination, the analyses were restarted and allowed to run until the critical energy
release rate was reached.
To build gures like Figs. C.4 and C.5, Fourier terms required to calculate the stabilization
ratios were extracted from Abaqus' message (MSG) les for each iteration within each cycle.
Fig. C.6 shows how
R0
converges to a stabilized solution for four arbitrarily chosen cycles in
the Mesh #1 VCCT analysis. This behavior was typical of all stabilization ratios for both
VCCT and XFEM-PN analyses. Once the cycle had stabilized, the nal ratio values were
selected and plotted in Figs. C.4 and C.5. For example, the nal data point is marked for
cycle number 2253.
103
Table C.1:
KI
Mesh
RMSE of
KI
Contour Values in Convergence Analysis.
Contour Value RMSE Relative to Analytical Results (%)
a = 0.1 in.
(2.54 mm)
a = 1.0 in.
(25.4 mm)
a = 3.0 in.
(76.2 mm)
0.20
0.33
0.97
0.13
0.31
0.96
0.10
0.31
0.96
1.02
Mesh
____
#1
#2
#3
Analytical
K/(( a)1/2)
1.015
1.01
1.005
a/b = 0.02
1
0.995
0.99
Fig. C.1:
0.1
0.2
FEQP normalized
0.3
KI
0.4
0.5
r/re
0.6
0.7
0.8
0.9
a = 0.1 in.
(2.54 mm)
1.02
a/b = 0.17
K/(( a)1/2)
1.019
1.018
1.017
Mesh
____
#1
#2
#3
Analytical
1.016
1.015
1.014
1.013
Fig. C.2:
0.1
0.2
FEQP normalized
0.3
KI
0.4
0.5
r/re
0.6
0.7
0.8
104
0.9
a = 1.0 in.
(25.4 mm)
K/(( a)1/2)
1.19
1.185
Mesh
____
#1
#2
#3
Analytical
a/b = 0.50
1.18
1.175
Fig. C.3:
0.1
0.2
FEQP normalized
0.3
KI
0.4
0.5
r/re
0.6
0.7
0.8
0.9
a = 3.0 in.
(76.2 mm).
1.2
Typ. Stabilization Ratio 10
15
Mesh
____
#1
#2
#3
#4
0.8
0.6
0.4
0.2
0.05
0.1
0.15
0.2
0.25
0.3
0.35
a/b
Fig. C.4:
Typical stabilization ratio for VCCT and XFEM-PN direct cyclic analyses during
fatigue crack propagation.
105
Constant Residual Ratio, R0
0.25
0.2
Mesh
____ #1
0.15
#2
#3
0.1
#4
0.05
0.05
0.1
0.15
0.2
0.25
0.3
0.35
a/b
Fig. C.5:
R0 , during fatigue crack propagation.
5
Constant Residual Ratio, R0
Cycle
_____
4
Fig. C.6:
R0
2253
3066
3491
3754
Cyc: 2253
Iter: 6
R0: 0.0183
3
4
Iteration Number within Cycle
labeled once cycle 2253 has been stabilized.
106
R0 ,
during four cycle iterations with
Appendix D
Monte Carlo Analysis Results
Distributions for the correlated Walker constants
C0 , m,
and
of Eq. (2.11) are provided in
Figs. D.1-D.3 for the 10,000 samples used in the analytical Monte Carlo analysis. Figs. D.4
and D.5 show the distributions of the 500 samples of
c3
and
c4
used in Abaqus' interpretation
of the Paris law, Eq. (4.19). A normal Probability Density Function (PDF) is superimposed
on the
m, ,
and
c4
histograms, and a log-normal PDF is superimposed on the
histograms. Each normal PDF was generated based on the mean,
of the data.
C0
and
c3
, and standard deviation,
All PDFs t well as expected since the correlated random numbers were
generated from a multivariate normal distribution in log-log space.
Walker Coefficient, C0 (x1015 mm/cycle/(MPamm1/2)m)

700
15
20
25
30
35
600
45
50
55
60
= 8.76e10 (2.70e14)
= 1.56e10 (4.82e15)
C.V. = 17.8%
500
Observations
40
400
300
200
100
0
0.4
Fig. D.1:
0.6
0.8
1
1.2
1.4
1.6
1.8
Walker Coefficient, C0 (x109 in./cycle/(ksiin1/2)m)
C0 , with 10,000 samples and a log-normal PDF.
107
700
600
= 3.84
= 0.05
C.V. = 1.3%
Observations
500
400
300
200
100
0
3.6
Fig. D.2:
3.65
3.7
3.75
3.8
3.85
3.9
3.95
Walker Coefficient, m (dimensionless)
m,
4.05
700
600
= 0.564
= 0.012
C.V. = 2.2%
Observations
500
400
300
200
100
0
0.51
Fig. D.3:
0.52
0.53
0.54 0.55 0.56 0.57 0.58 0.59

Walker Coefficient, (dimensionless)
108
0.6
0.61
Abaqus Coefficient, c (x109 mm/cycle/(MPamm)m/2)

3
100
650
700
750
800
850
900
1000
1050
= 2.80e02 (2.72e14)
= 2.07e03 (6.38e08)
C.V. = 7.4%
80
Observations
950
60
40
20
0
20
Fig. D.4:
25
30
Abaqus Coefficient, c3 (x103 in./cycle/(ksiin)m/2)
c3 ,
35
with 500 samples and a log-normal PDF.
100
= 1.92
= 0.03
C.V. = 1.4%
Observations
80
60
40
20
0
1.82
Fig. D.5:
1.84
1.86
1.88
1.9
1.92
1.94
1.96
Abaqus Coefficient, c4 (dimensionless)
109
c4 ,
1.98
with 500 samples and a normal PDF.
Appendix E
Compact Specimen Analysis Results
This appendix provides data used to model the three-dimensional crack growth in a compact
specimen. Table E.1 lists material data taken from Orringer et al. [78], Feddersen and Broek
[79], and supplemented by medium carbon steel properties found on www.matweb.com [74]
for those properties not found in the literature. Jeong, Tang, Orringer, and Perlman [80]
used a form of the Walker equation shown in Eq. (E.1). By comparing Jeong et al.'s equation
with the Walker equation, the Walker constant
constants
and
Q.
was solved for in terms of the empirical
Loading and boundary conditions on the reference points representing
the loading pins are listed in Table E.2. The loading applied to the top reference pin was
reduced by half since only half of the compact specimen was modeled. Required dimensions
are listed in Fig. E.1.
(K)m
da
= C0
dN
(1 R)Q
=1
The resulting deformation modes,
ux , uy ,
and
(E.1)
Q
m
uz
(E.2)
are presented in Figs. E.2, E.3, and E.4,
respectively. The maximum and minimum principal strains,

principal values are independent of coordinate system.
min ,
max
Figs.
and
min
are provided since
E.5 and E.7 show
respectively, for the entire compact specimen. Figs. E.6 and E.8 show
max
max
and
and
min ,
respectively, with the top portion of the compact specimen removed and zoomed in on the
crack tip. These last two gures illustrate the through-thickness variation of strain at the
crack tip. The jump in red contours in Fig. E.6 illustrate a non-uniform crack tip extension.
110
Table E.1:
Rail Steel Material
Properties
Rail Steel Material Properties
29 600 ksi
0.293
1 1011
C0
[8]
in./cycle
(ksi
1.63
in.
m/2
[80]
[80]
[80]
(1)
0.592
Kc
73.7 ksi
[74]
ult
133 ksi
[78]
c3
5.6 103
(2)
c4
(2)
(1)
(2)
Table E.2:
[79]
Calculated from Eq. (E.2)

Calculated from Eq. (4.19)
Compact Specimen Loading and Boundary Conditions
Location
Top Pin Reference Point
Bottom Pin Reference Point
Symmetry Plane
Boundary Conditions
ux = 0
Loading Conditions
P = 2.5
ux = uy = 0
uz = 0
111
kips
Fig. E.1:
Fig. E.2:
Compact specimen dimensions.
112
ux
with
a = 1.46 in. (R = 0.2).
Fig. E.3:
uy
with
a = 1.46 in. (R = 0.2).
Fig. E.4:
uz
with
a = 1.46 in. (R = 0.2).
113
Fig. E.5:
Fig. E.6:
max
with
114
max
a = 1.46 in. (R = 0.2).
with
a = 1.46 in. (R = 0.2).
Fig. E.7:
Fig. E.8:
min
with
115
min
a = 1.46 in. (R = 0.2).
with
a = 1.46 in. (R = 0.2).
Appendix F
MATLAB Code
The following code was developed in MATLAB verions R2013a and R2013b.
Function: CenterCrackF.m
function [ F] = CenterCrackF (a , b)
% Calculates dimensionless geometry factor , F , for a center
% craked plate based on equations from Figure 8.12 in Dowling pg
327
%
% a : crack half length
% b : specimen half width
%
% Rev
Date
Name
Changes
% 00
05 -29 -12
Josh Melson
- Created
%
alpha = a/b ;
F = (1 - 0.5* alpha + 0.326* alpha .^2) ./ sqrt (1 - alpha );
116
Function: CenterCrackFinalF.m
function [F , ac , alpha , fval , flag ] = CenterCrackFinalF ( Pmax , a , b ,

t , KIc )
% Created by Josh Melson 05 -29 -12 , modified 08 -01 -12
% Function solves for the final crack length based on KIc of a
center crack
% plate , equations from Figure 8.12 , Dowling pg 327
% x (1) : F
[ dimless ]
% x (2) : ac
[ length ]
% x (3) : alpha
[ dimless ]
%
% Rev
Date
Name
Changes
% 00
05 -29 -12
Josh Melson
- Created
% 01
08 -01 -12
Josh Melson
- Adjusted solving process
%
Smax = Pmax /(2* b*t) ;
% force / length ^2
G = @( x) [
x (1) - (1 - 0.5* x (3) + 0.326* x (3) ^2) / sqrt (1 - x (3) ) ;

KIc - x (1) * Smax * sqrt ( pi *x (2) ) ;
x (3) - x (2) /b ];
x0 = [ 1.12; a; a/b ];
option = optimset ( ' Display ' ,' off ') ;
[x , fval , flag ] = fsolve (G , x0 , option );
F = x (1) ;
ac = x (2) ;
alpha = x (3) ;
% length
117
Function: dadN_Regression.m (Use with experimental data)
function [ ai , dadN ] = dadN_Regression (a , N , pt )

% Calculates da / dN using a quadratic fit to data points of a cycles
% vs crack length curve
%
% Input :
a: Crack Length
%
N : Cycles
%
pt : Total number of points to use
% Output :
ai : Crack length evaluated from regression curve at
cycle Ni
%
dadN : Slope of regression curve at cycle Ni
%
% Rev
Date
Name
Changes
% 00
11 -06 -13
Josh Melson
- Created
%
%% Check for Errors
if nargin < 3;
error ( ' Not enough inputs ');
end
if nargin > 3;
error ( ' Too many inputs ');
end
if rem ( pt ,1) ~= 0; error ( ' Number of points must be an integer ');
end
if pt > 3;
warning ( ' Suggested points less than 4 ');
end
%% Initialize variables
ai
= zeros ( length (a ) -2* pt ,1) ;
dadN = zeros ( length (a) -2* pt ,1) ;
if length ( ai ) > length ( a); return ; end
%% Adjust magnitude of N to be silimar to a for a better fit
Nord = order ( N (2: end )) ;
Ordr = ceil ( mean ( Nord ) );
Adj = 10^ Ordr ;
Nadj = N/ Adj ;
%% Fit Points
for i = pt +1: length (N) -pt
NT = Nadj (i - pt :i+ pt ) ;
% Reduced set
of N
aT = a(i - pt :i + pt ) ;
% Reduced set
of a
b = polyfit ( NT ,aT ,2) ;
% Quad fit
db = polyder ( b);
% Derivative
Quad fit
ai (i - pt ) = polyval (b , Nadj ( i)) ;
% Crack length
at Ni
dadN (i - pt ) = polyval (db , Nadj (i) )/ Adj ;
% da / dN at Ni
end
118
Function: LifeIntegral_Walker.m (Use with Analytical model)
function [ dNda ] = LifeIntegral_Walker (a , b , deltaS , R , C0 , m , gamma

)
% Function used for integral of Walker fatigue crack growth model
% a : crack length , can be a vector
% b : specimen half width
% deltaS : stress range
% R : R - Ratio
% C0 , m , gamma : Walker coefficients
% dNda : inverse of the da / dN for integration of N from ai to ac
%
% Rev
Date
Name
Changes
% 00
10 -28 -13
Josh Melson
- Created
%
F = CenterCrackF (a , b);
C1 = C0 /(1 - R) ^( m *(1 - gamma ));
dNda = 1 ./ ( C1 *( F* deltaS .* sqrt ( pi *a )) .^ m) ;
119
Function: MonteCarlo_Walker.m (Use with Analytical model)
function [N , C0 , m , gam ] = MonteCarlo_Walker ( coeff , covar , num , ai ,

ac , b , dS , R)
% Monte Carlo simulation using Walker coefficients and covariance
matrix
%
% coeff : log space fit coefficients vector
% covar : log space covariance matrix
% num :
number of observations
%
% Fit based on Walker equation
% da / dN = C0 /(1 - R) ^( m *(1 - gamma )) * DeltaK ^m
%
% Linearize equation
% log ( da / dN ) = log ( C0 ) + m* log ( DeltaK ) - m *(1 - gamma )* log (1 - R)
%
Y = b + m1 * X1 + m2 * X2
%
% << log base 10 > >
%
% Input to fit :
%
X1 : log ( delK )
%
X2 : log (1 - R)
%
Y : log ( dadN )
% Output from fit into Monte Carlo :
%
coeff (1) : b , intercept
%
coeff (2) : m1 , slope1
%
coeff (3) : m2 , slope2
% Conversion to untransformed space
%
C0 = 10^ b
%
m = m1
%
gamma = 1 + m2 / m1
%
% Rev
Date
Name
Changes
% 00
03 -04 -14
Josh Melson
- Created
%
% covar = diag ( diag ( covar ) );
% Use for uncorrelated results
rndcoeff = mvnrnd ( coeff , covar , num ); % Correlated random numbers
C0
m
gam
N
=
=
=
=
10.^ rndcoeff (: ,1) ;

rndcoeff (: ,2) ;
1 + rndcoeff (: ,3) ./ rndcoeff (: ,2) ;
zeros ( size ( C0 ) );
for j = 1: num
Nfun = @(a ) LifeIntegral_Walker (a , b , dS , R , C0 (j) , m( j) , gam (j
)) ;
120
end
N( j) = integral ( Nfun , ai , ac );
121
Script: Abaqus_CrackSizeAnalysis.m (Use with FEQP)
% Creates data file for Python to read and Abaqus to execute

% Modified for Center Crack Plate , Hudson 1996 , 7075 - T6
%
% Rev
Date
Name
Changes
% 00
01 -15 -14
Josh Melson
- Created
%
close all ; clear all ; clc ;
format shortG
%% Load Data
dataPath = 'C :\ Users \ jmelson \ Documents \ Thesis \ ';
load ( strcat ( dataPath , ' Al_7075T6_Hudson '))
load ( strcat ( dataPath , ' Al_7075T6_Hudson_Dimensions '))
%% Define Models to Analyze
Rfld = { 'R0 ';
' R2 ';
' R5 ';
' R7 ' };
Sfld = { ' Sm15 '; ' Sm15 '; ' Sm15 '; ' Sm15 '};
Meth = { ' Walker '};
%% Set Number of Steps
j = 20;
% Must be even integer
%% Define CPU Properties

Memory = ' 18000 mb ';
CPUs
= '4 ';
%% Analyze in Abaqus
cmd1Val = ' abaqus python CenterCrack_CFE_AdjustValues . py % s %s %s ';
for m = 1: length ( Rfld )
S = Al .( Rfld { m }) .( Sfld {m }) ;
% Reduce calling
structure
S. CFE = struct ;
r = (S . ac / ai ) ^(1/ j) ;
jj = 1: j +1;
S. CFE .a = ai * r .^( jj -1) ;
S. CFE .a = S. CFE .a ';
S. CFE . K1
S. CFE . K2
S. CFE . JKs
S. CFE .J
Rvalue
pLoad
=
=
=
=
=
=
zeros ( length (S . CFE .a ) , 1) ;

zeros ( length (S . CFE .a ) , 1) ;
zeros ( length ( S. CFE . a) , 1) ;
zeros ( length (S. CFE .a) , 1) ;
Rfld {m };
strcat ( '-' , num2str (S. dS ) );
for n = 1: length (S. CFE .a)

122
% Current R
% Current Load
aLen
= num2str ( S. CFE . a(n) );
if length ( aLen ) < 5;

else q = 5; end
% Current Length
q = length ( aLen ) ;
jobName = strcat ( ' CFE_ ', Rvalue , '_S ', num2str ( round (S. dS )) ,'
_C ', aLen (1) , aLen (3: q)) ;
cmd1
= sprintf ( cmd1Val , Rvalue , pLoad , aLen );
cmd2
= sprintf ( ' abaqus cae noGUI =
CenterCrack_CFE_AdjustCrack . py ');
s1 = system ( cmd1 ) ;
file for new crack
s2 = system ( cmd2 ) ;
file for new crack
s3 = RunAbaqus ( jobName , CPUs , Memory ) ;
if s3 == 0
matName = grabDataFunc ( jobName ) ;
correctly , grab data
end
res = load ( matName ) ;
S. CFE . K1 (n) = mean ( res . K1 (3: end ));
contours 1 and 2
S. CFE . K2 (n) = mean ( res . K2 (3: end ));
S. CFE . JKs (n ) = mean ( res . JKs (3: end )) ;
S. CFE .J( n)
= mean ( res .J (3: end )) ;
% Adjust Python
% Adjust Abaqus
% Run Abaqus
% If executed
% Load Results
% Mean excluding
clear res
end
Al .( Rfld { m }) .( Sfld {m }) = S ;
structure
% Replace calling
end
AbaqusFileDelete ;
save ( ' Al_7075T6_Hudson_CFE . mat ', 'Al ')
123
Script: Abaqus_XFEM_terminate.m (Use with XFEM-PN)
% Terminates Abaqus when cyclic analysis increases by a single

cycle
% Must be run in a separate instance of Matlab while Abaqus is
running .
% Modified for Center Crack Plate , Hudson 1996 , 7075 - T6
%
% Rev
Date
Name
Changes
% 00
02 -17 -14
Josh Melson
- Created
%
format shortG
%% Define Variables and Strings
psTime = 10*60;
% Pause time in seconds
ckTime = 3*60;
cmdSta = ' abaqus python CenterCrack_XFEM_StatusCheck . py %s ';
cmdTer = ' abaqus job =% s terminate ';
%% Look for . lck file
endFile = [];
while isempty ( endFile ) == 1
lckFile = [];
while isempty ( lckFile ) == 1
lckFile = dir ( '*. lck ');
end
jobName = lckFile . name (: ,1: end -4) ;
fprintf ( ' Checking job = % s\n ', jobName )
staFile = [];
while isempty ( staFile ) == 1
staFile = dir ( strcat ( jobName , '. sta ')) ;
end
hltFile = [];
while isempty ( hltFile ) == 1
pause ( psTime )
s1 = system ( sprintf ( cmdSta , jobName ) );
hltFile = dir ( strcat ( jobName , '. hlt ')) ;
end
fprintf ( ' Terminating job = % s\n ', jobName )
s2 = system ( sprintf ( cmdTer , jobName ) );
pause ( ckTime )
endFile = dir ( ' JobComplete . end ');
124
if isempty ( endFile ) == 1
fprintf ( ' Continuing Abaqus Termination Process \ n\n ')
end
end
fprintf ( ' Program complete \n ')
delete ( ' JobComplete . end ');
125
Script: Abaqus_XFEM_MonteCarloAnalysis.m (Use with XFEM-PN)
% Creats new Input file , executes in Abaqus , and terminates at

single
% cycle steps . Modified for Center Crack Plate , Hudson 1996 , 7075 T6
%
% Rev
Date
Name
Changes
% 00
03 -11 -14
Josh Melson
- Created
%
format shortG
%% Load Data
dataPath = 'C :\ Users \ jmelson \ Documents \ Thesis \ ';
M = load ( strcat ( dataPath , ' Al_7075T6_Hudson_MatProp ')) ;
load ( strcat ( dataPath , ' Al_7075T6_Hudson '));
MC = load ( ' Al_7075T6_Hudson_MonteCarlo05base ');
%% Define Models to Analyze
Rfld = 'R2 ';
Sfld = ' Sm15 ';
Meth = ' Walker ';
S = Al .( Rfld ) .( Sfld );
W = M. Walker ;
dGc = (M. Kc ^2/ M .E) *(1 - S. R ^2) ;
%% Define CPU Properties
Memory = ' 18000 mb ';
CPUs
= '4 ';
%% Analyze in Abaqus
fprintf ( ' ****** RUN TERMINATION CODE IN PARALLEL ******\ n ')
fprintf ( ' Starting analysis of % d samples \ n ', length ( MC . C0 ))
cmd1Str = ' abaqus python CenterCrack_XFEM_MC_AdjustValues . py %.6 e
%.5 f %s ';
A
= struct ;
for m = 1: length ( MC . C0 )
[A .c3 , A. c4 ] = AbaqusGrowthConstants ( ...
Meth , S .R , MC . C0 ( m) , MC . m(m ) , M.E , MC . gam (m) );
MC . c3 (m ,1) = A. c3 ;
MC . c4 (m ,1) = A. c4 ;
if m < 10
runName = sprintf ( ' 00% d ' ,m);
else
runName = sprintf ( '0% d ',m );
126
% Define file name
end
cmd1 = sprintf ( cmd1Str , A. c3 , A.c4 , runName );
s1
= system ( cmd1 );
% Adjust Inp file for new
load
jobName = ( sprintf ( ' XFEM_R2_MC % s ' , runName )) ;
fprintf ( '\n \n ***** Starting Run %d (% s ) *****\ n ' , m , jobName )
s2 = RunAbaqus ( jobName , CPUs , Memory ) ; % Run Abaqus
matFile = grabODBdata_XFEM ( jobName );
res = load ( matFile ) ;
res . da = res .C (: ,1) ;
res . dG = res . G1 (: ,1) ;
res . c3 = A. c3 ;
res . c4 = A. c4 ;
res = rmfield ( res , { 'C '; ' G1 '} );
[ res .N , res . a] = LifeEstimate_XFEM_Walker ( res .da , res .dG , A.c3 ,
A. c4 , ' Poly ');
indx1 = find ( res . dG <= dGc , 1, ' last ') ;
% Reduce data to
dGc
if isempty ( indx1 ) == 0
MC . frm (m ,1) = res . frm ( indx1 );
indx2 = find ( res . a == res . da ( indx1 ));
if isempty ( indx2 ) == 1
indx2 = find ( res . a == res . da ( indx1 +1) );
end
MC . N(m ,1) = res . N( indx2 );
else
MC . frm (m ,1) = NaN ;
MC . N(m ,1)
= NaN ;
end
MC . mdl .( sprintf ( 'MC %s ', runName ) ) = res ;
clear res
jobKeep = dir ( '*. odb ');
AbaqusFileDelete ( jobKeep ,{ ' inp ' }) ;
end
save ( ' Al_7075T6_Hudson_MonteCarlo . mat ', '- struct ' , ' MC ')
save ( ' JobComplete . end ')
127
Function: LifeEstimate_XFEM_Walker.m (Use with XFEM-PN, Integrate Polynomial Method)
function [N , a , dN ] = LifeEstimate_XFEM_Walker ( ai , dG , c3 , c4 ,
Method )
% Estimate life based on integrated polynomial or trapezoidal
method
%
% Rev
Date
Name
Changes
% 00
02 -22 -14
Josh Melson
- Created
% 01
03 -14 -14
Josh Melson
- Add last increment w even
numbers
%
if nargin < 5;
Method = ' Poly ';
end
if max ( strcmp ( Method ,{ ' Poly '; ' Trapz '}) ) == 0

error ( ' Methods available are : Poly or Trapz ')
end
if rem ( length ( ai ) ,2) == 0
warning ( ' Vector length of a and dK must be even integers ')
aLast = 3;
else
aLast = 2;
end
if length ( ai ) ~= length ( dG )
error ( ' Vectors a and dK must be the same length ')
end
dadN_Abaq = @( delK ) c3 * delK .^ c4 ;
N = 0;
a = ai (1) ;
j = 2;
for n = 1:2: length ( ai ) - aLast
y = 1./ dadN_Abaq ( dG (n: n +2) );
x = ai ( n:n +2) ;
% da / dN function
% Cycles
% at crack length
% y vector
% x vector
switch Method
case ' Trapz '
dN (j -1) = trapz (x ,y) ;
% trapezoidal integration
case ' Poly '
pf = polyfit (x ,y ,2) ;
% polynomial fit
pfint = polyint ( pf );
% polynomial integration
dN (j -1) = polyval ( pfint ,x( end )) - polyval ( pfint ,x (1) );
end
128
end
N( j) = N(j -1) + dN (j -1) ;

a( j) = x( end ) ;
j = j + 1;
% Sum cycles
% Solve for last increment if even number comes in

if aLast == 3
y = 1./ dadN_Abaq ( dG ( end -1: end )) ;
x = ai ( end -1: end );
dN = trapz (x , y);
N( end +1) = N( end ) + dN ;
a( end +1) = x( end );
end
129
Appendix G
Python Code
The following code was developed in Python version 2.7.
CenterCrack_CFE_AdjustValues.py (Use with FEQP)
## Opens CenterCrack_CFE_AdjustCrack . py and replaces values for

crack analysis
#
# Rev
Date
Name
Changes
# 00
01 -15 -14
Josh Melson
- Created
#
# Modules
import os , sys
#- --- --- --- --- ---- --- --- --- -# Read in Variable names
Rvalue = sys . argv [1]
pLoad = sys . argv [2]
aLen
= sys . argv [3]
#- --- --- --- --- ---- --- --- --- -# Check for 3 input arguments
checkarg = len ( sys . argv ) - 1
numarg
= 3
if checkarg < numarg :
print " Fail "
sys . exit (" CenterCrack_CFE_AdjustCrack : Not enough input
arguments ")
elif checkarg > numarg :
print " Fail "
sys . exit (" CenterCrack_CFE_AdjustCrack : Too many input arguments
")
#- --- --- --- --- ---- --- --- --- -130
# Names to search for

rName = ' Rvalue = '
pName = ' pLoad = '
aName = ' aLen
= '
#- --- --- --- --- ---- --- --- --- -# Replace Names
rNameReplace = rName + '\ ' '+ Rvalue + '\ ' '
pNameReplace = pName + pLoad
aNameReplace = aName + aLen
#- --- --- --- --- ---- --- --- --- -# Read File and Replace
pyFile = open ( ' CenterCrack_CFE_AdjustCrack . py ') . read ()
newLine = '\ n '. join ( rNameReplace if line . startswith ( rName ) else
line for line in pyFile . splitlines () )
open ( ' CenterCrack_CFE_AdjustCrack . py ', 'w ') . write ( newLine )
newLine = '\ n '. join ( pNameReplace if line . startswith ( pName ) else
newLine = '\ n '. join ( aNameReplace if line . startswith ( aName ) else
131
CenterCrack_CFE_AdjustCrack.py (Use with FEQP)
# Adjust crack in CAE file

#
# Rev
Date
Name
# 00
01 -15 -14
Josh Melson
#
# Modules
from abaqus import *
from abaqusConstants import *
from caeModules import *
import mesh , sys
Changes
- Created
# Crack Lengths and Loading , Copied from Matlab

aLen
= 3.3386
pLoad = -6
Rvalue = 'R7 '
# Variable Names
fileName = ' CenterCrack_CFE '
mdlShort = ' CFE_short '
mdlLong = ' CFE_long '
partName = ' Plate '
instName = ' Plate -1 '
partionName = ' PartitionCrack '
loadName = ' StressLoad '
dimCrack = ' CrackLength '
dimRadius = ' EnrichRadius '
dimBuffer = ' BufferTop '
jobList
= []
# Update Abaqus Model
mdb = openMdb ( fileName )
if aLen <= 0.8:
mdlName = mdlShort
else :
mdlName = mdlLong
p = mdb . models [ mdlName ]. parts [ partName ]
s = p. features [ partionName ]. sketch
mdb . models [ mdlName ]. ConstrainedSketch ( name = ' __edit__ ', objectToCopy
=s )
s1 = mdb . models [ mdlName ]. sketches [ ' __edit__ ']
g , v , d , n = s1 . geometry , s1 . vertices , s1 . dimensions , s1 .
constraints
s1 . setPrimaryObject ( option = SUPERIMPOSE )
132
p. projectReferencesOntoSketch ( sketch =s1 , upToFeature = p. features [

partionName ] , filter = COPLANAR_EDGES )
s = mdb . models [ mdlName ]. sketches [ ' __edit__ ']
cBuffer = []
cRadius = []
if aLen <= 0.18:
cBuffer = 0.18
cRadius = 0.265625* aLen + 0.0534375
elif aLen > 0.18 and aLen <= 0.8:
cBuffer = aLen
cRadius = cBuffer *0.5625
else :
cBuffer = 0.8
cRadius = 0.45
s. parameters [ dimCrack ]. setValues ( expression = str ( aLen ))
s. parameters [ dimBuffer ]. setValues ( expression = str ( cBuffer ))
s. parameters [ dimRadius ]. setValues ( expression = str ( cRadius ))
s1 . unsetPrimaryObject ()
p. features [ partionName ]. setValues ( sketch = s1 )
del mdb . models [ mdlName ]. sketches [ ' __edit__ ']
p. regenerate ()
a = mdb . models [ mdlName ]. rootAssembly
a. regenerate ()
partInstances = ( a. instances [ instName ] , )
a. deleteMesh ( regions = partInstances )
e1 = a. instances [ instName ]. edges
seedRight = []
seedLeft = []
if aLen < 0.6:
seedRight = int ( (6 - aLen - cBuffer ) /5.72*13 )
seedCircm = 2
pickedEdges1 = e1 . getSequenceFromMask ( mask =( ' [#2800800 ] ',
), )
pickedEdges2 = e1 . getSequenceFromMask ( mask =( ' [#5 ] ', ) , )
pickedEdges3 = e1 . getSequenceFromMask ( mask =( ' [#2 c2a3550
#55995 ] ', ) , )
133
a. seedEdgeByBias ( biasMethod = SINGLE , end1Edges = pickedEdges1 ,

end2Edges = pickedEdges2 , ratio =4.0 , number = seedRight ,
constraint = FINER )
a. seedEdgeByNumber ( edges = pickedEdges3 , number = seedCircm ,
elif aLen >= 0.6 and aLen <= 0.8:
seedRight = int ( (6 - aLen - cBuffer ) /5.72*13 )
seedCircm = 3
pickedEdges1 = e1 . getSequenceFromMask ( mask =( ' [#2800800 ] ',
), )
pickedEdges3 = e1 . getSequenceFromMask ( mask =( ' [#2 c2a3550
#55995 ] ', ) , )
a. seedEdgeByNumber ( edges = pickedEdges3 , number = seedCircm ,
elif aLen > 0.8 and aLen <= 1.1:
seedLeft = 1
seedRight = int ( round ((6 - aLen - cBuffer ) /4.2*8 , 0 ) )
pickedEdges1 = e1 . getSequenceFromMask ( mask =( ' [#2000 ] ', ) ,
)
end2Edges = pickedEdges2 , ratio =4.0 , number = seedLeft ,
pickedEdges3 = e1 . getSequenceFromMask ( mask =( ' [#0 #2000 a0 ] '
, ), )
pickedEdges4 = e1 . getSequenceFromMask ( mask =( ' [#0 #60000 ] ',
), )
else :
seedRight = int ( round ((6 - aLen - cBuffer ) /4.2*8 , 0 ) )

if seedRight < 1:
seedRight = int ( 1 )
seedLeft = 10 - seedRight
pickedEdges1 = e1 . getSequenceFromMask ( mask =( ' [#2000 ] ', ) ,
)
134

end2Edges = pickedEdges2 , ratio =4.0 , number = seedLeft ,
pickedEdges3 = e1 . getSequenceFromMask ( mask =( ' [#0 #2000 a0 ] '
, ), )
pickedEdges4 = e1 . getSequenceFromMask ( mask =( ' [#0 #60000 ] ',
), )
a = mdb . models [ mdlName ]. rootAssembly
partInstances = ( a. instances [ instName ] , )
a. generateMesh ( regions = partInstances )
mdb . models [ mdlName ]. loads [ loadName ]. setValues ( magnitude = pLoad )
jobName = ' CFE_ '+ Rvalue + '_S '+ str ( int ( round ( pLoad ) )) [1:]+ '_C '+ str (
aLen ) [0]+ str ( aLen ) [2:5]
mdb . Job ( name = jobName , model = mdlName )
mdb . jobs [ jobName ]. writeInput ( consistencyChecking = OFF )
del mdb . jobs [ jobName ]
# jobList . append ( '%s ' % jobName )
135
CenterCrack_XFEM_MC_AdjustValues.py (Use with XFEM-PN Monte Carlo)
## Opens CenterCrack_CFE_AdjustCrack . py and replaces values for

crack analysis
#
# Rev
Date
Name
Changes
# 00
03 -14 -14
Josh Melson
- Created
#
# Modules
import os , sys , fileinput
#- --- --- --- --- ---- --- --- --- --- --- --- --- # Read in Variable names
c3 = sys . argv [1]
c4 = sys . argv [2]
num = sys . argv [3]
#- --- --- --- --- ---- --- --- --- --- --- --- --- # Check for 3 input arguments
checkarg = len ( sys . argv ) - 1
numarg
= 3
if checkarg < numarg :
print " Fail "
sys . exit (" CenterCrack_XFEM_MC_AdjustLoad : Not enough input
arguments ")
elif checkarg > numarg :
print " Fail "
sys . exit (" CenterCrack_XFEM_MC_AdjustLoad : Too many input
arguments ")
#- --- --- --- --- ---- --- --- --- --- --- --- --- # File Name
fileName = ' XFEM_R2_MC '
#- --- --- --- --- ---- --- --- --- --- --- --- --- # Names to search for
c3Name = 'c3 '
c4Name = 'c4 '
#- --- --- --- --- ---- --- --- --- --- --- --- --- # Replace Names
c3NameRepl = 'c3 = '+ c3
c4NameRepl = 'c4 = '+ c4
#- --- --- --- --- ---- --- --- --- --- --- --- --- # File Names
newFile = fileName + num
#- --- --- --- --- ---- --- --- --- --- --- --- --- # Read File , Replace , and Write
fIn = open ( fileName + '. inp '). read ()
fOut = open ( newFile + '. inp ', 'w ')
for line in fIn :
136
fOut . write ( line )

fOut . close ()
fLine = open ( newFile + '. inp '). read ()
fLine = '\n '. join ( c3NameRepl if line . startswith ( c3Name ) else line
for line in fLine . splitlines () )
fLine = '\n '. join ( c4NameRepl if line . startswith ( c4Name ) else line
for line in fLine . splitlines () )
open ( newFile + '. inp ', 'w '). write ( fLine )
137
CenterCrack_XFEM_StatusCheck.py (Use with all direct cyclic analyses
## Opens and examines status file for cycle increments

#
# Rev
Date
Name
Changes
# 00
02 -17 -14
Josh Melson
- Created
# 01
03 -11 -14
Josh Melson
- Add maxinc check
import numpy , fileinput , os , sys , fnmatch
## Read in values
fName = sys . argv [1]
## Define status file name
staName = fName + '. sta '
## Read data
staFile = fileinput . input ([ staName ])
## Matching String for Step 2, Any Increment Size
step2Str = '
2 *'
cycLast = 1
incAtOne = 0
maxInc
= 20
for line in staFile :
if fnmatch . fnmatch ( line , step2Str ):
lineNum = [ int ( r) for r in line . split () if r .
isdigit () ]
cycCurrent = lineNum [1]
if cycCurrent - cycLast == 1:
incAtOne = incAtOne + 1
cycLast = cycCurrent
if incAtOne >= maxInc :
open ( fName + '. hlt ' ,'w '). close () ;
else :
cycLast = cycCurrent
fileinput . close ()
138
msgextract_ConvergCriteria_2D.py (Use with all direct cyclic analyses)
## Opens and examines stabilization criteria for direct cyclic

#
# Rev
Date
Name
Changes
# 00
02 -10 -14
Josh Melson
- Created
#
import numpy , fileinput , os , sys , fnmatch , re
## Read in values
## Define msg file names
msgName = fName + '. msg '
## Read data
msgFile = fileinput . input ([ msgName ])
cycleStr = '* CYCLE * STARTS * '
iternStr = '* ITERATION * SUMMARY * '
avgFrcStr = '* AVERAGE FORCE * TIME AVG . FORCE * '
maxDispStr = '* MAX . COEFFICIENT OF DISP .* '
coefR0uStr = '* COEFF . OF RESI . FORCE ON CONST . TERM * '
coefRnuStr = '* COEFF . OF RESI . FORCE ON PERIODIC TERMS * '
corrU0uStr = '* CORR . TO COEFF . OF DISP .
ON CONST . TERM * '
corrUnuStr = '* CORR . TO COEFF . OF DISP .
ON PERIODIC TERMS * '
fourierStr = '* NUMBER OF FOURIER TERMS * '
allZeroStr = '* ALL * ZERO * '
directCycl = '*D I R E C T
C Y C L I C* '
staticStep = '*S T A T I C
A N A L Y S I S* '
cNumStr = []
iNumStr = []
avgFStr = []
avgMStr = []
maxUStr = []
maxRStr = []
R0uStr = []
RnuStr = []
R0rStr = []
RnrStr = []
U0uStr = []
UnuStr = []
U0rStr = []
UnrStr = []
frrStr = []
dirCycl = 0
for line in msgFile :
if fnmatch . fnmatch ( line , directCycl ):
dirCycl = 1
if fnmatch . fnmatch ( line , staticStep ):
dirCycl = 0
139
if dirCycl == 1:
# Cycle
if fnmatch . fnmatch ( line , cycleStr ):
cNum = [ int (s ) for s in line . split () if s. isdigit () ]
# Iteration
if fnmatch . fnmatch ( line , iternStr ):
iNum = [ int (r ) for r in line . split () if r. isdigit () ]
iNumStr . append ( str ( iNum [0]) )
cNumStr . append ( str ( cNum [0]) )
# Fourier
if fnmatch . fnmatch ( line , fourierStr ):
fNum = [ float (x) for x in re . findall ("[0 -9]+" , line )]
frrStr . append ( str ( fNum [0]) )
# Force
if fnmatch . fnmatch ( line , avgFrcStr ) :
avgF = [ float (x) for x in re . findall (" -?\ d +.?\ d *(?:[ Ee
] -?\+?\ d +) ?" , line )]
avgFStr . append ( str ( avgF [1]) )
if fnmatch . fnmatch ( line , maxDispStr ):
maxU = [ float (x) for x in re . findall (" -?\ d +.?\ d *(?:[ Ee
] -?\+?\ d +) ?" , line )]
maxUStr . append ( str ( maxU [0]) )
if fnmatch . fnmatch ( line , coefR0uStr ):
if fnmatch . fnmatch ( line , allZeroStr ):
R0u = [0]
else :
R0u = [ float (x) for x in re . findall (" -?\ d +.?\ d *(?:[
Ee ] -?\+?\ d +) ?" , line )]
R0uStr . append ( str ( R0u [0]) )
if fnmatch . fnmatch ( line , coefRnuStr ):
Rnu = [0]
else :
Rnu = [ float (x) for x in re . findall (" -?\ d +.?\ d *(?:[
Ee ] -?\+?\ d +) ?" , line )]
RnuStr . append ( str ( Rnu [0]) )
if fnmatch . fnmatch ( line , corrU0uStr ):
140
U0u = [0]
else :
U0u = [ float (x) for x in re . findall (" -?\ d +.?\ d *(?:[
Ee ] -?\+?\ d +) ?" , line )]
U0uStr . append ( str ( U0u [0]) )
if fnmatch . fnmatch ( line , corrUnuStr ):
Unu = [0]
else :
Unu = [ float (x) for x in re . findall (" -?\ d +.?\ d *(?:[
Ee ] -?\+?\ d +) ?" , line )]
UnuStr . append ( str ( Unu [0]) )
g = open ( fName + ' _Converge . txt ', 'w ')
for d in range ( len ( cNumStr )):
g. write ( '%s , %s , %s , %s , %s , %s , %s , %s , % s\n ' \
%( cNumStr [d ], iNumStr [d], frrStr [d] , \
avgFStr [ d], maxUStr [d ], R0uStr [ d], RnuStr [d] , U0uStr [d ], UnuStr
[d ]) )
g. close ()
141
msgextract_ConvergLimits_2D.py (Use with all direct cyclic analyses)
## Opens and examines stabilization limits for direct cyclic

#
# Rev
Date
Name
Changes
# 00
02 -10 -14
Josh Melson
- Created
#
import numpy , fileinput , os , sys , fnmatch , re
## Read in values
## Define msg file names
msgName = fName + '. msg '
## Read data
msgFile = fileinput . input ([ msgName ])
directCycl = '*D I R E C T
C Y C L I C* '
staticStep = '*S T A T I C
A N A L Y S I S* '
CR0 = '* CRITERION * RESI * CONST * '
CRn = '* CRITERION * RESI * ANY * '
CU0 = '* CRITERION * CORR * CONST * '
CUn = '* CRITERION * CORR * ANY * '
dirCycl = 0
for line in msgFile :
if fnmatch . fnmatch ( line , directCycl ):
dirCycl = 1
if fnmatch . fnmatch ( line , staticStep ):
dirCycl = 0
if dirCycl == 1:
# Limits
if fnmatch . fnmatch ( line , CR0 ) :
CR0val = [ float ( x) for x in re . findall (" -?\ d +.?\ d *(?:[
Ee ] -\ d +) ?" , line )]
CR0Str = str ( CR0val [0])
if fnmatch . fnmatch ( line , CRn ) :
CRnval = [ float ( x) for x in re . findall (" -?\ d +.?\ d *(?:[
Ee ] -\ d +) ?" , line )]
CRnStr = str ( CRnval [0])
if fnmatch . fnmatch ( line , CU0 ) :
CU0val = [ float ( x) for x in re . findall (" -?\ d +.?\ d *(?:[
Ee ] -\ d +) ?" , line )]
CU0Str = str ( CU0val [0])
142
if fnmatch . fnmatch ( line , CUn ) :

CUnval = [ float ( x) for x in re . findall (" -?\ d +.?\ d *(?:[
Ee ] -\ d +) ?" , line )]
CUnStr = str ( CUnval [0])
if all (( ' CR0val ' in locals () , ' CRnval ' in locals () , ' CU0val
' in locals () , ' CUnval ' in locals () )) :
break
g = open ( fName + ' _Limits . txt ' ,'w ')
g. write ( '%s , %s , %s , %s \n ' %( CR0Str , CRnStr , CU0Str , CUnStr ) )
g. close ()
143

Thesis Melson

Uploaded by

Copyright:

Available Formats

Thesis Melson

Uploaded by

Document Information

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Thesis Melson

Uploaded by

Copyright:

Available Formats

Fatigue Crack Growth Analysis

with Finite Element Methods

Thesis submitted to the Faculty of the

Robert L. West, Chair

Keywords: Fracture, eXtended Finite Element Method,

Copyright 2014, Joshua H. Melson

Fatigue Crack Growth Analysis

Low SIFs throughout growth resulted in life

and the large variability typically seen in growth

I would like to thank my adviser and committee chair, Dr.

Bob West, for his technical

guidance, availability to answer questions, and patience throughout my graduate career.

Norman Dowling, who's technical

John Kennedy, who helped guide me through

undergraduate and graduate school.

Overview and Statement of Need

Research Goals and Objectives . . . . . . . . . . . . . . . . . . . . . . . . . .

Mixed-Mode Critical Energy Release Rate

Fatigue Crack Growth Models . . . . . . . . . . . . . . . . . . . . . .

Onset of Crack Growth . . . . . . . . . . . . . . . . . . . . . . . . . .

Crack Extension Direction . . . . . . . . . . . . . . . . . . . . . . . .

Methods of Fracture Analysis

Finite Element Methods

Stochastic Analyses in Fracture Mechanics

Crack Growth Rates for 7075-T6 Aluminum Alloy

Fatigue Crack Growth Models . . . . . . . . . . . . . . . . . . . . . . . . . .

Paris Crack Growth Model . . . . . . . . . . . . . . . . . . . . . . . .

Walker Crack Growth Model . . . . . . . . . . . . . . . . . . . . . . .

Estimating Fatigue Fracture Life

Center Cracked Plate Model . . . . . . . . . . . . . . . . . . . . . . . . . . .

Forman Crack Growth Model

Finite Elements with Quarter Points (FEQP)

Virtual Crack Closure Technique (VCCT)

eXtended Finite Element Method with Phantom Nodes (XFEM-PN)

Monte Carlo Analysis

Random Numbers Selection

Conclusions and Future Work

Research Goals and Objectives . . . . . . . . . . . . . . . . . . . . . . . . . .

Appendix A Material Properties

Appendix B Center Cracked Plate Analysis Results

Appendix C Convergence Analysis for Finite Element Models

Appendix D Monte Carlo Analysis Results

Appendix E Compact Specimen Analysis Results

Appendix F MATLAB Code

Appendix G Python Code

Paris Equation Coecients . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Walker Equation Coecients . . . . . . . . . . . . . . . . . . . . . . . . . . .

Forman Equation Coecients

Analyzed Loading Conditions

Results from the Analytical Fatigue Fracture Analysis . . . . . . . . . . . . .

AFGROW Material Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

AFGROW Load Spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Results from the AFGROW Fatigue Fracture Analysis

Convergence analysis for the FEQP model

Results from the FEQP Fatigue Fracture Analysis . . . . . . . . . . . . . . .

Abaqus Crack Growth Constants. . . . . . . . . . . . . . . . . . . . . . . . .

(25.4 mm) in Fig. 4.7.

4.10 Stress Intensity Range Convergence Analysis.

4.11 Estimated Fatigue Fracture Life Convergence Analysis.

4.12 Results from the FEQP Fatigue Fracture Analysis . . . . . . . . . . . . . . .

Paris Equation Coecients . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Walker Equation Coecients . . . . . . . . . . . . . . . . . . . . . . . . . . .

Forman Equation Coecients

Schematic nite element representation of the VCCT method.

Paris equation t to Hudson's experimental data.

Walker equation t to Hudson's experimental data.

Forman equation t to Hudson's experimental data. . . . . . . . . . . . . . .

Simplied Center Crack Plate Model