Title: Computing Stress Intensity Factor Histories with the Finite Element
Method and Interfacing with AFGROW and NASGRO
Objective:
To illustrate the process of using a finite element
program, such as FRANC2D/L, to compute stress intensity factor histories
accurately, so that they can be transmitted to AFGROW or NASGRO for fatigue
crack growth rate (FCGR) and life predictions.
General
Description:
This problem details the processes of using the finite
element method to compute accurate relationships between stress intensity
factors and crack length, and of communicating these relationships to AFGROW or
NASGRO for further processing for FCGR and life predictions. Two problems are
used for illustration: a typical detail in the standard library of AFGROW or
NASGRO to confirm stress intensity factor accuracy, and a more complex problem
for which no standard library solution is available.
Topics Covered: Finite
element analysis, stress intensity factor computation, nonstandard geometry
and/or boundary conditions, geometry factors for AFGROW or NASGRO
Type of Structure: any
single or multiple layer planar structure amenable to 2D or axisymmetric
structural idealization
Relevant
Sections of Handbook: Sections 2, 5, 11
Author: Dr. A.
R. Ingraffea
Company Name: Fracture
Analysis Consultants, Inc.
121
Eastern Heights Drive
Ithaca,
NY 14850
6072574970
Contact Point: Dr. Paul Wawrzynek
Phone: 6072574970
eMail: wash@fracanalysis.com
Overview of Problem Description
An essential task in
predicting fatigue crack growth rate and remaining fatigue life is the accurate
computation of a relationship between crack length and stress intensity
factor. In Example 1 of this problem a
finite element code is used to compute such a relationship on a relatively
simple structural detail, Figure FAC1.1. The pin load, P, is distributed as a normal pressure following a cosine function
on the lower halfcircumference of the hole. To assess the accuracy of the
finite element calculations, they are compared to standard solutions stored in
NASGRO for this detail.
Figure FAC1.1.
Example 1 of this problem: a
structural detail from standard library in NASGRO (Problem TC03) and AFGROW
(Single Through Crack at HoleStandard Solution). Here W =10, D = 1, B = 5, t = 1, S_{o} = 10, and P = 10,
and c/D is varied from 0.09 to 1.35.
All in kips and inches. Figure taken from NASGRO 3.
Frequently, the structural detail of interest is not
contained in standard libraries. In
this case, either the actual detail must be idealized to resemble a standard
solution, or a finite or boundary element model must be constructed. The finite
element method can be used effectively in the role of general stress intensity
factor calculator. The method can
accommodate virtually any actual situation involving nonstandard geometry and
boundary conditions, and include multiple cracks, multiple materials, material
anisotropy, and initial stresses.
As an example of this situation, a variant of the detail
shown in Figure FAC1.1
is studied in Example 2 of this problem, Figure FAC1.2. This figure shows a simple lug
under nonsymmetric loading. A contactfit pin is inserted in the hole, and the
pin load, P, is distributed to the
hole by way of an elastic contact analysis.
A crack is initiated from a location of high tensile stress
concentration along the lug bore, and then allowed to propagate in mixed mode.
Figure
FAC1.2. Example 2 of this problem: a lug under unsymmetrical load. Lug and
pin are both steel, E= 29,000ksi, n=0.30, and frictionless, contactfit of the
pin in the lug is assumed. Thickness is
1.00 in.
Computational Models: Example 1
Two finite element models shown in Figures
FAC1.3 and FAC1.4 were used for Example 1 of this
problem. The purpose of using two
meshes is to show the degree of mesh refinement that is needed to achieve
accurate values of stress intensity factors using at least quadratic order
triangular and quadrilateral elements. There are about 2300 degreesoffreedom
in the initial, crackfree mesh shown in Figure FAC1.3, hereafter called the coarse
mesh.
There are about 9800 degreesoffreedom in the initial,
crackfree mesh shown in Figure FAC1.4, hereafter called the refined mesh. The number of degreesoffreedom in each
model increases as the mesh is updated to accommodate crack growth.
The mesh is refined
in the region of crack growth to the right of the hole in each case. A ruleofthumb in creating the initial mesh
in the region of expected crack growth is to make the average element size
there smaller than (say onehalf to onequarter) the specified increment in
crack length. Examples of this practice
are shown in Figures FAC1.3(b)
and FAC1.4(b) that show a detail of the
mesh around a crack tip. FRANC2D/L
automatically remeshes after each increment of crack growth. Figure FAC1.5
shows a detail of the refined mesh after a few steps of crack growth.
Figure FAC1.3. Example
1. (a) Coarse finite element model for structural detail shown in Figure
FAC1.1. (b) Detail of meshing near
initial crack tip, c/D=0.09.
Figures FAC1.3(b) and FAC1.4(b)
show some other important details of meshing around a crack tip. These include the relative size, type, and
distribution of elements immediately surrounding the crack tip. It is well known that mixedmode stress
intensity factors can be computed with excellent accuracy using the equivalent
domain formulation of the elastic JIntegral (Nikishkov 1987). FRANC2D/L can
automatically compute stress intensity factors with this technique. Other techniques, such as displacement
correlation (Chan 1970) and modified crack closure (Rybicki 1977), are also
available, but these have proven to be generally less accurate for a given
discretization. This formulation of the elastic JIntegral technique does not
require extreme levels of mesh refinement in the vicinity of a crack tip.
Figure
FAC1.4. Example 1. (a) Refined finite element model for structural detail
shown in Figure
FAC1.1. (b) Detail of meshing near initial crack
tip, c/D=0.09.
Experience has shown that the treatment shown in this
example, a rosette of eight sixnoded triangular elements immediately
surrounding the tip, quarterpoint versions of these elements, and element size
ranging from a few to as much as 25 percent of crack length, will produce very
accurate values of stress intensity factors.
The use of quarterpoint versions of these elements produces the
singular stress field known to exist in linear elastic fracture mechanics, is
automatically done in FRANC2D/L, and is well documented in the literature (Barsoum
1976).
Computational models like these require little effort to
generate. An experienced finite element
user can produce one in less than onehalf hour of persontime using codes such
as PATRAN^{tm} or ProE^{tm}.
Figure
FAC1.5. Example 1. Detail of
refined finite element model for structural detail shown in Figure FAC1.1
after a few steps of crack propagation, showing automatic remeshing.
Computational Results: Example 1
Table FAC1.1 shows a comparison between the stress intensity
factor histories from NASGRO and finite element analysis for the coarse
mesh. Each step of analysis took less
than one second of total computer time on a PC running Windows 2000 on a 300MHz
Pentium II processor. Shown are stress intensity factors due to the pin load, S_{3}, and the bypass load, S_{0}, and the total of these.
The average difference in total stress intensity factor across the range of
crack lengths analyzed was less than 1 percent.
Table FAC1.2 shows a comparison between the stress intensity
factor histories from NASGRO and finite element analysis for the refined
mesh. Each step of analysis took about
10 seconds of total computer time on a PC running Windows 2000 on a 300MHz
Pentium II processor. The average
difference in total stress intensity factor across the range of crack lengths
analyzed remained less than 1 percent, but did not improve over the results
from the coarse mesh. This implies that the results from finite element
analysis have converged to values acceptably close to those produced from the
boundary element method used to generate the values in NASGRO. It also implies
that relatively coarse meshes, such as those shown in Figure
FAC1.3, can be used to compute
stress intensity factors accurately when used with the equivalent domain
formulation of the elastic JIntegral. Finally, the personeffort and computer
time needed to produce an accurate stress intensity factor history for such 2D
details are not overwhelming.
Table FAC1.1. Example 1. Comparison
between NASGRO and finite element analysis stress intensity factor histories
for Example 1, using the coarse mesh shown in Figure FAC1.3.
Table FAC1.2. Example
1. Comparison between NASGRO and finite element analysis stress intensity
factor histories for Example 1, using the refined mesh shown in Figure FAC1.4.
Computational Models: Example 2
Figure FAC1.6 shows the initial finite element model for the lug
problem shown in Figure FAC1.2. The lug and its frictionless, contactfit pin are
both steel. The elastic contact problem between the pin and the lug is solved
in FRANC2D/L. Sixnoded, zerothickness
interface elements with nonlinear constitutive capability are inserted around
the contact surface. This process is
automated: the lug and pin are assigned different material property set numbers
(but the same properties), and a feature in FRANC2D/L finds the closed curve
defining the boundary between these two sets and automatically inserts
interface elements along the entire curve.
The constitutive model used to represent the normal
contact conditions is shown in Figure FAC1.7. No tension is allowed across the contact, and a high compressive
normal stiffness is assigned to minimize intrusion of the pin into the lug.
Figure FAC1.7. Example 2. Constitutive
model for normal stressdisplacement on the pin/lug contact surface.
A nonlinear solver is now needed, and FRANC2D/L uses a
dynamic relaxation procedure (Underwood 1983).
This is a relatively slow but extremely robust algorithm for solving
nonlinear equilibrium equations. Figure FAC1.8 presents the uncracked,
deformed shape and shows that separation of the pin from the lug has occurred
from about the 8 o’clock to the 2 o’clock positions, while nonoverlapping
contact occurs along the rest of the contact.
Figure FAC1.9
shows contours of major principal stress for the initial, uncracked
configuration. This solution involved
about 9200 DOF, and, with an error tolerance of 0.0005 on both equilibrium and
displacement change between time steps, required about 5200 time steps and 2.5
minutes on a PC running Windows 2000 on a 1GHz Pentium III processor.
Figure
FAC1.6. Example 2. Initial finite element mesh.
Figure FAC1.8. Example 2. Deformed shape
of uncracked configuration.
Amplification factor is 225.
Figure FAC1.9. Example 2. Contours of major principal stress in uncracked condition.
Computational Results: Example 2
Figure FAC1.8 indicates that there are two locations of high stress
concentration around the lug bore, as expected. A slightly higher concentration occurs in the lower left
quadrant, at about the 8 o’clock position, and a short crack (0.067 in.) is
initiated at the location of highest circumferential tensile stress near this
point. Growth of this crack is then simulated by:
1. Computing
K_{I} and K_{II} using the equivalent domain formulation of the
elastic JIntegral;
2. Calculating
the direction of crack growth using the maximum circumferential tensile stress
theory (Erdogan 1963);
3. Extending
the crack by 0.075 inch in this direction;
4. Resolving
the finite element problem, including the nonlinear contact between lug and
pin;
5. Repeating
steps 14 until the crack has extended about 2.5 inch at which point fatigue
life or residual strength limits are likely to have been reached.
Figure FAC1.10 shows the
corresponding amplified displaced shape, while Figure
FAC1.11 shows the resulting stress intensity factor histories. Figure FAC1.10 shows that the crack trajectory is
not quite radial, and is responding to the asymmetrical loading and geometry.
The trajectory is dictated here by the maximum circumferential tensile stress
theory that requires that K_{II }remain zero along the crack path. With a finite element model that discretizes
the trajectory into finite, straight segments, there will always be residual,
nonzero values of K_{II }computed at each crack tip location. If the
segments are short enough, these residuals should be small compared to the K_{I}
values. Figure
FAC1.12 shows that the values of K_{II }are oscillating_{ }around
zero, and are indeed small, in this case, never reaching more than 3.5% of K_{I. }The highest values usually occur early
in the trajectory while the finite element model is adjusting to the stress field
that is evolving as a result of crack growth, as shown in Figure FAC1.12.
The key practical issue suggested here is the length of crack growth
increment. This length should be
sufficiently short to accurately discretize a curvilinear trajectory and
provide enough data points for the accurate integration required for FCGR
calculations, while not being so short that excessive computation times
accrue. Here 32 increments were used. The number of DOF’s grew to nearly 14,000 at
the last increment, and a total of about 2 hours of computing time was
required.
Figure FAC1.10. Example 2. Final deformed shape of single cracked configuration. Amplification factor is 150.
As previously mentioned, a finite element code with
fracture mechanics features can be thought of as a general stress intensity
factor calculator. As such it can be
used to attack practically interesting variants of problems. For example, it is possible that two fatigue
cracks might initiate in this lug problem, one from each of the locations of
initially high stress concentration. This possibility is also simulated here,
under the assumptions that
Figure FAC1.11. Example 2. Computed stress
intensity factor histories. Single crack case.
Figure FAC1.12. Example 2. History of computed ratio of stress intensity factors. Single crack case.
the cracks initiate simultaneously, and that they have
equal rates of growth. The resulting
trajectories under these simplifying assumptions are shown in Figure FAC1.13. The corresponding mode I stress
intensity factor histories are given in Figure
FAC1.14. This figure shows that,
even if initiation were simultaneous, rate of growth would not be equal; the
left crack would have higher growth rate.
However, even under these simplifying assumptions, Figure FAC1.14 also shows that the growth rate of
the left crack would be higher than it would be if it were the only crack to
occur.
Figure FAC1.13. Example 2. Final deformed shape of multiply cracked configuration. Amplification factor is 100.
Figure FAC1.14. Example 2. Comparison of
stress intensity factor histories for single and multiple crack cases.
It is currently possible to couple the results of a stress intensity
factor history analysis from finite element analyses for a single crack to FCGR simulators like NASGRO and AFGROW. This
process is described next.
Interfacing with NASGRO and
AFGROW
This section describes the process of interfacing the
stress intensity factor history from finite element analyses with NASGRO
(NASGRO 2000) or AFGROW (AFGROW 2000) to do FCGR and remaining life
assessments.
Interfacing with NASGRO
In addition to its built in library of stress intensity
factor solutions for many standard structural geometries, NASGRO offers the
capability of performing crack growth analysis for nonstandard geometries
through its user defined data table (DT) approach. In this approach NASGRO requests stress
intensity factor input in the form of a onedimensional table (NASGRO:Reference
Manual, Appendix C):
a/D^{}


