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, non-standard 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
607-257-4970
Contact Point: Dr. Paul Wawrzynek
Phone: 607-257-4970
e-Mail: 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 FAC-1.1. The pin load, P, is distributed as a normal pressure following a cosine function
on the lower half-circumference 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 FAC-1.1.
Example 1 of this problem: a
structural detail from standard library in NASGRO (Problem TC03) and AFGROW
(Single Through Crack at Hole-Standard Solution). Here W =10, D = 1, B = 5, t = 1, So = 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 non-standard 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 FAC-1.1
is studied in Example 2 of this problem, Figure FAC-1.2. This figure shows a simple lug
under non-symmetric loading. A contact-fit 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
FAC-1.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, contact-fit of the
pin in the lug is assumed. Thickness is
1.00 in.
Computational Models: Example 1
Two finite element models shown in Figures
FAC-1.3 and FAC-1.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 degrees-of-freedom
in the initial, crack-free mesh shown in Figure FAC-1.3, hereafter called the coarse
mesh.
There are about 9800 degrees-of-freedom in the initial,
crack-free mesh shown in Figure FAC-1.4, hereafter called the refined mesh. The number of degrees-of-freedom 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 rule-of-thumb in creating the initial mesh
in the region of expected crack growth is to make the average element size
there smaller than (say one-half to one-quarter) the specified increment in
crack length. Examples of this practice
are shown in Figures FAC-1.3(b)
and FAC-1.4(b) that show a detail of the
mesh around a crack tip. FRANC2D/L
automatically remeshes after each increment of crack growth. Figure FAC-1.5
shows a detail of the refined mesh after a few steps of crack growth.
Figure FAC-1.3. Example
1. (a) Coarse finite element model for structural detail shown in Figure
FAC-1.1. (b) Detail of meshing near
initial crack tip, c/D=0.09.
Figures FAC-1.3(b) and FAC-1.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 mixed-mode stress
intensity factors can be computed with excellent accuracy using the equivalent
domain formulation of the elastic J-Integral (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 J-Integral technique does not
require extreme levels of mesh refinement in the vicinity of a crack tip.
Figure
FAC-1.4. Example 1. (a) Refined finite element model for structural detail
shown in Figure
FAC-1.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 six-noded triangular elements immediately
surrounding the tip, quarter-point 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 quarter-point 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 one-half hour of person-time using codes such
as PATRANtm or ProEtm.
Figure
FAC-1.5. Example 1. Detail of
refined finite element model for structural detail shown in Figure FAC-1.1
after a few steps of crack propagation, showing automatic remeshing.
Computational Results: Example 1
Table FAC-1.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, S3, and the bypass load, S0, 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 FAC-1.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
FAC-1.3, can be used to compute
stress intensity factors accurately when used with the equivalent domain
formulation of the elastic J-Integral. Finally, the person-effort and computer
time needed to produce an accurate stress intensity factor history for such 2D
details are not overwhelming.
Table FAC-1.1. Example 1. Comparison
between NASGRO and finite element analysis stress intensity factor histories
for Example 1, using the coarse mesh shown in Figure FAC-1.3.
Table FAC-1.2. Example
1. Comparison between NASGRO and finite element analysis stress intensity
factor histories for Example 1, using the refined mesh shown in Figure FAC-1.4.
Computational Models: Example 2
Figure FAC-1.6 shows the initial finite element model for the lug
problem shown in Figure FAC-1.2. The lug and its frictionless, contact-fit pin are
both steel. The elastic contact problem between the pin and the lug is solved
in FRANC2D/L. Six-noded, zero-thickness
interface elements with non-linear 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 FAC-1.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 FAC-1.7. Example 2. Constitutive
model for normal stress-displacement on the pin/lug contact surface.
A non-linear 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 FAC-1.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 non-overlapping
contact occurs along the rest of the contact.
Figure FAC-1.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
FAC-1.6. Example 2. Initial finite element mesh.
Figure FAC-1.8. Example 2. Deformed shape
of uncracked configuration.
Amplification factor is 225.
Figure FAC-1.9. Example 2. Contours of major principal stress in uncracked condition.
Computational Results: Example 2
Figure FAC-1.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
KI and KII using the equivalent domain formulation of the
elastic J-Integral;
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 non-linear contact between lug and
pin;
5. Repeating
steps 1-4 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 FAC-1.10 shows the
corresponding amplified displaced shape, while Figure
FAC-1.11 shows the resulting stress intensity factor histories. Figure FAC-1.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 KII remain zero along the crack path. With a finite element model that discretizes
the trajectory into finite, straight segments, there will always be residual,
non-zero values of KII computed at each crack tip location. If the
segments are short enough, these residuals should be small compared to the KI
values. Figure
FAC-1.12 shows that the values of KII are oscillating around
zero, and are indeed small, in this case, never reaching more than 3.5% of KI. 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 FAC-1.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 FAC-1.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 FAC-1.11. Example 2. Computed stress
intensity factor histories. Single crack case.
Figure FAC-1.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 FAC-1.13. The corresponding mode I stress
intensity factor histories are given in Figure
FAC-1.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 FAC-1.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 FAC-1.13. Example 2. Final deformed shape of multiply cracked configuration. Amplification factor is 100.
Figure FAC-1.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 non-standard geometries
through its user defined data table (DT) approach. In this approach NASGRO requests stress
intensity factor input in the form of a one-dimensional 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 FAC-1.3. Portion of NASGRO stress
intensity factor data table for a non-NASGRO-standard geometry. D= 4.0 in. and
S3 = 7.5 ksi.
NASGRO’s general expression for stress intensity factors
is,
|
|
NASGRO Eqtn. (2.2)
|
where the stress quantities S0, S1 and S2, S3, and S4 are for applied tension/compression, bending in the
thickness and width directions, pin bearing pressures, and special cases,
respectively. For Example 2 only S3
is non-zero, so NASGRO expects a relationship of the form
where F3
is the geometry correction factor for the present lug
problem. The values of F3 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 FAC-1.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 S3 is the pin bearing stress
S3 = 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 non-standard geometries
through its user-defined Beta Table approach (AFGROW: Users Guide). In this approach for a single through-crack,
AFGROW requests stress intensity factor input in the form of a one-dimensional
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 Fi values are obtained for NASGRO. Figure FAC-1.15 shows the User-Input 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 FAC-1.15. AFGROW dialogue box for
creating user-input 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, AFRL-VA-WP-TR-2000-XXXX, 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 25-38.
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 1-17.
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 519-525.
NASGRO: Reference Manual, Fatigue
Crack Growth Computer Program “NASGRO” Version 3.0, JSC-22267B, March 2000.
Nikishkov, G.P. and Atluri, S.N. (1987), "Calculation
of Fracture Mechanics Parameters from an Arbitrary Three-Dimensional Crack by
the Equivalent Domain Integration Method," International Journal for Numerical Methods in Engineering, Volume
24 (1987) pp 1801-1821.
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. North-Holland,
Amsterdam, 1983.