F3


0.067


0.52


0.142


0.96


0.217


1.34


0.292


1.66


0.367


1.97


0.442


2.23


Table FAC1.3. Portion of NASGRO stress
intensity factor data table for a nonNASGROstandard geometry. D= 4.0 in. and
S_{3 }= 7.5 ksi.
NASGRO’s general expression for stress intensity factors
is,
where the stress quantities S_{0}, S_{1} and S_{2}, S_{3}, and S_{4} are for applied tension/compression, bending in the
thickness and width directions, pin bearing pressures, and special cases,
respectively. For Example 2 only S_{3}
is nonzero, so NASGRO expects a relationship of the form
_{}_{}
where F_{3}_{
}is the geometry correction factor_{ }for the present lug
problem. The values of F_{3}_{ }needed for the
NASGRO data table are contained within the stress intensity factor values
computed with FRANC2D/L, and can be explicitly obtained through
_{}
Table FAC1.3 contains a
portion of the geometry correction factor history obtained from FRANC2D/L for
Example 2 with a single crack. Herein D is the pin diameter, 4.0 in., and S_{3} is the pin bearing stress
S_{3} = P/Dt = 30 kips/((4
in.)(1.0 in)) = 7.5 ksi
Once stress intensity factor history data in this form has
been supplied, all standard NASGRO fatigue crack growth capabilities become
available for any cracked structural model for which FRANC2D/L is an
appropriate modeling system.
Interfacing with AFGROW
In addition to its built in library of stress intensity
factor solutions for many standard structural geometries, AFGROW offers the
capability of performing crack growth analysis for nonstandard geometries
through its userdefined Beta Table approach (AFGROW: Users Guide). In this approach for a single throughcrack,
AFGROW requests stress intensity factor input in the form of a onedimensional
table of crack length versus stress intensity factor geometry correction, or
beta, values, where beta is defined by
_{}
and where s is the relevant applied stress
and c is crack length. The beta values are obtained from finite
element results in the same way as the equivalent F_{i} values are obtained for NASGRO. Figure FAC1.15 shows the UserInput Beta Table
dialogue box in AFGROW with the first few beta values obtained for Example 2
with a single crack, and with s set to unity.
Once stress intensity factor history data in this form has
been supplied, all standard AFGROW fatigue crack growth capabilities become
available for any cracked structural model for which FRANC2D/L is an
appropriate modeling system.
Figure FAC1.15. AFGROW dialogue box for
creating userinput beta table, showing portion of table for Example 2 with a
single crack. Here s has been set to unity.
References
AFGROW: Users Guide and Technical Manual, AFGROW for
Win95/98/NT4/2K, Version 4.0001011.8, AFRLVAWPTR2000XXXX, May 2000.
Barsoum, R.S. (1976), "On the use of Isoparametric
Finite Elements in Linear Fracture Mechanics," International Journal for Numerical Methods in Engineering, Volume
10 (1976) pp 2538.
Chan, S.K, Tuba, I.S., and Wilson, W.K., "On the
Finite Element Method in Linear Fracture Mechanics," Engineering Fracture Mechanics, Volume 2 (1970) pp 117.
Erdogan, F., and Sih, G.C., "On the Crack Extension
in Plates Under Plane Loading and Transverse Shear," Journal of Basic Engineering, Volume 85 (1963) pp 519525.
NASGRO: Reference Manual, Fatigue
Crack Growth Computer Program “NASGRO” Version 3.0, JSC22267B, March 2000.
Nikishkov, G.P. and Atluri, S.N. (1987), "Calculation
of Fracture Mechanics Parameters from an Arbitrary ThreeDimensional Crack by
the Equivalent Domain Integration Method," International Journal for Numerical Methods in Engineering, Volume
24 (1987) pp 18011821.
Rybicki, E.R., and Kanninen, M. (1977), "A Finite
Element Calculation of Stress Intensity Factors by a Modified Crack Closure
Integral," AIAA Journal Volume
2, 1977.
Underwood, P., (1983)"Dynamic Relaxation", in Computational Methods for Transient
Analysis, Volume 1. Belytchko T., and Hughes, J.R. Eds. NorthHolland,
Amsterdam, 1983.