
Citation 
 Permanent Link:
 https://ufdc.ufl.edu/UFE0015560/00001
Material Information
 Title:
 Anisotropic Fracture Analysis of the BX265 Foam Insulation Material Under MixedMode Loading
 Creator:
 KNUDSEN, ERIK ( Author, Primary )
 Copyright Date:
 2008
Subjects
 Subjects / Keywords:
 Coordinate systems ( jstor )
Engineering ( jstor ) Foams ( jstor ) Fracture mechanics ( jstor ) Fracture strength ( jstor ) Geometric angles ( jstor ) Geometric lines ( jstor ) Insulation ( jstor ) Room temperature ( jstor ) Specimens ( jstor )
Record Information
 Source Institution:
 University of Florida
 Holding Location:
 University of Florida
 Rights Management:
 Copyright Erik Knudsen. Permission granted to University of Florida to digitize and display this item for nonprofit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
 Embargo Date:
 8/31/2006
 Resource Identifier:
 649815428 ( OCLC )

Downloads 
This item has the following downloads:

Full Text 
ANISOTROPIC FRACTURE ANALYSIS OF THE BX265 FOAM INSULATION
MATERIAL UNDER MIXEDMODE LOADING
By
ERIK KNUDSEN
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
2006
Copyright 2006
by
Erik Knudsen
This thesis is dedicated to my parents, friends, family, lab mates and everyone else who
was there for me when I needed them the most. I will never forget you!
ACKNOWLEDGMENTS
First and foremost I would like to thank my parents for their unwavering love and
support. My advisor, Dr. Nagaraj Arakere, has been in my corner since I first came to
graduate school and his contacts at NASA paved the way for the two degrees I now hold
in engineering. In the summer of 2002 I accompanied Dr. Arakere to Huntsville as part
of the Summer Faculty Fellowship Program. During our ten week stay, I met many
experts and scientists who work at Marshall, particularly Dr. Greg Swanson, Jeff
Raybum, Greg Duke, Dr. Gilda Battista, Stan Oliver, Preston McGill, Doug Wells, Dr.
Eric Poole, Dr. Philip Allen and Denise Duncan (in a later trip), and many, many others.
They have helped me grow as a student and a person, and I am forever indebted.
After my trip to NASA, approximately one year later, I earned a NASA Fellowship
to sponsor this work. Ed Adams of the Education Directorate at Marshall deserves
special praise here, and I thank him for blessing me with the opportunity to, once more,
work with NASA. I shared an on office Dr. Philip Allen and one could not ask for a
more pleasant person to work with and learn from.
The Comell Fracture Group was also a tremendous help in my quest to better
understand and learn the principles of fracture mechanics. Drs. Paul "Wash" Wawryznek
and Bruce Carter are two of the nicest and most patient people I have ever known and I
owe them many thanks for taking the time to help me.
Last, but certainly not least, I wish to thank my labmates over the years: Srikant,
Sangeet, Shadab, Niraj, Shannon, Yogen, Jeff, TJ, and George. Without their insight,
jokes, and patience I would not be where I am today. I now consider many of them to be
my best friends and I hope we can stay in touch.
TABLE OF CONTENTS
A C K N O W L E D G M E N T S ................................................................................................. iv
LIST OF TABLES .... ...... .. ........................... ........... .................... .. viii
LIST OF FIGURES ......... ......................... ...... ........ ............ ix
ABSTRACT .............. ..................... .......... .............. xii
CHAPTER
1 BA CK GR OU N D ......... ....................... .... ........ ..... ............................ .
Problem M motivation ......... .. ...... ...... ......... .. ....... .......... ...............
Program Scope and Focus ................................................. ............................... 2
2 FOAM TESTING AND MATERIAL CLASSIFICATION ...............................6
Cellular Solid B background Inform action ............................................ .....................6
Initial F oam M material T testing ................................................................ ....................7
Determination of Elastic Constants and Material Classification..........................7
Fracture Testing for the BX265 Foam Material......................................20
Additional Testing: Divot Test Specim ens.............................................................. 24
Sum m ary ...................................... ................. ................. .......... 28
3 AN INTRODUCTION TO FRACTURE MECHANICS..........................................30
M material D efin ition s ......................................................................... ..................... 3 0
Isotropic Fracture M echanics ...................... ....................... ................... 31
Linear Elastic Fracture Mechanics for Anisotropic Materials............................... 33
Determining the Stress Intensity FactorK............... ..... ............... 37
FRANC3D Next Generation Software ............................................ ............... 38
Computing SIFs Using Displacement Correlation............................................45
Effects of Temperature on the SIF Solution...........................................................52
Sum m ary ...................................... ................. ................. .......... 54
4 CRACK TURNING THEORY AND FINITE ELEMENT MODELING ...............56
C rack T u rning T h eory ...................................................................... ................ .. 56
FE M odeling of the M (T) Specim en ................................. ...................................... 64
5 RESULTS AND DISCU SSION ........................................... .......................... 68
Effect of Material Orientation on the SIF Solution ............................... ...............68
Crack Turning Predictions ....... ................ .... ..... ... ....... ...................76
M(T) Specimen Fabrication and Determination of Local Knit Line Axes .........77
Local Crack Tip Com putations ........................................ ....................... 78
Influence of Thermal Loads on the SIF Solution ....................................... .......... 85
Sum m ary ............... .... ..... ......... ......................................93
6 CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE WORK ..............95
C o n c lu sio n s............................... .......................................................................... 9 5
Recom m endations for Future W ork ........................................ ....................... 96
APPENDIX
GENERAL SOLUTION FOR THE PLANE PROBLEM OF A LINEAR ELASTIC
A N ISO TR O PIC B O D Y ................................................................... .....................98
L IST O F R E FE R E N C E S ........................................................................ ................... 103
BIOGRAPHICAL SKETCH ............................................................. ............... 108
LIST OF TABLES
Table page
21 Measured material property values from experiment ........................ ...........18
31 Various m ethods for determ ining the SIF ..................................... .................38
41 Direction cosines .......... ....... ........ .................... ................... 66
51 Direction cosines for the case one material orientation .......................................73
52 KI for cases 08, 0 degree inclination ........................................... ............... 76
53 KI for cases 916, 0 degree crack inclination ................ .................................. 76
54 Measured and predicted turning angles (in degrees)................... ... .............83
55 Sum m ary of m idplane K IK II......... ............... ........................... .. ...................91
56 Sum m ary of m idplane K ill......... ............... .................................. ............... 91
LIST OF FIGURES
Figure page
11 Foam loss during Shuttle ascension ........................................ ........ ............... 2
12 Sprayon foam insulation (SO FI)................................................................... ......3
21 Idealized foam microstructure...................................................... ........................ 8
22 Square specimen tensile test used to determine the Poisson ratio..........................15
23 F oam test specim ens........... ............................................................ ...... .... ... ... 15
24 Summary of stressstrain data .............................................. 16
25 Coordinate system for the transversely isotropic material............................... 17
26 Young's moduli vs. temperature .................... ........................... 18
27 Shear m oduli vs. tem perature...................... ................................. ............... 19
28 Coefficient of thermal expansion vs. temperature ...........................................19
29 Thermal conductivity vs. temperature....... .......... ................... ..............20
210 From left to right: M(T), C(T), and SNE(B) fracture test specimens ....................22
211 Summary Kvalues for various fracture specimens ...........................................23
212 Clip gauge on C(T) specimen used to measure CMOD .....................................24
213 Typical load vs CMOD curve for a brittle material ...........................................25
2 14 D iv ot test p an el............ .... ........................................................................... .. ...... .. 2 6
215 2D view of divot test................................................................. ............... 27
216 2D view of test panel after heat flux has been applied .......................................27
217 Foam blow out from divot test ............................... ............... 28
31 Coordinate system used in equation (31)................................... ............... 31
32 From left to right: mode I (opening), mode II (shearing), mode III (tearing) .........32
33 Path of the Jintegral .............................................. .... .... .. ........ .... 39
35 D definition of Irwin's crack closure integral ............................ .............................43
36 Correlation points typically used for displacement correlation ............................46
37 Triangular quarter point elem ents ........................................ ........................ 48
38 Isoparametric element and degenerated element with midside nodes moved to
quarter point locations .............................................................................48
41 D definition of crack turning angle ........................................ ......................... 57
42 Elliptical relationship for T ............................................................................58
43 D definition of hoop stress ................................... .... ............................ ...............59
44 Orthotropic toughness values ........................................................ ............. 62
45 Definition of the a and n vectors ................................................... ............... 62
47 M (T) m odel created in AN SY S ........................................ .......................... 65
48 Sam ple transform action ................................................. ................................ 66
49 Fractured M(T) specimen: the dotted line is added to show the location of the
k n it lin e ............................................................................ 6 7
51 Definition of crack inclination angle.................................................................. 69
52 D definition of crack front distance....................................... .......................... 69
53 M (T) m odel built in A N SY S ........................................................................ .. .... 70
54 D definition of m material orientation..........................................................................70
55 Top view of the cones shown in Figure 54......................................................71
56 Hypothetical material orientation relative to the ET..............................................71
57 Definition of case one material orientation for the M(T) FE model ......................72
58 Mode I SIF for cases 14; 0 degree crack inclination ...........................................74
59 Mode I SIF for cases 14; 30 degree crack inclination..................... ...............74
510 Mode II SIF for cases 14; 30 degree crack inclination .......................................75
511 Mode III SIF for cases 14; 30 degree crack inclination.............................75
512 From left to right: front, left, right, and rear sides of the M(T) specimen..............77
513 D eterm nation of knit line plane..............................................................................78
514 M (T) m odel with inclined crack ..................................................... ............. 79
515 Cartesian and cylindrical stresses................................ ................... ...... ........ 81
516 Crack turning angle .................. ............................... ..... .. ............ 82
517 Application of therm al loads .............................................................................87
518 Room temperature KIKIII (F3D) vs. KIKIII (DC) with T1 = 0F and T2 =
10 0 F ............................................................................... 8 8
519 Room temperature KI (F3D) vs. KI (DC) with T1 = 2000F and T2 = 100F..........88
520 (p = 30 degrees: room temperature KI (F3D) vs. KI (DC) with T1 = 0F and T2 =
10 0 F ............................................................................... 8 9
521 (p = 30 degrees: room temperature KI (F3D) vs. KI (DC) with T1 = 2000F and
T 2 = 10 0 F ........................................................................ 9 0
522 Temperature distribution along the crack front..................................................... 92
523 KIKIII vs. normalized crack front distance; TG2 applied in thickness direction...92
Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy
ANISOTROPIC FRACTURE ANALYSIS OF THE BX265 FOAM INSULATION
MATERIAL UNDER MIXEDMODE LOADING
By
Erik Knudsen
August 2006
Chair: Nagaraj K. Arakere
Major Department: Mechanical and Aerospace Engineering
The National Aeronautics and Space Administration (NASA) uses a closedcell
polyurethane foam to insulate the external tank (ET) which contains the liquid oxygen
and hydrogen for the Space Shuttle main engines. This is a type of sprayon foam
insulation (SOFI), similar to the material used to insulate attics in residential
construction. In February of 2003, the Shuttle Columbia suffered a catastrophic accident
during reentry. Debris from the ET impacting the Shuttle's thermal protection tiles
during liftoff is believed to have caused the Space Shuttle Columbia failure during re
entry.
NASA engineers are very interested in understanding the processes that govern the
breakup/fracture of this complex material from the shuttle ET. The foam is anisotropic in
nature and the required stress and fracture mechanics analysis must include the effects of
the direction dependence on material properties. Over smooth, flat areas of the ET the
foam can be sprayed down in a very uniform fashion. However, near bolts and fitting
points it is possible for voids and other defects to be present after the foam is applied.
Also, the orientation of the foam, as it rises from nonuniform surfaces, can be
significantly different from the orientation over typical acreage sprays. NASA believes
that air present in these small voids is liquefied and then turned into a gas during liftoff.
While the Shuttle is ascending to space, the pressure in these cavities can become large
enough to force a subsurface crack toward the exterior of the tank, thus freeing portions
of foam insulation.
As a first step toward understanding the fracture mechanics of this complex
material, a general theoretical and numerical framework is presented for computing stress
intensity factors (SIFs), under mixedmode loading conditions, taking into account the
material anisotropy. The effects of material orientation and mode mixity on the
anisotropic SIF solution are analyzed. Crack turning predictions under mixed mode
loading are presented. Furthermore, the influence of temperature gradients on the SIF
solution is studied, in view of the thermal gradient present through the foam thickness on
the ET. The results presented represent a quantitative basis for evaluating the strength
and fracture properties of anisotropic BX265 foam insulation material.
CHAPTER 1
BACKGROUND
Problem Motivation
On February 1st, 2003, the Space Shuttle Columbia suffered a catastrophic failure
during reentry. NASA has conducted an exhaustive investigation of the failure and the
consensus now is the breakup was caused by a segment of BX265 foam insulation,
roughly the size of a suitcase, striking the wing during liftoff. The foam impact, in prior
Shuttle launches, has been known to cause impact damage to the thermal protection tiles,
but was not considered to be a serious problem until the Columbia disaster. The damage
to the thermal protection tiles on the leading edge of the wing is thought to have allowed
gases to penetrate the wing during reentry, triggering a cascading series of catastrophic
events that led to the loss of the Shuttle.
It is believed that when the foam insulation is applied to the external tank (ET) that
voids are created in certain areas where geometric discontinuities such as bolts, flanges,
and fittings are encountered. Since the tank is filled with liquid oxygen and hydrogen
the foam is exposed to cryogenic temperatures at the ET's surface. The air inside these
voids can be liquefied and nitrogen can condense into these small cavities. During liftoff
the outer surface of the foam is exposed to aerodynamic heating. This heating raises the
temperature of the liquid nitrogen, turning it into a gas. The pressure difference
associated with the formation of the gas can cause pieces of foam to be blown out during
liftoff( see Figure 11).
Figurel1. Foam loss during Shuttle ascension
Program Scope and Focus
Prior to the Columbia disaster, the foam that insulates the tank was not understood
at a fundamental level. Quantities like the elastic moduli, and thermal properties such as
coefficient of thermal expansion are needed for a meaningful failure analysis of this
material. This prompted NASA to initiate an intensive testing program to ascertain
various properties at cryogenic, room, and elevated temperatures. The pictures in Figure
11 show the foam breaking off areas near fitting points. At these locations, insulation is
sprayed around various parts, such as bolts, and it is here that the foam's orientation,
relative to the tank, can change by quite a bit. For typical acreage sprays, there is very
little difference between the local coordinate axes of the foam and those of the substrate
(Figure 12).
Pictures of the Shuttle during liftoff, such as those shown in Figure 11, show us
that large pieces of foam are not really breaking away over smooth, open, areas of the
tank. The problem seems to be confined to areas where the foam is applied over fitting
points, bolts, and other areas where various parts can protrude the surface of the ET.
Spray
Direction
" 11
Material
coordinates
Knitines 
Knitlines Note: For typical acreage sprays,
SubstrateS hoop
Substrate H
Coordinates A=axial
R=radial
If theta
and sub state coordinates are aligned
Figure 12. Sprayon foam insulation (SOFI)
Right away we notice that two potential problems come to mind at these locations: it is
possible for voids to be created when the foam is being laid down, and the orientation of
the foam can change by a large degree.
The material tests conducted by NASA have shown that this type of foam is brittle.
Perhaps the best way to analyze sharp, cracklike, defects for this sort of material is
through linear elastic fracture mechanics (LEFM). For isotropic materials, the neartip
stresses are completely characterized by Kthe stress intensity factor (SIF). The
analytical Ksolutions for anisotropic materials take a similar form, but have additional
terms that come from the various elastic moduli. The formulations of neartip stresses for
isotropic materials are well established and are now incorporated in many finite element
(FE) software packages. Some Ksolutions for certain kinds of directiondependent
materials started appearing in the 1960s and more generalized solutions were presented in
the 1980s. Since many engineering materials are not isotropic (such as rolled aluminum
in aircraft fuselages and singlecrystal superalloys used in high performance gas turbine
engines), more attention is being paid toward the implementation of anisotropic K
solutions in various FE schemes. Still, many of the most popular FE software packages
do not allow accurate computations of SIFs for simulations involving anisotropic
materials.
Foams are inherently anisotropic, and NASA, after extensive testing, has classified
this material as transversely isotropic. Such materials also have fracture characteristics
that are direction dependent and both of these traits must be accounted for in the
analytical and numerical analyses. Although anisotropic Ksolutions are not readily
available, the Cornell Fracture Group has developed a software package, FRANC3D, that
simplifies the meshing process necessary for FE modeling of cracklike defects in elastic
solids and also computes a Ksolution for anisotropic materials. Concurrently, an
alternative method for computing the SIFs, based on the work of Hoenig (1982), for
anisotropic materials has been employed and extended to account for thermal stresses and
strains.
It is hoped that the present study will enhance NASA's understanding of this
material, since the Ksolutions used in this analysis account for both the direction
dependent properties and the impact of thermal stresses and strains. The following points
summarize the scope of this study in a broad sense:
1. Through the use of the middle tension, M(T), specimen and the FE method a
systematic rotation of material properties and its impact on the Ksolution is
examined
2. The impact of mode mixity on the Ksolution is studied
3. Crack turning predictions that take into account both the anisotropic stiffness and
fracture properties are presented
5
4. The effect of a thermal gradient through the M(T) specimen and its impact on the
Ksolution are presented
CHAPTER 2
FOAM TESTING AND MATERIAL CLASSIFICATION
Cellular Solid Background Information
Cellular solids are generally comprised into two sets of structural categories: 1)
honeycomb layouts where the microstructure is comprised of twodimensional arrays of
polygons that fill a planar area and 2) foams, the threedimensional counterpart to
honeycombs where the cells take the form of polyhedra. Within Gibson and Ashby's
(1988) definitive text on cellular solids, we see how many types of materials can be
foamed: metals, glasses, and ceramics can be fabricated into cells that take on a variety of
shapes and sizes. Among the most common types of foam are polymers which are used
for many applications including insulation for residential homes or disposable coffee
cups. Polymer foams are usually very light and they are popular in industries where
weight is of prime importance, such as in the aerospace industry.
Foams are typically made by bubbling a gas into the liquid monomer or hot
polymer. The gas can be introduced mechanically by stirring or by using a blowing
agent. Physical blowing agents are normally gases like nitrogen or carbon dioxide.
Chemical blowing agents are special additives that decompose when they are heated or
react in such a fashion where gas is released. The bubbles grow and eventually settle and
then, as the liquid is cooling down, become solid.
In general, polymer foams have considerably lower strength and density than solid
metals. However, their low thermal conductivity coupled with their much lighter weight
make them an ideal insulator for the ET. Table 21 (Gibson and Ashby, 1988) lists
several properties of polymer foams and true solids (solid metals and ceramics).
Table 21. Properties of solids and polymer foams
Density (lb/in3) Conductivity (BTU/hr in F) Young's modulus (psi)
Solids 3.6 x 102 3.6 x 10 4.8 x 10 48 3 30 x 106
Polymer 3.6 x 104 4.8 x 103 3 x 103
foams
Initial Foam Material Testing
The Columbia disaster forced NASA and all those associated with the application
of the foam insulation material to the ET to see it in an entirely different light. Numerous
tests were conducted at cryogenic, room, and elevated temperatures to ascertain the
properties of this material. In addition to standard torsion and tension tests to obtain the
elastic moduli, various tests were performed to evaluate the failure characteristics, such
as the plane strain fracture toughness, KIc. While these tests are an important first step in
understanding how the foam behaves on a fundamental level, NASA also took extra steps
to examine how the foam sheds from the tank through a specialized test covered in a later
section.
Determination of Elastic Constants and Material Classification
Some of the initial work toward understanding the mechanical properties of foams
took places in the 1960s and 1970s. Gent and Thomas (1959), and Patel and Finnie
(1970) among others were among the first to explore the properties of cellular materials.
The model developed by Gibson and Ashby (1982), however, is perhaps the most widely
accepted way to analyze the properties of cellular solids. Their analysis is rooted in a
unit cell comprised of a network of interconnected struts (Figure 21) The equations that
arise from this mathematical treatment are in terms of the strut dimensions, which in turn
Figure 21: Idealized foam microstructure
allow some material properties to be related to the relative density. The relative density
is the density of cellular material, p*, divided by the density of the solid from which the
cell walls are made, ps.
Generally, foams can be divided into two main categories: those with open or
closed cells. By closed we mean materials where each cell is partitioned from its
neighbor by membranelike faces. Within opencelled structures the individual members,
struts say, are individually connected. With Gibson and Ashby's (1982) model, many of
the equations that determine various properties, such as the elastic moduli or the fracture
toughness are derived considering the configuration of the cell walls.
Understanding the mechanical and metallurgical aspects of these struts at a
microscopic level is very important. After all, it is these struts that transmit the applied
mechanical loads and transfer heat energy. But from a practical standpoint, models that
contain these microscopic parameters are not very useful to the engineering community.
The engineer may not have access to the equipment and facilities necessary to measure
things like cell wall thickness. Generally, the analyst, or engineer, may only be privy to a
quantity like the relative density. This parameter is one of the most important, and
useful, concepts that help define and understand the properties of cellular materials. The
relative density can be related to the dimensions of the cell walls. For honeycombs we
have
p* t
=c (21)
P, I
where t is the cell wall thickness and the 1 is the edgelength. For opencell foams
= C2 (22)
P, 1
Finally, for the foams of the closedcell variety we have
p* t
=C3 (23)
P, I
The C terms are numerical constants that depend on the cell shape, normally taking a
value of unity.
As mentioned beforehand, certain opencelled foams can be modeled as a network
of beams. Using an idealized array of connected beams, the theory in many solid
mechanics texts such as Timoshenko and Goodier (1982) is adequate to determine the
deflections, strains, and stresses. It can be shown that the Young's modulus for the open
celled foam is given by
E* C1EJ (24)
E 14
where E* is the modulus of the foam material. The relative density is related to the
dimensions, t and 1 as shown in (22). The second moment of the area, I, can also be
related to the cell wall dimensions
I oc t4 (25)
Now (24) can be rewritten as
E C,P (26)
Using the same reasoning, we can write down the relationship for the shear moduli
P =C2 (27)
Lastly, we need v*, the Poisson ratio, which is defined as the ratio between the lateral and
axial strains. The strains are proportional to the bending deflection and their ratio is
constant. Thus, v* depends only on the cell wall geometry. We can show this for
isotropic foams where [as, the shear modulus is
E
p ) (28)
2(1+ v)
Using (28) along with (26) and (27) we have
v* C= 1 = C (29)
2C,
The above analysis is for opencelled foams. When examining closedcell foams,
the derivations are more complex because most foams are made from a liquid and surface
tension draws the material to the edges which might leave only thin membrane across the
cell. For sufficiently thin membranes, the opencelled formulas can be applied. But if an
appreciable quantity of foam is present in these membranes, the amount, or fraction, of
this material contributes to the stiffness. There are three potential contributions that sum
to the total stiffness of a closedcell material. The aforementioned cell walls make a
contribution when they are stretched. The second contribution is present due to a fluid,
usually air, trapped within the cells, and finally the cell walls, and beams, can also be
bent when a load is applied.
The contribution from the stretching is derived by considering how the applied load
bends and stretches the cell face. The cell edge bending is proportional to 12S62 where S
is the stiffness of the cell edge, and S o Esl/13, and 6 is the displacement that arises from
the applied load. The contribution from the stretching of the face is proportional to
/2EsE2Vf. The Vf term is the volume of solid in the cell face and P oc 6/1 and Vf o 12tf.
The thickness of the edges and faces are denoted as te and tf, respectively. In the linear
elastic regime, we can define the work done by the applied force as
1 aE I32 CL J12
I 13 +E. 12tI (210)
2 /f
Equation (210) can be rewritten when we note that I t4e and E* ocx (F/12)/(6/1) which
yields
=a' +/f' (211)
E, I4 1
For many foams the edge and face thickness is related to the relative density by the
following empirical relations
t =1.4(1 ~b)p*
1 p,
(212)
t 0.9301/2 0P*
I P
If we substitute (212) into (211) we obtain the ratio of the elastic moduli in terms of the
relative density and volume fraction of material in the cell edges is 0
E*= C 2 +C' (1: ), (213)
E, P C P,
where a, p, a', 3', C1, and Cl'are all constants of proportionality. The above equations
take care of the contribution from the stretching of the cell walls. The next contribution
comes from accounting for the fluid trapped in the cell walls, E*g. From Gent and
Thomas (1963) we have
Eg* p(12v*) (214)
where po is the initial gas pressure. Gibson and Ashby (1988) note that if this pressure
happens to be equal to the atmospheric pressure, this contribution is small. Lastly, shear
modulus for the closedcell foam is shown below
G=C22 p +2(1) (215)
E P, P,
Finally the last contribution to consider comes from the bending on the cell struts, and
that analysis is identical to the opencelled formulation presented earlier.
We note that the above relations are for foams where the cell walls are equiaxed.
Most manmade foams are anisotropic. When this particular foam is sprayed down,
foaming agents cause it to rise and the cells are stretched in the rise direction. Cell shape,
then, can significantly impact the material properties. The shape anisotropy ratio, R, is
R1 
L2
(216)
R13
T3
where L, L2, and, L3 are the lengths of principal cell dimensions and L, is the largest of
the three. Following a similar analysis for the bending stresses within the beams for the
isotropic lattice Gibson and Ashby (1988) go onto define a quantity called the Young's
modulus anisotropy ratio for opencelled foams as
E *3 2R2
3 2(217)
E*1 (1+(1/R)3)
2R
For closed cells an additional term (1 +) appears in (217). For the
1+ (1 / R)
anisotropic shear modulus, we have
3 (218)
G*1 1+R
The Poisson ratio is once again the ratio of lateral and normal strains and, just like in the
equiaxed case, is independent of the relative density. As such this elastic constant is
dependent solely on the cell geometry. Unfortunately, for the anisotropic foams,
employing the same type of dimensional analysis used to determine the elastic moduli for
the Poisson ratio will not offer any insight into its dependence on cell wall geometry and
Gibson and Ashby (1988) do not discuss this in their text.
So there appear to be two ways to try and determine the properties of the foam
insulation material. The properties can be estimated using the relative density and shape
anisotropy ratio or they can be experimentally evaluated. The BX265 foam contains
closed cells, so the equations presented earlier for this particular cell wall geometry along
with corrections necessary to account for the anisotropy could be utilized for evaluating
various material properties in terms of the relative density.
Evaluating the allimportant relative density, though, may not be as straightforward
as it seems. While it is not difficult to measure the density of the foamed material (one
simply weighs the sample and divide it by the measured volume), getting a handle on the
unfoamed density for this particular kind of polyurethane foam is not so simple unless
the exact same foam is available in the literature somewhere. Also, when NASA was
preparing various foam test panels (from which the test specimens are machined) they
determined the density could vary within the specimen. Also, there are several ways in
which the insulation can be sprayed down. The foam can be applied via a mechanical
process where a machine sprays it down, or it is laid down by hand using special
equipment that delivers the insulation to the testing surface (normally an aluminum
plate).
Thus, since the density could change throughout the test panel, and even vary
depending on the process by which it is sprayed, NASA developed a test program where
a vast database of information was to be established for this particular foam material. It
was decided that standard test procedures for evaluating both elastic properties and
fracture response were to be employed, instead of estimating the properties based on
relations that require the relative density to be known beforehand.
Figures 22 and 23 display some of the various test specimens used in various
tension and torsion tests. A summary of typical stressstrain data is shown in Figure 24.
Within this chart the results of several tests performed with different foam orientations
(denoted as local rise, spray, and axial) at a few temperatures (RT denotes room
temperature and LN2 is the temperature of liquid nitrogen).
Figure 22. Square specimen tensile test used to determine the Poisson ratio
As expected the colder the temperature, the more 'stiff the foam becomes and in general
the foam fractures with little deformation. It should be pointed out that the foam is not a
material that behaves exactly like a generic isotropic steel alloy. Most engineering
Figure 23. Foam test specimens
16
Summary of Typical Tensile Stress Strain Curves
BX265 / Category 4
002 004 006 008 01 012
Strain (inlin)
Figure 24. Summary of stressstrain data t
metals, when placed in a uniaxial tension test, have a linear relationship between stress
and strain up until the yield point. Cellular materials generally have a nonlinear
relationship between stress and strain and the constitutive relations defined in this chapter
assume that this material has a linear stress strain relationship because a model that
accurately describes the response for foam has not yet been developed.
The aforementioned tension and torsion tests have led NASA to classify the foam
as a transversely isotropic material. These materials are sometimes called hexagonal
materials and five independent elastic constants relate the stresses to the strains in the
constitutive matrix. The 1122 plane in Figure 25 is a plane of isotropy.
Source: personal communication, Doug Wells (MSFC)
Figure 25. Coordinate system for the transversely isotropic material
Using the coordinate axes from Figure 24 we can define the constitutive relations for a
transversely isotropic material. The Young's moduli in the 11 and 22 directions are
defined as being equal as are the shear moduli in the 33, or rise, direction. The other
shear modulus, G12, is determined through a relation with Ell and v12: G12= 2(1+ v12)/Ell.
Table 21 lists the room temperature values of the elastic constants for the orientation
shown in Figure 25.
Since one of the objectives of this study is to analyze the effect of temperature, or
thermal loads, on the Ksolution the material properties at various temperatures are
needed for a meaningful analysis. NASA evaluated the foam's elastic moduli, thermal
conductivity, and coefficient of thermal expansion.
E E2 E3
21 1 0 0 0
Si E E,2 3 011
o22 31 32 0 0 0 (219)
033 El E2 E3 033
723 0 0 0 1 0 0
713 G23 13
712 0 0 0 0 1 0 "12
G13
0 0 0 0 0
G,,
where the five independent elastic constants are
E. = E ;E
E11 E22; E33
G23 G13 (220)
V13 V31
,3 12 21
E3 El
Table 21. Measured material property values from experiment
En = 950 psi v12 = 0.45 G12 = 328 psi
E22 = 950 psi v31 = 0.3 G23 = 231 psi
E33 = 2400 psi v3 = 0.3 G31 =231 psi
These values are plotted in Figures 26 through 29. The tests performed by NASA
indicated that the Poisson ratios did not vary substantially with temperature. As of this
writing, only one value for k, the thermal conductivity, is available. It is possible for that
value to vary with direction as well as temperature. It will be assumed, however, that k is
isotropicc,' or has no direction dependency.
0Ell
Young's moduli  E22
 E33
6000
5000
4000
'. 3000
2000
1000
0
550 450 350 250 150 50 50 150 250
temperature (F)
Figure 26. Young's moduli vs. temperature
Source: personal communication, Doug Wells (MSFC)
19
Shear moduli
500 400
300 200 100
temperature (F)
0 100 200 300
Figure 27. Shear moduli vs. temperature
Secant coefficient of thermal expansion
0.00016
0.00014
0.00012
0.0001
I 0.00008
0.00006
0.00004
0.00002
500 400 300 200 100 0
temperature (F)
100 200 300
Figure 28. Coefficient of thermal expansion vs. temperature
Source: personal communication, Doug Wells (MSFC)
600
500
400
u' 300
200
100
 G12
G13
 G23
as11
 as22
 as33
Thermal conductivity
0.0018
0.0016
0.0014
0.0012
0.001
3 0.0008
I
0.0006
0.0004
0.0002
0
0 
500 400 300 200 100 0 100 200 300
temperature (F)
Figure 29. Thermal conductivity vs. temperaturet
Fracture Testing for the BX265 Foam Material
While many foam materials do not exhibit linearelastic response before yielding,
brittle foams normally exhibit linearelastic behavior (while in tension) up until fracture
takes place. Thus, we would like to use LEFM concepts to compute the neartip stresses.
These analytical expressions assume that we are dealing with a continuum. So if the
crack is very small, cell size could influence these computations. Brittle fracture of
foams has been studied by Fowlkes (1974), McIntyre and Anderton (1978), Morgan et al.
(1981), Maiti et al. (1984), and Huang and Gibson (1991) among others. Perhaps the
most widely accepted theory for the brittle facture of foam is the Maiti et al. (1984)
study. They present a relation to compute the toughness for brittle foams as
KIC = C8f (Tl) 2 (p*/ )3 2 (221)
where Cs is a constant of proportionality and cfs is the failure strength of the struts.
Equation (221) can also be rewritten as
KI, =C, KIJs (p/p* )3/2 (222)
where KIcs is the toughness of the solid material. Equation (222) assumes that KIcs is
proportional to Cfsl1/2 and once more the desired property is related to the relative density.
The concern, now, is at what length scale is such an equation valid.
Huang and Gibson (1991) analyze the fracture toughness of brittle honeycombs.
The primary aim of their paper is to examine the effect of short cracks and determine
some sort of correction factor for (221). They go on to show that cell size does indeed
impact the fracture toughness calculations. Using fractographic analysis, Brezny and
Green (1991) study short cracks in oxidized carbon foams. Their results dovetail with
Huang and Gibson's (1991) in that cell size can influence the fracture response.
However, in both studies short cracks are used in the various experiments. By
short, we mean crack lengths less than an order of magnitude greater than the average cell
size. But in Huang and Gibson's (1991) paper, the wellknown relations based on
continuum assumptions seem to be valid as long as the crack length is an order of
magnitude larger than the cell size. Thus, if one is considering cracks much greater than
the size of the cell lattice, it seams reasonable to model the foam as a continuum.
To that end NASA conducted numerous singleedge bending SEN(B), compact
tension C(T), and middle tension M(T) tests to examine how the foam fractures under
load and used the regular, continuumbased relations to determine fracture properties,
such as the plane strain toughness. Pictures of these test specimens are shown in Figure
210. Since this material is anisotropic, each time the fracture toughness is calculated,
the material's orientation is supposed to be reported along with it. Figure 211
summarizes some of the tested geometries at various temperatures (left and right sides of
the vertical lines), for a few orientations, denoted by the labels along the bottommost
horizontal axis.
*0
Figure 210. From left to right: M(T), C(T), and SNE(B) fracture test specimens
While NASA utilized many geometries to test the fracture toughness, the C(T) test
was the most extensively tested specimen used to estimate a toughness. The plane strain
fracture toughness is usually obtained with slow loading rates (Barsom and Rolfe, 1999).
By plane strain, we mean thick plates, or test specimens, with deep cracks. The foam is
brittle and fracture is sudden with little or no stable crack growth. Wells (1966)
suggested that the fracture behavior near the crack tip can be characterized by crack tip
opening displacement (CTOD). He showed that CTOD is analogous to the crack
extension force (sometimes referred to as Go in the literature) and CTOD could be related
to the plain strain fracture toughnessKIc.
If the material is brittle and the subsequent load vs crack mouth opening
displacement (CMOD) curve (Figure 212) is linearelastic and the C(T) specimen meets
the standards mandated by the ASTM, KI, can be calculated from this test method and
this approach is that NASA used to calculate the fracture toughness for the foam. In
order to determine KIc, a preliminary fracture toughness, KQ, is determined by the
following relationship (Anderson, 1991).
Where PQ (Figure 213) is maximum load applied to the C(T) specimen before
failure, B is the specimen thickness, W is the width of the specimen measured from the
centerline of the pin holes to the rear of the specimen, and f(a/W) is a dimensionless
polynomial function. The ASTM E399 standard denotes which function to use
depending on what specimen one is using. If the a/W ratio and B (thickness) are correct
per ASTM procedure and as long as Pmax is < 1.10PQ then KQ = KIc.
BX265 / Geometry and Constraint Overview I CAT 4
30.0 320F RT 320F R RT 320F 20F RT
25.0 
20.0   
I I
a <* i 0.
$ g
15.0   
0.0
SEN(B)(AS) M(T)(A) M(T)(33S) C(T)(A.S) C(T)(33S)
Figure 211. Summary Kvalues for various fracture specimens
P
K B f (a/W) (25)
Clip
A CMOD
Figure 212. Clip gauge on C(T) specimen used to measure CMOD
Additional Testing: Divot Test Specimens
While NASA was conducting tests to evaluate the material properties at various
temperatures, additional experiments were devised to examine how and why the foam
becomes detached from the tank. Gibson and Ashby (1988) discuss the design of
sandwich panels with foam cores. These types of structural members are normally
comprised of two, stiff, skins separated by a lightweight core. These parts are utilized
heavily by the aircraft industry, particularly in applications where resistance to bending
and buckling loads are important, such as helicopter rotor blades. One mode of failure
for these types of structural members is decohesion of the adhesive bond between the
skin and the core.
Load, P
Pmax = PQ
/
A CMOD
Figure 213. Typical load vs CMOD curve for a brittle material
Often the epoxy bond is stronger than the core material itself. So if the interface
between the skin and core is defectfree, delamination is usually not a concern. The
situation changes if there are defects within the interface, however. This type of analysis
is complicated by the difference in elastic constants between skin, adhesives and core.
The delamination described by Gibson and Ashby can be likened to 'weakening' the
bond between the core and the outer skin. When the strength of the bond becomes too
weak, the foam core and skin can peal apart. The failure mode that NASA is
encountering, however, is more of an explosive and sudden 'blowout' of foam from the
ET.
Since the foam loss seems to be the greatest near areas where the surface of the
tank is somewhat uneven (from bolts and fittings), NASA believes that voids created
when the foam is first being sprayed down are the primary reason why the foam is able to
break off. These are indeed defects between the ET, or skin, and the foam, or core
material. To model this phenomenon, NASA devised an experiment to examine how
various voids within foam test panels could contribute to large pieces of foam being
blown out under certain conditions. Rectangular test panels like the one shown in Figure
214 are used for the 'divot testing.'
Figure 214. Divot test panel
In this experiment a cylindrical bore is machined into the panel. At the top of that
bore, a sharp notch is created to simulate a crack, or defect, near this void. Liquid
nitrogen is then pumped into the bored hole and a heat flux is applied across the other
side of the foam panel. Twodimensional schematics of this process are shown in Figures
215 and 216.
The applied flux warms the liquid nitrogen and eventually a phase change takes
place. The experiment takes place in a thermovacuum chamber to simulate the
environment the foam is exposed to during liftoff. As the Shuttle ascends to space, the
atmospheric pressure is dropping, but the gaseous N2 applies pressure on the walls of the
void and the crack faces. With enough force, the flaw turns and propagates toward the
surface.
Figure 215. 2D view of divot test
In Figure 217, we see the aftermath of one of these experiments. A frustum
shaped piece of foam has been blown out and the crater left behind resembles a 'divot,'
similar to the removal of a small piece of turf after a golf shot has been played. It is here
that predicting the crack turning angle could have a practical application.
Heat Flux
O\\\ \,\\
X Crack
Figure 216. 2D view of test panel after heat flux has been applied
For a given void and flaw size, along with the applied tractions, it might be possible to
calculate this angle, denoted as coc in Figure 216, by calculating neartip stresses and
applying a turning theory such as the maximum tangential stress criterion. The concept
of the crack turning angle and methods to compute it will be covered in a later chapter.
Sharp Notch
Void
Liquid N,
Figure 217. Foam blow out from divot test
Summary
The pertinent theory (Gibson and Ashby, 1988) regarding cellular solids is
presented and various formulas are available depending on what type (open or close
celled) of foam one is analyzing. Formulae for anisotropic foams are presented as well.
The elastic constants, for both isotropic and anisotropic constants can be estimated
through the relations presented in the early sections of this chapter or via experiment.
While there is extensive literature available on foams, NASA decided to perform tests to
evaluate the elastic constants and fracture properties at various temperatures. This
enabled NASA to have a large collection of experimental data for several material
orientations, over a wide range of geometries.
Once the tests were conducted, determining the fracture properties can be done
using the wellknown empirical relations for standard test specimens. One of the main
assumptions behind these equations is that they are to be applied to a continuum. The
foam is a porous material and it has been shown that cell size can impact fracture
properties if one is conducting analyses with short cracks. However, continuum
assumptions are acceptable as long as the crack is an order of magnitude longer than the
cell length.
Lastly, additional and highly specialized tests were performed to investigate how
the foam is able to free itself from the tank. NASA believes that voids near fittings and
bolts allow liquefied air to collect in these cavities prior to liftoff. During liftoff, the air
rushing over the outside surface of the insulation warms the liquid nitrogen which
converts it into a gas. The pressure, while small, appears to be sufficient to drive a crack
toward the exterior and this seems to account for the foam loss outlined in chapter one.
An idealized case of this process is examined through the divot test. These
experiments entail machining oval or cylindricalshaped voids into a foam test panel and
pumping in liquid nitrogen. The cryogen is heated when a flux is applied to the top part
of the foam panel. At the top of the void, a sharp notch is inserted and when the nitrogen
is turned into a gas, the pressures exerted on the void walls and crack faces are sufficient
to drive a crack toward the surface. It is here that a determination of the crack turning
angle could be of some use to NASA and there a few theories prevalent in the literature
on this topic. Most of them require the evaluation of neartip stresses and/or strains.
Since this an anisotropic material, the general equations for isotropic materials are not
applicable in most cases. In the next chapter a detailed discussion covering linear,
elastic, anisotropic fracture mechanics is covered and these concepts are applied in
chapter four when crack turning theory is presented.
CHAPTER 3
AN INTRODUCTION TO FRACTURE MECHANICS
Material Definitions
For local crack tip calculations there are material definitions, or classifications, that
need to be defined before moving on with discussions on how to compute neartip
stresses for materials that are anisotropic, or have directiondependent properties. In
general, most texts on fracture mechanics focus on materials that are isotropic. Here the
constitutive matrix contains three elastic constants, two of which are independent, E and
v. For isotropic materials, there are an infinite number of planes of material symmetry
and the strains and stresses are related to each other via equation (31).
E 1/E v/E v/E 0 0 0o ~ x
Eyy 1E v/E 0 0 0 a
S1/E 0 0 0 C (31)
y" 2(1+v)/E 0 0 a
,, 2(1+v)/E 0
7, sym 2(1+v)/E_ c
Many other types of materials can be characterized by the number of planes of
internal symmetry (Dieter 1976). Cubic materials, for example, have three independent
material parameters and nine planes of symmetry. Materials with three internal planes of
symmetry are known as orthotropic materials and many engineering metals fall into this
category, such as rolled aluminum. A special class of orthotropic materials has a single
plane of symmetry that is also isotropic. These are known as transversely isotropic
materials and the constitutive relations for this material were defined in the previous
chapter. Monoclinic materials possess a single plane of material symmetry, but the
behavior in that plane is not isotropic.
Being cognizant of the constitutive relations for the particular material in question
(be it isotropic or fully anisotropic) is important because many of the initial derivations of
neartip stress fields for anisotropic materials tried to take advantage of material
symmetry to decouple 'inplane,' and 'outofplane' displacements, which, in turn, makes
the mathematical formulation a little less complicated.
Isotropic Fracture Mechanics
For isotropic materials, the neartip stress fields have been analyzed by
Westergaard (1939), Sneddon (1946), Williams (1957), and Irwin (1960) among others.
These solutions pertain to specific cracked configurations with applied farfield loads and
are governed by a single parameter, K, the stress intensity factor (Anderson, 1991).
K
0, (r, o) = f,j (o) + higher order terms (31)
y
P(x,y)
Crack faces
Figure 31. Coordinate system used in equation (31)
The coordinate system for equation (31) is shown in Figure 31, where the fij term is a
trigonometric function. We note the Greek letter, theta, is normally used to define the
angle with the horizontal in the polar coordinate system presented above, but theta will be
used later on to define the material orientation. A different variable name is selected to
avoid confusion. The key points about (31) are that the neartip stresses have no
dependence on the elastic constants and the first term is proportional to 1//r and,
therefore, the stresses approach infinity as r is made smaller and smaller. Also, as r
approaches zero the 'higher order terms' in (31) get smaller or remain at some finite
value. Thus, in many texts and papers the analysis is restricted to small values of r such
that the higher order terms are neglected.
There is typically a roman numeral written next to the stress intensity factor to
denote the mode of loading and Figure 32 displays the three modes of deformation.
Figure 32. From left to right: mode I (opening), mode II (shearing), mode III (tearing)
While the stresses are proportional to 1/Nr, the displacements are related to
(32).
0) 2 0) 2 0
fu2 2 2 2 KI
x 2, rcos K 1+2sin2 0) sin K+1+2c05 C o3
uy =s sin fK+l2cos2 cos2 Kl2sin2j 0 KII (32)
u 2, 2 2s 2J 2l 2( 0
Zu KIII
0 0 sin
2
where K = 34v for plane strain and K = (3v)/(1+v) for plane stress and ps is the shear
modulus.
Linear Elastic Fracture Mechanics for Anisotropic Materials
In general, the equations in the previous section cannot be applied to this problem
directly because the material properties are directiondependent. The analytical work
geared at developing the neartip stress fields for elastic anisotropic materials has been
examined by Sih et al. (1965), Bowie and Freese (1972), Barnett and Asaro (1972), Bogy
(1972), Tupholme (1974), Kuo and Bogy (1974), Rathod (1979), Hoenig (1982), and
Dhondt (2002) among other researchers.
The important work of Sih et al. is one of the most cited references for determining
the stress fields near the crack tip for anisotropic materials. Following Westergaard's use
of stress functions and complex variables to determine neartip stresses for isotropic
materials, Sih et al. use a similar approach for materials with directiondependent
properties. Whenever stress functions are utilized to solve this kind of problem, there is
normally a corresponding characteristic equation (a detailed explanation of the
characteristic equation can be found in Appendix A). For the most general case, the
order of the characteristic equation is six. In Sih et al.'s paper the material in question is
monoclinic and the crack is resting in the symmetry plane. The constitutive relations for
monoclinic materials are shown in equation (33).
EX S1 S,, S3 0 0 S16 x
EY S22 S23 0 0 S26 (y
e, S33 0 0 S3 (33)
Yyz S44 S45 0 O
7Y55 0 O
7, sym S66
This configuration is advantageous because the inplane and outofplane
displacements become decoupled, which in turn reduces the order of the characteristic
equation to four.
Re _L 2 LRe[2 1 Re
Re 17r2 Rer i r1 1
Re I 1 _1) Re [_ 1 1 P
A 0 /2 0 Q1 Q A lP2)Q1 Q2
0 0
C. = (S31 + S32, + S36 )/S33
where the elastic constants in compliance form are defined as
, = Sc j
Re[ r (/ipP2Q P2 )] Re( 2 ('2p, P1)]
Re 1 '2 2q q qQ)] Re [ (q2Q2 ql)
0 0
0
Re Q3
C45 + P3C44
KI (37)
KII
[Kill]
The roots l1i, [l2, and [l3 for this particular formulation come from the special case for
monoclinic materials listed in Appendix A.
P, =S',11 u S',16 u +S ', (38)
q, = S',12 / S'26 + S'22/ (39)
Sih et al. enforce a plane strain condition and the reduced compliance matrix becomes
1
: 2 r
IyyI
y0([7
(7,
CyY
0
0
0
Re ii
Re[1
.a3
KI
KII (34)
KIJJ
(35)
U
/dzI
(36)
S' = S, (310)
S Y S33
Q, = Jcos ) +/ sin (311)
The C45 and C44 terms in (37) are the elastic constants which are in stiffness form
cr = Cec (312)
where
C, = S, (313)
The aforementioned decouplingg' of the antiplane displacement is now evident in
equation (37). The Uz displacement does not depend on KI and KII. This is similar to
the isotropic case in equation (32); the z component of displacement only depends on
Kill.
When analyzing the foam problem for NASA, a crack, or defect, can be at any
orientation relative to the insulation. In that sense, it does not seem very practical to
employ the Sih et al. solution unless we know the crack always lies in the symmetry
plane. More general solutions, ones that do not require a decoupling of the antiplane
displacements, will be used instead.
The asymptotic solution derived by Hoenig is the one we will adopt to get the near
tip stress fields. Materials with virtually any orientation, irrespective of the location of
the crack front, can be modeled with this method. Hoenig uses the pioneering work of
Lekhnitksii (1963) to derive the stresses and displacements given below by using
complex variables and stress functions to derive a coupled set of differential equations
whose solution entails determining the roots of the ensuing sixth order characteristic
equation. We notice that K once again arises as a scalar multiplier just as it did for the
isotropic solution.
1 A3 2NJ 1K
x = Re /cos + u sin Ro
1 3 'K
0 = ReZ
,!2,r ,=1 /cos o) +/ sin a
1 3 / A,uN, Z K
Se= uNRe J'K (314)
v 2;r ,=1 /cos cw + sin o
1 3 N
" Re AN, 1Kj
/2,r ,=1 /cos o0 + ,u sin oa
1 e N, 1KK
r Re
x 2 r ,=1 ,/cos 0 + /, sin o
S= (31 + S32yy + S34yz +S350, + S36 x )/S33 (315)
u = Re m NJlK, N (cos0 + j sin)) (316)
Here an important distinction is made between the general case and Sih et al's
formulation. Notice in (316) how the mode IIII SIFs are to be summed over all three
modes of displacement. Thus, all three components of displacement are dependent on
KIKIII.
n 1, = S1 ', 2 S'16+S' +A, (S'15 S'14)
=' S\ 24q (317)
m,, = S', S'26 S2 +A, S S,24(317)
m3, = S 41 /, S'46 42 S45 S'44
/A /A
The u, are the roots that occur in conjugate pairs and they arise from the sixth order
polynomial, equation (318). The coordinate system used in equations (314) and (316)
is the same as the one used in (31).
14 (/) ()32 (/)=0 (318)
Where
1 ) = SLp2 2S'4p + S44
13 ()= S 3 (S'4 + S6)2 + (S5 + S6) S4
4 (p) = S[4 2S'6p3 + (2S2 + 66)2 2S6p + S
The N matrix in equation (314) is defined as
N, = p 2 J3 (320)
Where j is defined as 13 (02)/l(u,).
There are two situations, also called degenerate cases, where this solution is
invalid, however. If antiplane shear and plane strain displacements becomes decoupled
for monoclinic materials, sayHoenig's formulation falls apart. He does present an
alternate solution for monoclinic materials which turns out to be consistent with the Sih et
al solution mentioned earlier. His solution also encounters difficulties if the crack is
lying in a plane that is isotropic. The Nmatrix, equation (320), becomes singular
because l1 and 02 = i. However, we have neartip solutions for the isotropic cases.
These are shown in equations (31) and (32).
Determining the Stress Intensity FactorK
In the above sections, it is necessary to develop solutions for determining the stress
intensity factor, K. Over the past 50 years, numerous methods have been developed to
determine this parameter. For simple geometries, a handbook (Tada et al, 1973) of K
solutions could be consulted. Alaibadi and Rooke (1991) present an excellent summary
of the literature on the subject of numerical evaluations of SIFs. They divide many of the
available methods into three main categories, denoted as stages in Table 31.
Table 31. Various methods for determining the SIF
Stage 1 Stage 2 Stage 3
Handbooks Superposition Collocating, or mapping
Stress concentration Integral transform
Stress distribution Body force method
Green's function Method of lines
Approximate weight function Finite element method
___Compounding Boundary element
Alternating technique
As one might have guessed, the procedures listed in the Stage 3 category are the
most accurate. Since computers have become less expensive and more powerful in recent
years, many of the current studies geared toward Ksolutions incorporate the use of the
finite element (FE) method or boundary element (BE) techniques. Boundary element
methods are particularly attractive because the dimensionality of the numerical problem
is reduced, i.e. 2D problems involve discretizing the lineboundary of the domain and
for 3D problems, just the surface of the domain is discretized.
This analysis will call on the FE method to help compute K, and there are a few
ways of determining it with this method, such as: extrapolation of neartip stress and
displacement fields, Rice's Jintegral (1968), strain energy, and the virtual crack
extension method.
FRANC3D Next Generation Software
The Cornell Fracture Group has developed software capable of creating crack
meshes and computing an anisotropic Ksolution. This particular method of solution
would fall under Stage 3 in Table 31, and in this case Rice's Jcontour integral along
with displacements output from FE software will be used to determine K along the crack
front. This software is very versatile in that it is not restricted to analyzing any one kind
of material; Ksolutions can be generated for isotropic and fully anisotropic materials
alike.
Since many FE models were built and the subsequent Ksolutions computed with
this software, it would be prudent to discuss the pertinent theory (BanksSills et al, 2005)
that governs this particular way of calculating K. Rice's J integral, a conservation
integral that measures the energy flux into the crack tip, is equivalent to the strain energy
release rate (denoted as G) for small scale yielding. Furthermore, it can be shown for a
plane problem that only has KI and KII that
G = J = a(K2 + KII2) (321)
Where a = (v2)/E for plane strain and a = 1/E for plane stress. One can see that the J
integral equation, as is, cannot give us the SIFs on an individual basis.
J = Wdy T ds (322)
with W = (1/2)oij8ij, also known as the strain energy density, T = ijnj and ui is the
displacement vector.
The path, F, is defined in Figure 33
y
n
~~=== Vx
Figure 33. Path of the Jintegral
The Jintegral was initially formulated for 2D problems, however since this
material has properties that are direction dependent, our analysis requires the use of 3D
FE models. The Jintegral, then, needs to be modified so it can be used in integration
that take place over volumes. Following Li et al (1985), the Jintegral can be rewritten as
iJ= cO, W8 q ds (323)
In this formulation, the crack is undergoing a 'virtual extension,' and only a small portion
of the crack front is being advanced. When using the formulation from the Li et al
paper, the Jintegral is now evaluated in a cylindrical domain centered over the crack
front.
Crack front
X2
Virtual crack advance
X3
Figure 34. Cylindrical domain and virtual crack advance for a 3D crack front
Following Figure 34 we can define incremental area of virtual crack advance as
A = Aama q (s) ds (324)
AL
The definition of J in equation (322) requires F to be very small. It is difficult to
compute stresses and strains about a vanishingly small path about the crack tip. Instead,
area integrals (Moran and Shih, 1987) are used (for 2D problems) because they are
easier to implement numerically. So, the q term in (323) is simply a mathematical tool
that enables us to recast the Jintegral into a slightly more usable form.
This virtual extension method is the modemday approach to solve computational
fracture mechanics problems. One of the main advantages of such a method is the ability
to compute the energy release rates across the crack front for 3D problems. Since the
energy release rates can vary over the front that implies K can also change with respect to
crack front position as well
JJ(s)q,(s)ds j
J = (325)
Sq,(s)ds Aq
where qt(s) is the crack tip perturbation.
Rice's Jintegral does not allow us to look at the SIFs on an individual basis for a
mixed mode crack problem. This hurdle was first surpassed by Ishikawa (1980) where
he decomposes the stress, strain, and displacement fields into symmetric and anti
symmetric parts. This allows him to use the Jintegral to compute the Ks individually.
Ishikawa's method, though, has not been extended for anisotropic materials. Yau et al.
(1980) uses the Mintegral developed by Freund (1978) to determine KI and KII,
individually, by using the idea that two separate and independent equilibrium states
(denoted by superscripts (1) and (2)) for a linear elastic, isotropic material can be related
to a third equilibrium state, denoted with no superscript. That is to say the stresses,
strains, displacements, and SIFs from these separate equilibrium states can be superposed
()= 1 + (2 KI= KI( + Kij(
S (1) + E(2) KI = KII(1) +KIY(2) (326)
S= (1) u (2) Kill =KIII) + KIII2
The Jintegral can now take the form
J= J1 + (2 +M1,2) (327)
The M(1,2) term is known as the Mintegral initially derived by Freund and it can be recast
in terms of the aforementioned equilibrium states
F (2) FUn(1) O1 q
M/(1,2) = 0(1,2) +j(2)t (2)(I'2)IJ (328)
F a 9x, a9 ] }9x
This integral describes the interaction between the two equilibrium states, and sometimes
the Mintegral is called the 'interaction' integral in the literature.
The first step in the procedure to determine the SIFs is to define the strain energy
release rate using Irwin's (1957) crack closure integral. Consider the crack in Figure 35.
Now a compressive, or closure stress, is applied in such a fashion to clamp the crack
faces down along the length 6r. The work it takes to perform this closure can be related
to the energy release rate, G. Irwin goes on to use the force displacement curve and
equate that to the work required to close the crack. Substitutions are made for neartip
displacements and stresses and eventually an integral expression is obtained that relates
the energy release rate to the stress intensity factor.
In a similar manner, Hoenig (via Equation 330) provides us an analytical
expression for the energy release for general anisotropy in terms of the SIFs and material
constants. Irwin's work was initially developed for brittle, isotropic materials; his
expression for the energy release rate includes the isotropic Poisson ratio and the Young
and shear moduli.
1
G lima (3 r, 0).u (r, r)dr (329)
0
IX
Compressive stress
Figure 35. Definition of Irwin's crack closure integral
Now, if the stress and displacement fields are substituted in (329), we have
G= = [KI m(m,,N, 1K,)+KIIIm(ml,N, 1K,)+KIIIm(m,,N, 1K)] (330)
When we use the relations for SIFs in (326), we get
1 (KI( + KI(2) 1M(K2 N ) + (KII() + KII(2))m(m2 N, 1K,)
G 2 +(KIII(1) + KIII2))m(mN,1K) (331)
The terms in (327) have bars over them because a certain definition of J is used in (3
23). Let us define an alternate form of (327) where J is not normalized with respect to
the extended area as
J = (1) + (2) + M(1,2) (332)
Now we substitute (331) into (332) and for clarity the individual terms of that equation
are listed below
2 (333)
(2) = [K(2) m(2,N (2))+ KI(2) Im(m,N1K(2)) + () I(m3,N 1K(2))]
M(1,2) = [KI Im( 1K(2)) +KI(2) Im 1K(1))
+KII) Im (mN, 1K 2)) + KI(2) Im (m1,NJ K(1)) (334)
+KIIm ( K 2)) + KI (2) Im (m3N 1K 2)
It is now evident how the decomposition of displacements, strains, and stresses into
separate equilibrium states is helpful. One can see that we have two definitions of the M
integral, equations (328) and (334). If they are equated, the individual SIFs (KI(1),
KII'), KIII1)) can be computed. But before that step, we must define 'auxiliary'
solutions for the second equilibrium state, denoted by superscripts (2a, 2b, 2c). The
terms with the (1) superscript are going be computed numerically via FE software.
The three auxiliary solutions are defined as: 2a, 2b, and 2c. For case 2a, KI(2a)= 1
and KII(2a) = KIII(2a) = 0. For case 2b, KII(2b) = 1 and KI(2b) = KIII(2b) = 0, and for case 2c,
KIII (2c) = 1 and KI(2c) = KII(2c) = 0. Since the SIFs are being prescribed beforehand, the
stresses, for those equilibrium states, can be computed from (314). Now the two
relations for the Mintegral are equated
KI(1) Im(m2,Nj1K (2))+
"(2) d(1)(2) zN1K))+
x 01 0 1 KII() m(mlN 'K (2))
(335)
A, 2 KII(2) Im(1,N 1K (l))+
KIII() im(m3, N 1Kj (2)+
KIIIl(2) m(m3, N1K (2)
From (335) we can assemble a system of equations that contain the three unknown SIFs
we are after
22m (m2A1) I (mN1+ m, +AN) I \(M +11 m2N31 )KI (1 H 'M1/A
Im (M A +m 1 2 Im (mA,) Im ( i3 KIJ(1) ={ (1,2b)/Ad (336)
Im m2^N 1 +ml^N1) =Imhn^N^ + N )
Im m2A31 +M3I1) Im MilA31 +M3 21 2Im (M3 331) KI(1) (1,2) A
Finally, we note that (336) involved the extraction of stresses near the crack tip. These
stresses are taken from the Gauss point locations within the elements that encompass the
crack front. Also, the integral in (336) is evaluated numerically using Gauss quadrature.
Computing SIFs Using Displacement Correlation
Displacement correlation (DC) is perhaps the oldest and simplest method to
compute SIFs via FE displacements. Chan et al (1970) were among the first group of
researchers to employ such a method to obtain SIFs this way. In general, the correlation
point is selected where the crack tip displacements are the largest so the relative error in
displacements is small. Another advantage in using such a method is the SIFs are
separated into the three modes defined in Figure 32. In its most basic form, crack tip
opening displacement (CTOD) can be used to determine the SIF by equating the
displacements at the points in question to some analytical solution.
From equation (32) with co = t7 and Figure 36 for a plane strain problem dealing
with an isotropic material, we have
K b (uyb ,) 2u (337)
K+1 V r
(u(uxb)2l/ 2z (338)
K+1 V r
where [s, again, is the shear modulus determined by E/2(1+v) and K = 34v for plane
strain We note that Uy and Ux are the opening and sliding modes of displacement. The
above method is sometimes called a 'onepoint' matching technique since it involves
extracting displacements from one point near the crack tip. The downside of such a
a
y
Figure 36. Correlation points typically used for displacement correlation
method is that very refined meshes are needed for accurate results. Improved accuracy is
obtained through the use of a special element specifically designed to capture the
singularity present at the crack tip.
From equation (31), we see how the stresses have a 1/'r singularity at the crack
tip. For accurate computations, we would like our FE models to also have this
singularity. Quarter point elements are one way to do this. The first attempts at
modeling cracks were done so using quadrilateral elements with the midside nodes
moved to quarter point locations (Barsoum, 1976) When this is done the desired strain
singularity is achieved. However, one downside of using this type of element is that the
singularity is only present on the edges that contain the quarter point. If triangular
shaped elements are used instead, the singularity exists both along the edges and within
the element (Anderson, 1991).
In Figure 37, we see how an eight node quadrilateral element can be transformed
to a triangle and then certain nodes are moved to the quarter point locations. Let us
derive how the required strain singularity is obtained along the edge of the quadrilateral
element with the midside nodes moved to the quarter point locations. Consider the
bottom edge of the element on the right in Figure 38. Nodes one, five, and two lie along
this boundary. In order to show how the strain singularity comes about, we will need to
list the shape functions for this particular element. The shape functions (Logan, 2002) for
the isoparametric element at the corner nodes are
N, = (1 + ,)(1 + 1,)(1 +r, 1) (339)
where i is the number of the shape functions at the corresponding node (i.e. N1 is the
shape function for node 1) and
S=1,1,1,1 (i = 1,2,3,4)
77, =1,1,1,1 (i = 1,2,3,4)
We also have shape functions for the midside nodes
N, = 1 ( 2)(1+ ) t=1, (i= 5,7)
2 (341)
N, = I(I )(l+7 72) =1,1 (i= 6,8)
The global origin is placed at the corer of the element. We will need the shape functions
at nodes one, two, and five. Setting r = 1 we have
N1 = (1 ~)
2
N, = (1+) (342)
=2 )
N5 =(I 2)
Isoparametric
quadrilateral element
Degenerated isoparametric
quadrilateral element /
SNodes moved to
quarter point locations
1/4 pt elements placed circumferentially
about the crack tip
Figure 37. Triangular quarter point elements
7 3 4
L 6
1 x 5 2
7 3
8
/ 3L/4
L/4
Figure 38. Isoparametric element and degenerated element with midside nodes moved
to quarter point locations
The global x coordinate is related to the isoparametric coordinate system via
(343)
where the repeated index is an implied summation. Setting xl
0, x2 = L, and x5 = L/4
we have
1 )L
2 4
x=Nx,
(344)
We can now solve for
= 2J 1 (345)
In the same way that we can relate the global x coordinate to the shape functions, we can
also determine the displacements using the shape functions
u = Nu, (346)
Again determining the displacements along the edge containing nodes one, two, and five
we write
2 2
u= 1( l + (l++)u,+(1_ (347)
Now we can substitute (345) into (347) to obtain
u 1+2 22 u+ 1+2 2 u2+4 5 (348)
We note that the displacement shown in (348) has a /lx term present. This is consistent
with equation (32) where the displacements vary with /r. We need to differentiate
equation (348) to obtain the strain in the x direction
u 1 3 4 1 4 U 2 4
SCJ7 q. 4 I+ /_+ u2+ U5 (349)
a x 2[xL L 2 f Zj 2L fxL L
Finally, it is shown in (349) how the strain displays the desired singularity (x1/2).
One can see why the quarter point element has been the standard for modeling cracks for
the past 30 years. No special programming, or internal tampering, with the FE code is
required. Any FE program that carries quadratic, or higher order, elements can support
these special ones with the nodes moved to the quarter points.
Let us now return to our discussion of using displacement correlation to extract
SIFs from FE models. If quarter point elements are used, equations (337) and (338) are
modified because r = L/4
KI (ub )4p 2 (350)
K+1 L
KI= (U Ux,)4p 27+l (351)
K+1 L
The SIFs extraction via the DC method need not be restricted to just one set of nodes, or
correlation points. One can also consider the displacements over the whole element. The
form of these equations depends on the shape functions of the elements. Ingraffea and
Manu (1980) derive a DC extraction method using 3D elements. This formulation is
particularly useful to the present study since 3D simulations are used to model the
fracture response of the foam material. Ingraffea and Manu use a special labeling
convention shown in Figure 310 for a collapsed 20 node brick element (Figure 39), also
known as a 15 node wedge.
3 20 7 3 20 7
515
4 10 17 1 14 *10 14
12" 2 19 16 11 2 19 15
9 ___13 9 y
1 A 18 5 1,4,12 17.18 5,8, 16
Figure 39. 20 node brick and 15 node wedge elements
The shape functions for a 20 node brick element are listed below. The local (S,r, )
coordinate system is placed at the center of the brick element. For the corner nodes at
positions i = 1,2... 8 and i, ri, (i = +1
N, = )( )(l (a, + +C, 2) (352)
8
For the nodes at positions i = 17, 18, 19, 20 and ji = 0, ri = 1, (i = 1
N, = (353)
4
At i = 10, 12, 14, 16 and i = +1, ri = 0, (i = +1
N, = (354)
4
Finally, for the nodes at positions i = 9, 11, 13, 15 and ji = 1, ri = 1, i = 0
(1+ C, ) (1 + ~ ) )(1+ ;2)
N, = (355)
4
For anisotropic materials, Saouma and Sikiotis (1986) extend Ingraffea's work for
anisotropic materials. They use the displacement relations developed by Sih et al. A
similar procedure is followed in this study, but instead of using displacement equations
from Sih et al's paper, the relations in (316) are to be employed. The equations to
determine the anisotropic Ks, using Hoenig's displacement equations, are as follows
KI
KII =[B] [A] (358)
KIII 2L
where [B] is
mIN1, (cos a + p, sin C)
[B]= Re m2 Nl (cos m+ ,u sin ) (359)
m3ljN (cosu +, sin o)
C'
x, u A BC'
r I
2C C
+12 (vF +v> Dv v u, 
C E B '
.e
'' l
Figure 310. Collapsed 20node brick element
2v, v, +2vE +v, 2v, +vcv, 2vE +Fv, VD + (4v, +v, +4v, v, + 4v, ve, 4v,, + v,,)
+ (V2 +VC 2VD F v, 2vD)
2 f T o t F (357)
HC + 2" +u, 1 .+uc,  .+uC,uD,+UD (4u,+uc+4u UF +4u,u,,4uE,+U,')
+2 (2F +UC ".UF'C C )
11 wc+ +w, 1 2.+w,,1 +wp, wD, + (4w, +wc +4w, +4w, w, 4wE, +FW,
+^ ( +W F:. '+WC' 2.. F CF
Once again, the repeated indices are an implied summation. The terms in the B matrix
only depend on the roots of the sixth order characteristic equation.
Effects of Temperature on the SIF Solution
In chapter two it is pointed out that temperature plays no small role in causing the
foam to break off from the ET. For example, when Shuttle is resting on the launching
pad, one side of the foam is exposed to the temperature of liquid hydrogen and oxygen;
roughly 3000F. The other side of the foam is exposed to air, perhaps at 75 oF prior to
liftoff. Over typical acreage sprays, the foam is approximately three inches thick. This
small thickness coupled with a large thermal gradient suggests that thermal strains could
be significant. As shown with the divot testing in chapter two, even though there are no
applied farfield mechanical loads, the thermal gradient seems to be sufficient enough to
impart stresses capable of driving a crack toward the surface of the ET resulting in a loss
of material.
To account for thermal strains, we simply add this term to our constitutive matrix.
Following Hyer (1998), the stressstrain relations can be written as
E T, (T, T7Y) SS,1 S1 3 0 0 0 a0
E2 Th (T, 7f) 21 S2 S 23 0 0 0 a22
E33 E33 (T,ef) 31 32 S 33 0 0 0 33 (360)
723 44 0 0 023
713 S5 0 13
712 sym S66 '12
where the total strain is {fll 822 833 723 713 712 Th11(T, Tref) is the free thermal strain in
the 11direction, and Tref is the reference temperature. The mechanical strains are
1nmech 11 Thll(T,Tref)
mech Th
Ell 11 \ IJTre/)
.22mech E22 E Th22 (T, Tf )
E33ech 33 Th33(T Tf ) (361)
mech
723 Y23
mech
713 Y13
mech
Y12 Y712
or, stated more concisely the total strain is
S mech thermal (362)
Total= me + (362)
If the thermal expansion is linear with the temperature change, we can rewrite (360) as
11 a AT S11 S12 S13 0 0 0 o011
E22 a2AT S2 S22 S3 0 0 0 022
33 a3AT S31 S32 S33 0 0 0 033 a
723 S44 0 0 023
713 55 0 13
12sym S66 012
The ai term is the coefficient of thermal expansion and AT is the difference of
temperature between the state of interest and the reference state.
The above equations allow us to compute thermal strains in the numerical models.
Incorporating thermal and mechanical loads will involve two separate steps because a
thermal gradient implies that a temperature distribution must be obtained before solving
the model that contains both the thermal and mechanical boundary conditions. The first
step is to solve the conduction problem using the thermal conductivity coefficient listed
in chapter two. The nodal temperature distribution is saved for later use. In a subsequent
run, now with mechanical loads, the nodal temperatures are carried over as body forces
along with any applied mechanical loads.
Summary
Summing up, the analytical expressions from Hoenig's neartip solution is
presented. His equations depend solely on the roots of the characteristic equation and the
stress intensity factor, K. Hoenig's solution is not as well known as Sih et al's paper and
as such the vast majority of studies dealing with neartip solutions for cracked anisotropic
bodies use the Sih et al formulae to compute the stresses. This is acceptable as long as
one realizes those particular equations are valid for certain configurations where a
decoupling of inplane and outofplane displacements takes place. Hoenig's solution is
more general and his formulation of the energy release rate is the one chosen to compute
the SIFs within the FRANC3D software. The DC technique presented also uses Hoenig's
equations, but only the ones that pertain to displacements.
These two methods, interaction integral and DC, are presented to show how the
mode IIII SIFs can be computed. These formulations are different in a few ways.
Computing K with the Jintegral is inherently more accurate because only a moderate
level of mesh refinement is needed for the solution. However, the implementation of
such a method is very involved. Alternatively, the DC approach is much easier to apply
but the accuracy of the solution does depend on the level of mesh refinement.
Finally, the relations dealing with handling thermal strains are presented. These
equations assume an expansion that is linear with temperature change. Handling thermal
loads is straightforward by first solving the conduction problem to obtain the nodal
temperature distribution. With that, those values are forwarded as body forces that are
applied alongside the mechanical loads.
CHAPTER 4
CRACK TURNING THEORY AND FINITE ELEMENT MODELING
Crack Turning Theory
There are three theories prevalent in the literature to predict incipient crack turning
angles for isotropic materials: maximum energy release rate, minimum strain energy
density, and maximum hoop stress. All of these theories are based on LEFM concepts
and were initially developed from experiments involving brittle plates.
The turning theory proposed by Hussain et al. (1974) seeks the direction where the
energy release rate, G, will be a maximum. The minimum strain energy density theory
(Sih, 1974) postulates crack turning in the direction where the strain energy density is at a
minimum. The maximum hoop stress theory proposed by Erdogan and Sih (1963)
predicts crack propagation in the direction whereoo, defined in equation (41), is a
maximum relative to the crack tip. Crack turning, in this study, denotes the incipient
turning angle. We define the crack turning angle, wc, below in Figure 41.
For isotropic materials, the neartip stresses (in the r, co polar coordinate system)
have been derived by Williams (1957). If the derivative with respect to theta of equation
(41) is taken and set to zero, the critical angle, oc, can be determined; equation (44).
1 C 3 a) 3
) = COS KI COS KIIsin (41)
2irr 2 2 2
S= cos KI 2+sin +KI sino2KIItan (42)
2 r 2 2 2 2
Figure 41. Definition of crack turning angle
o = 2 r cos [KIsinco+KII(3cos 1)] (43)
) = 2 tan 1l 1 + J (44)
4(KII/KI)
The maximum hoop stress theory is the easiest to implement and that is perhaps the
main reason why it is widely used for turning angle predictions. The minimum strain
energy density criterion is also very popular and there is some debate (Gdoutos, 1984) as
to which one is superior. Some might argue that it is not correct to use just a single
component of stress to predict the incipient turning angle, whereas a quantity such as
strain energy density seems to be more comprehensive as all components of stress and
strain are included. Maccagno and Knott (1989, 1992) study mixedmode facture
behavior of PMMA and lightly tempered steel alloys and in both cases the crack turning
angles are well approximated by the maximum tangential stress theory.
There have been some topical investigations for crack turning within materials that
have direction dependent properties. Buzcek and Herakovich (1985) predict the crack
extension angle for orthotropic materials and some of the recent work in this area mirrors
their ideas. They assume the tensile strength, T, of the orthotropic material varies with
direction, r. Since it is difficult to measure the strength of a material in all possible
orientations, they assume that T = f(r, 0, XT, YT) where r is the angle in a polar
coordinate system, 0 is the material orientation and XT and YT are the strengths of the
material in the axial and fiber directions respectively (they are analyzing orthotropic
composites). The equation for T, also denoted as a fracture resistance parameter, must be
independent of r if the material is isotropic. The lowest order function that satisfies these
requirements is an ellipse. This function is plotted in Figure 42 in (T,,, r) polar
coordinates. Now T, can be expressed as
T7 = X, sin2 7 + cos2 77 (45)
YT
XT XT
YT
Figure 42. Elliptical relationship for T,
The crack will turn, according to Buzek and Herakovich, in the direction where the
ratio of the hoop stress to T is a maximum. The hoop stress, sometimes called the
tangential stress (or oy), is defined in Figure 43. Consider a point P and a vector
connecting that point to the crack tip. The hoop stress is normal to this vector that
connects P to the crack tip.
Y
P
Figure 43. Definition of hoop stress
Mixed mode crack propagation in anisotropic materials is also investigated by
Saouma et al. (1987). The neartip stresses are obtained via Sih et al. (1965) and in their
analysis the maximum tangential stress theory is used to predict the angle of propagation.
Saouma et al. use a slightly different criterion for crack turning, but it is of similar flavor
to what Buzek and Herakovich propose. Instead of using strength as a fracture resistance
parameter, the plane strain fracture toughness is utilized.
Boone et al. 1987) use FE models to study crack propagation in orthotropic
materials. The idea of maximizing the ratio of the hoop stress to some strength parameter
is once again utilized to predict the crack turning angle. Chen (1999) analyzes crack
propagation in aircraft fuselages. He, too, makes use of a ratio of the hoop stress to the
plane strain fracture toughness. In this case, the material in question has isotropic elastic
moduli, but the fracture properties are direction dependent. A similar analysis is
conducted by Pettit (2000). He examines crack turning in rolled aluminum, which again,
has isotropic stiffness properties but anisotropic fracture characteristics.
Carloni et al. (2003) and Nobile and Carloni (2005) consider incipient crack
turning for an orthotropic plate under biaxial loading. The former researchers utilize
essentially the same fracture resistance parameter as Buzek and Herakovich, whereas the
latter study uses a slightly modified version of the Saouma et al. turning criterion.
Carloni et al. and Nobile and Carloni also call on complex variables and stress functions
in their formulations for neartip stresses. Both studies, however, make assumptions that
reduce the order of the characteristic equation that arises in such derivations (see chapter
three), down to order four.
While Pettit's analysis is for a material with isotropic stiffness properties, he
generalizes the fracture resistance parameter to three dimensions. His analysis accounts
for an arbitrary location of the crack with respect to the material. Thus, we would like to
use the method to obtain stress intensity factors given by BanksSills et al. (2005) in
conjunction with Pettit's formulation for a fracture resistance parameter. This approach
is a bit more generalized than the aforementioned studies because their work seems to be
confined to specific configurations that usually decouple the inplane and antiplane
displacements.
Pettit's method for determining the fracture resistance is formulated for orthotropic
materials which possess cubic symmetry, i.e. three orthogonal planes of symmetry and
within each plane lay two principal toughness values resulting in six principal toughness
values: K12, K21, K23, K32, K13, and K31. The first number in that nomenclature is the
vector normal to the crack plane and the second number corresponds to the direction of
propagation. It is difficult to obtain a toughness for all possible orientations for materials
that exhibit some form direction dependence. Essentially we are interpolating between
the six principal toughness values for the orthotropic material. This gives us an estimate
of the toughness in any direction we desire. Consider Figure 44 and the material
orientation as shown. What we see in that figure is block of foam with various M(T)
specimens oriented within it. If the M(T) specimens are aligned with the principal axes
of the large block of foam, there are six possible orientations and with each orientation
there is a corresponding toughness.
Pettit's interpolation equation, then, requires six unique toughness values for the
orientations shown in Figure 44. In our case, only three toughness values are needed,
however. In Figure 44, consider the M(T) specimens that have the K31 and K32 labels.
Here the load is applied in the 3direction and the crack is propagating in the 1 and 2
directions, respectively. Since the 3direction is the rise direction and the 12 plane is a
plane of isotropy, it is reasonable to assume that K31 = K32. That is to say, the crack is
resting in the plane of the knit line for those two cases. For the other four cases, the
fracture plane is normal to the knit line plane. Thus, K21 = K12 and K13 = K23. So the
total number of required toughness values has been reduced to three and NASA did
conduct tests, at room temperature, to obtain them: K32 = 22.4 psi/in, K12 = 17.4 psi/in,
and K23 = 19.5 psi/in.
Pettit's equation that will determine the effective K needs several variables besides
the six principal toughness values. Remember, the goal here is to try and predict the
direction of propagation. We can conveniently describe this direction with a vector,
denoted as a. That vector lies in a plane tangent to the fracture surface, whose normal
vector will be defined as n. These two vectors are shown in Figure 45.
Figure 44. Orthotropic toughness values
Crack faces
Figure 45. Definition of the a and n vectors
To develop the equation for the effective K, we first define the trace angles of the
vector a makes with the principal planes. These are
tan (,) = a3
a2
tan (o) = a
a3
tan(o93) = a
a1
The subscript in the theta term denotes the axis normal to the principal plane. Using
equation (45) we can write a 2D interpolation function of the form
(46)
Kk (,k) = K, cos2 (ok) + K, sin2 (ok) (47)
We also observe the following trigonometric identities
cos2 tan b  C
c b 2 + (48)
sin2 tan b b22
c b +c
and the components of unit vector a must satisfy the relation
a1 +a2 +a2 =1 (49)
Pettit goes on to define equations that can be interpreted as the fracture resistance
components of a in the principal planes. Consider Figure 44 once more where the 1axis
is the loading direction. To estimate an effective K on this plane, we will interpolate
between two toughness values, K12 and K13. Using (47) through (49) we can write K1
as
K, (a) (K2a22 + K,3a32) (410)
1 a2
In a similar fashion we can determine K2 and K3
K2(a)_1 1 2 (K23a32 +K21a12)
K, (a) =  Ka + Ka
a2 (411)
K3 (a) 1 2 (K31a12 + K32a22)
1 a3
Ia
Finally, Pettit sums the relations in (410) through (411). The effective K, denoted as
Kp, can now be written as
K,(a,n) = K1i2 +K2n22 +K3n32 = (K 2 2 +K,13a2)
1 a 2
S2 ( (412)
+ I22 (K23a32 +K2a2)+ n2(K3,2 + K32a22)
1a 2 1a3
FE Modeling of the M(T) Specimen
While the importance of crack turning was made evident by the divot test, this
experiment is highly specialized and is not recognized as a standard fracture test
specimen. A more conventional way of analyzing crack turning can be performed with
the middle tension, M(T), test specimen. Material orientation, mode mixity, and
boundary conditions can all be carefully controlled with the M(T) specimen.
The idea behind creating FE models of M(T) specimens is twofold. One of the
objectives of this study is to examine the effect of mode mixity on the Ksolution, or
stated differently, as the material orientation is changed how does this impact KIKIII
along the crack front? Also, several foam M(T) specimens were fabricated and the
resulting crack turning angles measured. Numerical predictions, based on Ks extracted
from the FE models are used to make predictions via the maximum hoop stress theory
discussed earlier.
M(T) FE models are built using ANSYS FE software. Crack meshing is handled
by FRANC3D, a FE software package developed by the Cornell Fracture Group. Figure
46 shows a schematic of the M(T) specimen and Figure 47 displays one of the FE
models used in this analysis. The FE model has a unit applied pressure at the top and
bottom faces and is simply supported. The model is comprised of entirely of SOLID92
and SOLID95 elements (ANSYSTM, 1999) which are used to make 3D models and have
the capability to handle anisotropic material properties. The model is twelve inches high,
fiveandahalf inches wide, and an inchandhalf thick. The boxed in area in Figure 47
highlights the densely meshed region containing the crack.
i i
2a
Figure 46. Schematic of M(T) specimen
Figure 47. M(T) model created in ANSYS
Since the foam has direction dependent properties, this needs to be accounted for in
the FE models. This is best done by using direction cosines. Using Figure 48, let us
look at simple transformation of axes to define a new material orientation. Consider two
sets of axes: x, y, z and x', y', z'. Initially, the both sets of axes are aligned and then a
rotation, anticlockwise (of angle zeta) about the x axis, takes place. Table 41 lists the
direction cosines, i.e. the cosines of the angles between the primed and unprimed axes.
Here aq is the cosine of the angle between the x' and x axis, U2 is the cosine of the angle
between the y' and x axis, and so on. We will examine several material orientations with
different initial crack inclinations. In this way, mode mixity is introduced into the
problem and we also examine the effect of material orientation on the predicted
propagation path.
K Z
Z' Y
X, X'
Figure 48. Sample transformation
Table 41. Direction cosines
x y z
x' ai Pi Yi
Y' Ga2 32 Y2
z' X3 P3 73
As mentioned in an earlier section, for typical acreage sprays there is little offset
between the material and substrate coordinate systems except when the foam is applied
near fittings and/or bolts where the foam can to rise at angles as high as 30 degrees
relative to the ET. Twelve M(T) fracture test samples with varying material orientation
and crack inclination were fabricated and tested. A broken M(T) specimen is shown in
Figure 49 and a comparison between the numerical prediction and the measured turning
67
angle, as well as how the foam M(T) specimens are fabricated, will be made in chapter
five.
Figure 49. Fractured M(T) specimen: the dotted line is added to show the location of the
knit line
CHAPTER 5
RESULTS AND DISCUSSION
Effect of Material Orientation on the SIF Solution
On most sections of the ET, the foam is sprayed on in such a fashion that it rises
normal to the surface. In that sense, there is very little offset between the tank (substrate)
coordinate system and the foam (material) coordinate system. Also, when modeling
these specimens, we should note that the curvature of the ET is rather large; some 28 ft in
diameter and our specimens tend to include small defects (usually two inches). Thus, it
seems reasonable to assume that curvature effects are minimal. But when the spraying
process takes place near a bolt or fitting point, the foam must be applied around, and over
top of, these parts. In this situation it is possible to have the foam rising at different
angles relative to the surface.
The anisotropic nature of the foam requires 3D FE models to be employed so the
SIFs can be evaluated along the crack front. To that end, the effect of material
orientation on the Ksolution will be examined within models that have both straight and
diagonal cracks. Before moving on with this discussion, let us define a few important
terms that we will use extensively throughout this chapter. The crack inclination angle,
(p, is the angle the crack makes with the horizontal as shown in Figure 51. M(T) models
use thoughcracks, or flaws that are placed through the thickness of the specimen, or
model, in question. The crack front distance defined as the length of the crack in the
thickness direction, shown in Figure 52.
a O aCT a
Figure 51. Definition of crack inclination angle
To analyze the effect of material orientation of the SIF solution, 17 material
orientation cases are defined; see Figures 54 and 55. M(T) models are constructed
using ANSYSTM FE software and the flaw is inserted using FRANC3D Next Generation.
Figure 52. Definition of crack front distance
""YAZO
A unit pressure is applied to the simply supported M(T) models and in all cases, the
crack length is two inches. A picture of the ANSYS model (with dimensions) is shown
below. The material orientations are also input into the FE models along with the loads
Figure 53. M(T) model built in ANSYS
z
x
Figure 54. Definition of material orientation
Figure 55. Top view of the cones shown in Figure 54
and boundary conditions. There are many different orientations at which the foam can
rise relative to the tank. But when the foam is sprayed down, at some locations, the rise
direction relative to the tank can be as high as thirty degrees. An example of this
happening is shown in Figure 55.
Portion of External Tank
Figure 56. Hypothetical material orientation relative to the ET
One can see how the foam can rise at angle when it sprayed out the bolt shown in the
above in Figure 56. As such we would like to have a reasonable 'tolerance' that the rise
direction can fall into. This would give rise to a range, or set, of possible orientations that
could occur when the foam is being applied. Since the foam can rise at thirty degrees, it
might be interesting to examine the possibility of the foam rising at another, larger, angle;
45 degrees, for example. Thus we now have two possible cases to examine; a rise angle
that is typical near fitting points and bolts, and one that is slightly larger.
The material orientations used in this study can be visualized by setting up two
concentric cones. The points on the 'rims' of the cones denote where the z'axis (local
rise direction) will be for each material case. Let us consider an example. Point one in
Figure 54 lies on the inner, 30 degree, cone. To define the material axes for this case, a
single transformation of 30 degrees, denoted by zeta, in Figure 57 is performed. This
can be likened to the knit lines being oriented within the M(T) specimen shown in Figure
57. Table 51 lists the direction cosines, i.e. the cosines of the angles, between the
primed and unprimed axes. Here ua1 is the cosine of the angle between the x' and x axis,
c2 is the cosine of the angle between the y' and x axis, and so on.
z
Case 1 ..
Y 30
X, X
Figure 57. Definition of case one material orientation for the M(T) FE model
Table 51. Direction cosines for the case one material orientation
x y z
x' i = 1 pi = 0 Y = 0
y' C2 = 0 P2 = 0.866 Y2 = 0.5
z' 3 = 0 P3 = 0.5 y3 = 0.866
Figure 58 is a plot of the mode I SIF against normalized crack front distance for
cases 1 4. The mode I SIF for the case zero orientation, the case where the material and
substrate axes are coincident, is also plotted to make a comparison. Figure 58 pertains to
a straight crack that is parallel to the horizontal (no inclination). Even though this
material is anisotropic, for this geometry, KII and KIII remain small compared to KI.
This is not always the case, however, particularly when the crack is inclined. To examine
the effect of mode mixity, the crack in the M(T) model is inclined by 30 degrees. All
cases are rerun and in Figures 59 through 511, the mode IIII SIF is plotted versus the
normalized crack front distance for cases one through four. As expected, when the crack
is inclined KI decreases, and KIIKIII increase. This is consistent with isotropic
assumptions. For the configuration shown in Figure 61, the KI and KII are related to p,
or the crack inclination angle (Anderson, 1991)
KI= o cos2 ((P) Q.a
(51)
KII = sin (P) cos (P) /a
We can see that KII is zero when p = 0 and that KI will decrease for increasing values of
(p.
KI vs normalized crack front distance
(p =0 degrees
Case 0
Case 1
 Case 2
e Case 3
 Case 4
normalized crack front distance
Figure 58. Mode I SIF for cases 14; 0 degree crack inclination
KI vs normalized crack front distance
 Case 0
 Case 1
Case 2
Case 3
 Case 4
(P =30 degrees
normalized crack front distance
Figure 59. Mode I SIF for cases 14; 30 degree crack inclination
KII vs normalized crack front distance
(p =30 degrees
Case 0
A Case 1
 CCase 2
c Case 3
 Case 4
normalized crack front distance
Figure 510. Mode II SIF for cases 14; 30 degree crack inclination
Kill vs normalized crack front distance
Case 0
ACase 1
CCase 2
 Case 3
 Case 4
(P =30 degrees
normalized crack front distance
Figure 511. Mode III SIF for cases 14; 30 degree crack inclination
The mode I SIF for (p = 0 is listed in tables 52 and 53. Here we can see the variation of
the mode I SIF as the material orientation is changed. For these cases, the ones that
resulted in the largest KI are: 11 and 15 (2.27 psi /in). When you compare that number
to the KI for case zero, the orientation for typical acreage sprays, the largest KI is just
14% higher. This is an interesting result in that even if we 'lean' the zaxis by quite a bit,
the resulting mode I SIF is not much greater than the mode I SIF for the case 0
orientation. However, it is possible to orient this material in an infinite number of ways
and it could be possible to obtain a KI value for an arbitrary orientation that is greater
than the largest value determined in this section (2.27 psi /in).
Table 52. KI for cases 08, 0 degree inclination
Case 0 1 2 3 4 5 6 7 8
KI (psi Vin) 1.99 2.01 2.14 2.25 2.14 2.01 2.14 2.25 2.14
Table 53. KI for cases 916, 0 degree crack inclination
Case 9 10 11 12 13 14 15 16
KI (si in) 2.07 2.19 2.27 2.19 2.07 2.19 2.27 2.19
Crack Turning Predictions
We can numerically predict the crack turning angles using the maximum tangential
stress criteria since we have SIFs from FE models and the neartip stress solutions from
Hoenig (1982). As discussed in chapter four, we will use a ratio of the tangential stress,
CGo, to Kp for turning angle estimations. Let us denote the ratio of a>m/Kp as R and
where R is a maximum, that is the direction of predicted propagation. To make a
comparison, NASA has donated several broken M(T) specimens with various material
orientations. The turning angles within these specimens will be measured and compared
to their numerical counterparts.
M(T) Specimen Fabrication and Determination of Local Knit Line Axes
Fabrication of the M(T) specimens involves spraying the foam on a metal plate, it
rises and eventually a rind forms. When the foam has settled, the next layer is applied
until the desired thickness is achieved. From here, the foam is cut, or machined, from
this parent block so the knit lines are running in the desired directions. One of the foam
samples is shown in Figures 512.
Three M(T) foam test specimens are analyzed in this section. These samples have
various orientations and crack inclinations ((p values). The first step in this analysis is to
determine what the material orientation is for each specimen and we do this by measuring
Figure 512. From left to right: front, left, right, and rear sides of the M(T) specimen
the angles of the knit lines on the specimen's surface. Consider Figure 513. We need to
determine the orientation of the x', y' and z' vectors in order to define the material
orientation. This information is used to define the element coordinate system in ANSYS.
The intersection of two planes is determined by taking the crossproduct of their
normal vectors. If one wants the 'trace' angle on a particular face of the M(T) specimen,
the dot product can be used. In our case, the normal to the knit line plane is not known
beforehand. Instead we will measure two angles, 0 and F, and use them to define the knit
line plane shown in Figure 513. What we are interested in, is the trace, or intersection,
the knit line plane makes on the x and y = 0 faces. Along the x = 0 face, we need the
angle of the knit line measured from the horizontal, or yaxis. In this case the knit lines
on the x = 0 face are relatively straight; 0 is approximately zero. On the y = 0 face the
angle, F, is also measured relative to the horizontal, or the xaxis. Thus, the traces on the
y and x = 0 faces define the x' and y' material axes respectively. Using the cross product,
the axis normal to knit line plane is determined, i.e. z' = x' x y'.
Knit line plane
r X
Ie I / Iz
Y Side view
X
Figure 513. Determination of knit line plane
Local Crack Tip Computations
For crack tip computations involving stresses, there are a few subtleties that should
be addressed in this section because the material classifications discussed in chapter three
are in a 'global' sense. When looking at a foam M(T) specimen, one would say it is
indeed a transversely isotropic material. However, the postprocessing of neartip
stresses requires material properties rotated into a local crack front coordinate system.
Let us consider an example. First, assume we have the case one orientation; a material
rotation (0 = 300) takes place about the global x axis shown in Figure 513. If the
inclined crack is oriented in such a fashion so that it rests in the plane of symmetry, the
outofplane and inplane displacements become decoupled and this is a socalled
degenerate case described by Hoenig. His formulation cannot be applied for this special
case. As described in chapter three, when this situation arises (crack lying in a symmetry
plane) the usual isotropic neartip formulas are utilized.
Knit Line
Figure 514. M(T) model with inclined crack
Now set
us denote this as a 'general' case and a rotation is needed to set the properties in the local
crack front system. The way the constitutive matrices are transformed depends on how
they are defined. If they are in compliance form
{} = [S]{a} (52)
then the transformation formulas (Ting, 1996) for the elastic constants are
[S,']= [Q][S][QS] (52)
where [S'] is the matrix of elastic constants in the rotated system and
a,2 AP2 '112 12y a,/71 a/,P,
a2 2 /22 72 22/2 a2 2 2)2
Q a,2 A32 732 A3'3 a3 3 3) (53
[es]= (53)
2a2a3 2/623 2 273 /7273 + 2/73 y2 a3 +a23 a282 + 2a3
2a3a, 2/,3/31 27371 ,37"1 + 73y1 a/zy + y3a, a,31 + 3,a,
2a1a2 2/1/32 22172 712 + 2 2 a y2 + y1a2 a12 + P1,a1
The order of the terms in and (52) and (53) is based on the convention commonly used
by the composites and elasticity community where the stresses and strains are ordered as
{ ox, o, % z, T Yz, y, T, xT. If the constitutive matrix is in stiffness form
{) = [C]{}) (54)
there is an alternate definition of the [Qs] matrix
a'12 2 "12 2/7121 2a,/1 2a,/81
a2 2 /22 "22 2/2212 2a/222 2a232
[C a3 / 3 2 7232 27323 20a323 2a33 (55)
a2,a3 /2/33 2,223 /3213 +2/233 22a3 + a2213 a2/62+ 2a3P,
a3a1 /33/1 7371 A/31 + y731 3a,2 + 2"y,3 a/31 + PAa
a1a2 /01/2 7172 A/322+7/372 a'12 + 2y1a a/1A2+ A11a2
The transformed elastic constants in stiffness form are therefore determined by
C'= [Qc [C]e [ ] (56)
The two [Q] matrices are related
[Q] = [ l] (57)
Equations (52) or (56) allow us to transform the material properties from one Cartesian
coordinate system to another. The only information we need for such a transformation is
the direction cosines of the angles between the coordinate axes where we start, or the
system where the material properties are initially defined, and the system where we want
to set the properties.
The incipient turning angles analyzed in this program use the maximum tangential
stress theory. In order to apply it, the first thing we need is the hoop, or Go,, component
shown in Figure 515. The neartip stresses derived by Hoenig are in a (x, y, z)
Cartesian coordinate system. To obtain the hoop stress, we transform axes, in this case a
rotation anticlockwise about the zaxis by co degrees. This is done using standard
transformation equations (Boresi et al 1993).
Y
j ic
Figure 515. Cartesian and cylindrical stresses
Calculating co is done via
o)o = aO22 0x 22 + 22 + 2 2/ + 22yz + 2/202"zx + 2a2i2oxy (54)
Here the (X2, 32, Y2 terms are the direction cosines of the angles between the old and new
axes
Now we can calculate the hoop stress at any r and co of our choosing. The next step
is to divide the hoop stress by the effective K, Kp. We do this because not only are the
elastic constants direction dependent, but the fracture properties are as well. The crack is
predicted to turn in the direction where R is a maximum. The crack turning angle,
defined as co, is shown in Figure 516. Here we have a closeup view of a M(T)
specimen containing an inclined crack. The inclined crack offset by angle, yp, and the
direction of incipient turning is denoted by the dashed line.
Figure 516. Crack turning angle
Let us pause briefly and quickly summarize the procedure needed to compute near
tip hoop stresses and the predicted turning angle
* Enter material orientation (denoted by case number in Figure 52) into a FE model
coordinate system
* Run the model and extract the mode IIII SIFs via FRANC3D
* Rotate material properties into a local crack front coordinate system
* Compute neartip stresses using Hoenig's formulas
* Transform yy to cm via the transformation equations found in Boresi et al
* Divide yy by Kp and determine where this ratio is a maximum to predict the
turning angle, oc
Table 54. Measured and predicted turning angles (in degrees)
co, (measured) oc (predicted) o, isotropicc prediction)
specimen A 27 32 43
specimen B 9.7 21 27
specimen C 12 19 27
For the three test specimens, one of them had a crack inclination of 30 degrees
(specimen A). The other two specimens (B and C) had crack inclinations of 10 degrees.
All three specimens have different material orientations. Table 54 lists the measured and
predicted incipient turning angles. The SIFs can vary along the crack front and since the
free surface can influence the numerical Ksolution, the measured and predicted angles
are computed at the midplane. For comparison's sake, isotropic predictions for the same
M(T) geometries are also tabulated. The isotropic angle is the critical angle, 0o, from
equation (44). The experimental tests were carried out at 750F and as such, the elastic
properties at that temperature are used to predict the turning angle.
The turning angle was measured using a Brown and Sharpe coordinate measuring
machine. These devices are used to inspect various geometries such as cylinder heads or
gear boxes. This instrument uses a probe that physically touches the work piece and the
points are recorded digitally. In this case, we wish to use the points to define the fracture
planes. The angle between these planes is the crack turning angle. Within each specimen
are four fracture surfaces. The fracture angle was measured at a location near the middle
and an average of the four values is taken and listed in Table 54.
Specimen A had the best agreement with experiment whereas the other two
predictions fall somewhere between the isotropic prediction and actual values. Even
though both the direction dependence of both the fracture toughness and neartip stresses
are accounted for in these predictions, the results presented here still vary by large
amount. One possible source for the difference is how the material orientation is defined
and modeled. When tasked to predict the turning angle, NASA furnished us with the
three test specimens but the orientations were not specified beforehand. These particular
specimens were initially fabricated to analyze mode mixity, not the effects of material
orientation. As such, the orientations of the knit lines vary throughout the specimens.
When defining the knit line plane, the knit lines (denoted with black dots in Figure 512)
closest to the crack were used to define the orientation, even though the material
orientation can vary over the specimen.
Based on the results seen here, it seems that perhaps the crack is seeking a pure
mode I condition, or that the crack is seeking a direction that is normal to the highest
principal stress. In such cases, if a crack, for example, is orientated at 30 degrees relative
to the horizontal, the crack will turn 30 degrees to it is running straight across or in a
direction normal to the maximum principal stress. However, only three specimens are
analyzed in this section and to some extent, the material orientation can be more carefully
controlled. Obviously more samples would have to be analyzed to determine if the
cracks favor a direction that is normal to the highest principal stress.
Influence of Thermal Loads on the SIF Solution
FRANC3D is very versatile in that it can generate different types of crack meshes
and also compute a Ksolution for both isotropic and fully anisotropic materials but at the
time of this writing, it is unable to handle problems with a thermal gradient. This is main
reason why the DC technique outlined in chapter three is utilized. To examine the effect
of a thermal gradient throughout the M(T) specimen we will examine the case two
orientation with a straight and inclined crack. This varies the amount of mode mixity and
two different sets of temperatures are applied to the top and bottom faces of the M(T)
specimen shown in Figure 517.
In the figures that follow the primary aim is to see how the SIFs change when
temperatures are applied to the M(T) specimen. Since the Ks are evaluated using a
different method, it is useful to compare the results of these two approaches when the
same set of boundary conditions is applied. We have results tabulated for several cases at
room temperature. Selecting one of them, case 1 for example, we can compare the SIF
results for a straight crack. The midplane KI at room temperature is 2.01 psi /in. Using
the DC method we obtain a value of 1.8 psi /lin, a difference of approximately 10%. The
quality and density of the mesh does impact the alternate solution presented here. In
BanksSills et al. (2005) study a similar routine is used to check the Mintegral
computations. Recall that the interaction, or M, integral is used to compute the
anisotropic Ksolution in FRANC3D. Excellent agreement is shown between their
various schemes, but much denser meshes are utilized within their FE models.
In this study it is possible to refine the mesh, perhaps by changing the location of
the quarter point elements by some degree, but this may not be feasible because there is a
node and element limitation within this FE software. BanksSills et al. use 2D models
and that allows one to change the mesh density extensively without adding a lot of
computational cost in terms of physical model size or run time.
Consider Figure 517 with (p = 0 degrees (straight crack). Here T1 = OoF and T2 =
1000F. The same unit pressure is also applied. All SIF data from FRANC3D is at room
temperature (750F). The DC method used to compute the SIFs is able to determine these
parameters at any temperature. To examine the influence of thermal loads on the K
solution, the mode IIII SIF for case two at both room temperature and with a thermal
gradient is plotted in Figure 518. From Table 52, the room temperature midplane
value for KI is approximately 2.14 psi /in. With the prescribed thermal boundary
conditions, the midplane KI drops to 1.44 psi vin; a 33% reduction.
Also, in a few of the figures, some odd activity takes place near the free surface.
When looking at Figure 518, and the line that pertains the to KI, for the DC solution, the
third point drops down to approximately 0.946 psi /in. This happens again in the rest of
the figures that pertain to the DC solution. One possible explanation is that the free
surface along with the boundary conditions is influencing the displacement results at that
point.
Sometimes quadratic elements are prone to spurious oscillations. Within Dhondt's
(2002) paper on computing mixedmode anisotropic SIFs, he performs a smoothing
operation via an averaging technique to remove the unwanted oscillations. Likewise Man
et al (1993) encounter a similar problem when using boundary elements to model a
contact problem. Their solution involves removing the higher order elements in the
problem area, and replacing them with linear elements instead. No such averaging
technique or element modification is employed here.
There is a marked change in KI, and the other SIFs (KII and KIII) also had fairly
substantial changes for this particular orientation and geometry. In this case KIII remains
small across the crack front. But KII increases by quite a bit, just like KI does. KII at the
midplane is roughly 0.04 psi /in at room temperature. With the addition of temperature
to this model, KII rises to 0.315 psi /Yin. Granted for this geometry and set of loads, KII
and KIII still remain small when compared to KI. Let us now consider a case where we
have a larger gradient applied to the M(T) model in Figure 517. The next set of results
for this material orientation are for T1 = 200oF and T2 = 1000F.
In Figure 519 we see how the larger thermal gradient impacts the mode I SIF.
With these boundary conditions, KI, at the midplane, is 1.78 psi /in. But this time, KII
and KIII both show a fairly large increase after this gradient is applied. Here KII at the
midplane reaches 0.122 psi /in and KIII= 0.355 psi /in.
///////Figure 517. Application of thermal loads
( T
Figure 517. Application of thermal loads

Full Text 
PAGE 1
ANISOTROPIC FRACTURE ANALYSIS OF THE BX265 FOAM INSULATION MATERIAL UNDER MIXE DMODE LOADING By ERIK KNUDSEN A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2006
PAGE 2
Copyright 2006 by Erik Knudsen
PAGE 3
This thesis is dedicated to my parents, fr iends, family, lab mates and everyone else who was there for me when I needed them the most. I will never forget you!
PAGE 4
iv ACKNOWLEDGMENTS First and foremost I would like to thank my parents for their unwavering love and support. My advisor, Dr. Nagaraj Arakere, ha s been in my corner since I first came to graduate school and his contacts at NASA pa ved the way for the two degrees I now hold in engineering. In the summer of 2002 I acco mpanied Dr. Arakere to Huntsville as part of the Summer Faculty Fellowship Program During our ten week stay, I met many experts and scientists who work at Marshall, particul arly Dr. Greg Swanson, Jeff Rayburn, Greg Duke, Dr. Gilda Battista, Stan Oliver, Preston McGill, Doug Wells, Dr. Eric Poole, Dr. Philip Allen and Denise Duncan (in a later trip), and many, many others. They have helped me grow as a student and a person, and I am forever indebted. After my trip to NASA, approximately one year later, I earned a NASA Fellowship to sponsor this work. Ed Adams of the E ducation Directorate at Marshall deserves special praise here, and I thank him for bl essing me with the opportunity to, once more, work with NASA. I shared an on office Dr. Philip Allen and one could not ask for a more pleasant person to work with and learn from. The Cornell Fracture Group was also a tr emendous help in my quest to better understand and learn the principl es of fracture mechanics. Drs. Paul Wash Wawryznek and Bruce Carter are two of the nicest and most patient people I have ever known and I owe them many thanks for taking the time to help me. Last, but certainly not least, I wish to thank my labmat es over the years: Srikant, Sangeet, Shadab, Niraj, Shannon, Yogen, Jeff, TJ, and George. Without their insight,
PAGE 5
v jokes, and patience I would not be where I am today. I now consider many of them to be my best friends and I hope we can stay in touch.
PAGE 6
vi TABLE OF CONTENTS page ACKNOWLEDGMENTS.................................................................................................iv LIST OF TABLES...........................................................................................................viii LIST OF FIGURES...........................................................................................................ix ABSTRACT......................................................................................................................xii CHAPTER 1 BACKGROUND..........................................................................................................1 Problem Motivation......................................................................................................1 Program Scope and Focus............................................................................................2 2 FOAM TESTING AND MATER IAL CLASSIFICATION........................................6 Cellular Solid Background Information.......................................................................6 Initial Foam Material Testing.......................................................................................7 Determination of Elastic Constant s and Material Classification...........................7 Fracture Testing for the BX265 Foam Material.................................................20 Additional Testing: Divot Test Specimens.................................................................24 Summary.....................................................................................................................28 3 AN INTRODUCTION TO FRACTURE MECHANICS...........................................30 Material Definitions....................................................................................................30 Isotropic Fracture Mechanics.....................................................................................31 Linear Elastic Fracture Mechan ics for Anisotropic Materials....................................33 Determining the Stress Intensity Factor K ...............................................................37 FRANC3D Next Generation Software.......................................................................38 Computing SIFs Using Displacement Correlation..............................................45 Effects of Temperature on the SIF Solution...............................................................52 Summary.....................................................................................................................54 4 CRACK TURNING THEORY AND FINITE ELEMENT MODELING.................56 Crack Turning Theory................................................................................................56
PAGE 7
vii FE Modeling of the M(T) Specimen..........................................................................64 5 RESULTS AND DISCUSSION.................................................................................68 Effect of Material Orient ation on the SIF Solution....................................................68 Crack Turning Predictions..........................................................................................76 M(T) Specimen Fabrication and Determ ination of Local Knit Line Axes.........77 Local Crack Tip Computations...........................................................................78 Influence of Thermal Loads on the SIF Solution.......................................................85 Summary.....................................................................................................................93 6 CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE WORK..............95 Conclusions.................................................................................................................95 Recommendations for Future Work...........................................................................96 APPENDIX GENERAL SOLUTION FOR THE PLANE PROBLEM OF A LINEAR ELASTIC ANISOTROPIC BODY..............................................................................................98 LIST OF REFERENCES.................................................................................................103 BIOGRAPHICAL SKETCH...........................................................................................108
PAGE 8
viii LIST OF TABLES Table page 21 Measured material prope rty values from experiment...............................................18 31 Various methods for determining the SIF................................................................38 41 Direction cosines......................................................................................................66 51 Direction cosines for the cas e one material orientation...........................................73 52 KI for cases 08, 0 degree inclination......................................................................76 53 KI for cases 916, 0 degree crack inclination...........................................................76 54 Measured and predicted tu rning angles (in degrees)................................................83 55 Summary of midplane KIKII.................................................................................91 56 Summary of midplane KIII.....................................................................................91
PAGE 9
ix LIST OF FIGURES Figure page 11 Foam loss during Shuttle ascension...........................................................................2 12 Sprayon foam insulation (SOFI)...............................................................................3 21 Idealized foam microstructure....................................................................................8 22 Square specimen tensile test used to determine the Poisson ratio............................15 23 Foam test specimens.................................................................................................15 24 Summary of stre ssstrain data .................................................................................16 25 Coordinate system for the tr ansversely isotropic material.......................................17 26 Youngs moduli vs. temperature..............................................................................18 27 Shear moduli vs. temperature...................................................................................19 28 Coefficient of thermal expansion vs. temperature...................................................19 29 Thermal conductivity vs. temperature......................................................................20 210 From left to right: M(T), C(T), and SNE(B) fracture test specimens......................22 211 Summary Kvalues for va rious fracture specimens.................................................23 212 Clip gauge on C(T) specimen used to measure CMOD...........................................24 213 Typical load vs CMOD curve for a brittle material.................................................25 214 Divot test panel.........................................................................................................26 215 2D view of divot test...............................................................................................27 216 2D view of test panel after heat flux has been applied...........................................27 217 Foam blow out from divot test.................................................................................28 31 Coordinate system us ed in equation (31)................................................................31
PAGE 10
x 32 From left to right: mode I (opening), mode II (shearing), mode III (tearing).........32 33 Path of the Jintegral................................................................................................39 35 Definition of Irwins crack closure integral.............................................................43 36 Correlation points typically used for displacement correlation...............................46 37 Triangular quarter point elements............................................................................48 38 Isoparametric element and degenerated element with midside nodes moved to quarter point locations..............................................................................................48 41 Definition of crack turning angle.............................................................................57 42 Elliptical relationship for T...................................................................................58 43 Definition of hoop stress..........................................................................................59 44 Orthotropic toughness values...................................................................................62 45 Definition of the a and n vectors..............................................................................62 47 M(T) model created in ANSYS...............................................................................65 48 Sample transformation.............................................................................................66 49 Fractured M(T) specimen: the dotted lin e is added to show the location of the knit line.....................................................................................................................67 51 Definition of crack inclination angle........................................................................69 52 Definition of crack front distance.............................................................................69 53 M(T) model built in ANSYS....................................................................................70 54 Definition of material orientation.............................................................................70 55 Top view of the cones shown in Figure 54.............................................................71 56 Hypothetical material orie ntation relative to the ET................................................71 57 Definition of case one material or ientation for the M(T) FE model........................72 58 Mode I SIF for cases 14; 0 degree crack inclination..............................................74 59 Mode I SIF for cases 14; 30 degree crack inclination.............................................74 510 Mode II SIF for cases 14; 30 degree crack inclination...........................................75
PAGE 11
xi 511 Mode III SIF for cases 14; 30 degree crack inclination..........................................75 512 From left to right: front, left, righ t, and rear sides of the M(T) specimen...............77 513 Determination of knit line plane...............................................................................78 514 M(T) model with inclined crack...............................................................................79 515 Cartesian and cylindrical stresses.............................................................................81 516 Crack turning angle..................................................................................................82 517 Application of thermal loads....................................................................................87 518 Room temperature KIKIII (F 3D) vs. KIKIII (DC) with T1 = 0oF and T2 = 100oF........................................................................................................................88 519 Room temperature KI (F3D) vs. KI (DC) with T1 = 200oF and T2 = 100oF...........88 520 = 30 degrees: room temperature KI (F3D) vs. KI (DC) with T1 = 0oF and T2 = 100oF........................................................................................................................89 521 = 30 degrees: room temperature KI (F3D) vs. KI (DC) with T1 = 200oF and T2 = 100oF................................................................................................................90 522 Temperature distribut ion along the crack front........................................................92 523 KIKIII vs. normalized crack front distan ce; TG2 applied in thickness direction...92
PAGE 12
xii Abstract of Dissertation Pres ented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy ANISOTROPIC FRACTURE ANALYSIS OF THE BX265 FOAM INSULATION MATERIAL UNDER MIXE DMODE LOADING By Erik Knudsen August 2006 Chair: Nagaraj K. Arakere Major Department: Mechanic al and Aerospace Engineering The National Aeronautics and Space Ad ministration (NASA) uses a closedcell polyurethane foam to insulate the external tank (ET) which cont ains the liquid oxygen and hydrogen for the Space Shuttle main engine s. This is a type of sprayon foam insulation (SOFI), similar to the material used to insulate attics in residential construction. In February of 2003, the Shuttle Columbia suffered a catastrophic accident during reentry. Debris from the ET impac ting the Shuttles thermal protection tiles during liftoff is believed to have caused the Space Shuttle Columbia failure during reentry. NASA engineers are very interested in unde rstanding the processes that govern the breakup/fracture of this complex material from the shuttle ET. The foam is anisotropic in nature and the required stress and fracture mech anics analysis must include the effects of the direction dependence on mate rial properties. Over smoot h, flat areas of the ET the foam can be sprayed down in a very uniform fashion. However, near bolts and fitting
PAGE 13
xiii points it is possible for voids and other defects to be present after the foam is applied. Also, the orientation of the foam, as it rises from nonuniform surfaces, can be significantly different from the orientation over typical acreage sprays. NASA believes that air present in these small voids is liquefied and then turned into a gas during liftoff. While the Shuttle is ascending to space, the pressure in these cavities can become large enough to force a subsurface crack toward the exterior of the tank, thus freeing portions of foam insulation. As a first step toward understanding the fracture mechanics of this complex material, a general theoretical and numerical framework is presented for computing stress intensity factors (SIFs), under mixedmode loading conditions, taking into account the material anisotropy. The effects of mate rial orientation and mode mixity on the anisotropic SIF solution are analyzed. Cr ack turning predictions under mixed mode loading are presented. Furthermore, the in fluence of temperature gradients on the SIF solution is studied, in view of the thermal gradient present thro ugh the foam thickness on the ET. The results presented represent a quantitative basis for evaluating the strength and fracture properties of anisotropi c BX265 foam insulation material.
PAGE 14
1 CHAPTER 1 BACKGROUND Problem Motivation On February 1st, 2003, the Space Shuttle Columbia suffered a catastrophic failure during reentry. NASA has conduc ted an exhaustive investiga tion of the failure and the consensus now is the breakup was caused by a segment of BX265 foam insulation, roughly the size of a suitcase, striking the wing during liftoff. The foam impact, in prior Shuttle launches, has been known to cause imp act damage to the thermal protection tiles, but was not considered to be a serious probl em until the Columbia disaster. The damage to the thermal protection tiles on the leading edge of the wing is thought to have allowed gases to penetrate the wing during reentry, triggering a cascading se ries of catastrophic events that led to the loss of the Shuttle. It is believed that when the foam insulation is applied to the external tank (ET) that voids are created in certain areas where geomet ric discontinuities such as bolts, flanges, and fittings are encountered. Since the tank is filled with liquid oxygen and hydrogen the foam is exposed to cryogenic temperatures at the ETs surface. The air inside these voids can be liquefied and nitrogen can conde nse into these small cavities. During liftoff the outer surface of the foam is exposed to ae rodynamic heating. This heating raises the temperature of the liquid nitrogen, turning it into a gas. The pressure difference associated with the formation of the gas can cause pieces of foam to be blown out during liftoff( see Figure11).
PAGE 15
2 Figure11. Foam loss during Shuttle ascension Program Scope and Focus Prior to the Columbia disaster, the foam that insulates the tank was not understood at a fundamental level. Quantities like the el astic moduli, and thermal properties such as coefficient of thermal expansion are needed for a meaningful failure analysis of this material. This prompted NASA to initiate an intensive testing program to ascertain various properties at cryogenic, room, and elevated temperatur es. The pictures in Figure 11 show the foam breaking off areas near fit ting points. At these locations, insulation is sprayed around various parts, such as bolts, a nd it is here that th e foams orientation, relative to the tank, can change by quite a bit. For typical ac reage sprays, there is very little difference between the local coordinate ax es of the foam and those of the substrate (Figure 12). Pictures of the Shuttle during liftoff, su ch as those shown in Figure 11, show us that large pieces of foam are not really breaking away over smooth, open, areas of the tank. The problem seems to be confined to areas where the foam is applied over fitting points, bolts, and other areas where various pa rts can protrude the su rface of the ET.
PAGE 16
3 Figure 12. Sprayon foam insulation (SOFI) Right away we notice that tw o potential problems come to mi nd at these locations: it is possible for voids to be created when the foam is being laid down, and the orientation of the foam can change by a large degree. The material tests conducted by NASA have show n that this type of foam is brittle. Perhaps the best way to analyze sharp, cracklik e, defects for this sort of material is through linear elastic fracture mechanics (LEF M). For isotropic materials, the neartip stresses are completely characterized by K the stress intensity factor (SIF). The analytical Ksolutions for anisotropic materi als take a similar form, but have additional terms that come from the various elastic moduli. The formulations of neartip stresses for isotropic materials are well established and are now incorporated in many finite element (FE) software packages. Some Ksolutions for certain kinds of directiondependent materials started appearing in the 1960s and mo re generalized solutions were presented in the 1980s. Since many engineeri ng materials are not isotropic (such as rolled aluminum
PAGE 17
4 in aircraft fuselages and singl ecrystal superalloys used in high performance gas turbine engines), more attention is being paid to ward the implementation of anisotropic Ksolutions in various FE schemes. Still, many of the most popular FE software packages do not allow accurate computations of SI Fs for simulations involving anisotropic materials. Foams are inherently anisotropic, and NASA, after extensive testing, has classified this material as transversely isotropic. Su ch materials also have fracture characteristics that are direction dependent and both of th ese traits must be accounted for in the analytical and numerical analyses. Alt hough anisotropic Ksoluti ons are not readily available, the Cornell Fracture Group has deve loped a software package, FRANC3D, that simplifies the meshing process necessary for FE modeling of cracklike defects in elastic solids and also computes a Ksolution for anisotropic materials. Concurrently, an alternative method for computing the SIFs, based on the work of Hoenig (1982), for anisotropic materials has been employed and extended to account for thermal stresses and strains. It is hoped that the present study will enhance NASAs understanding of this material, since the Ksolutions used in th is analysis account fo r both the directiondependent properties and the impact of therma l stresses and strains. The following points summarize the scope of this study in a broad sense: 1. Through the use of the middle tension, M(T), specimen and the FE method a systematic rotation of material propertie s and its impact on the Ksolution is examined 2. The impact of mode mixity on the Ksolution is studied 3. Crack turning predictions that take into account both the anisotropic stiffness and fracture properties are presented
PAGE 18
5 4. The effect of a thermal gradient through the M(T) specimen and its impact on the Ksolution are presented
PAGE 19
6 CHAPTER 2 FOAM TESTING AND MATER IAL CLASSIFICATION Cellular Solid Background Information Cellular solids are generally comprised into two sets of structur al categories: 1) honeycomb layouts where the microstructure is comprised of twodimensional arrays of polygons that fill a planar area and 2) foam s, the threedimensional counterpart to honeycombs where the cells take the form of polyhedra. Within Gibson and Ashbys (1988) definitive text on cellular solids, we see how many types of materials can be foamed: metals, glasses, and ceramics can be fabr icated into cells that take on a variety of shapes and sizes. Among the most common t ypes of foam are polymers which are used for many applications including insulation fo r residential homes or disposable coffee cups. Polymer foams are usually very light and they are popular in industries where weight is of prime importance, such as in the aerospace industry. Foams are typically made by bubbling a gas into the liquid monomer or hot polymer. The gas can be introduced mechan ically by stirring or by using a blowing agent. Physical blowing agents are norma lly gases like nitrogen or carbon dioxide. Chemical blowing agents are special additive s that decompose when they are heated or react in such a fashion where gas is released. The bubbles gr ow and eventually settle and then, as the liquid is cooling down, become solid. In general, polymer foams have considerab ly lower strength a nd density than solid metals. However, their low thermal conductivity coupled with their much lighter weight
PAGE 20
7 make them an ideal insulator for the ET. Table 21 (Gibson and Ashby, 1988) lists several properties of polymer foams and true solids (solid metals and ceramics). Table 21. Properties of solids and polymer foams Density (lb/in3) Conductivity (BTU/hr in F)Youngs modulus (psi) Solids 3.6 x 102 3.6 x 101 4.8 x 103 48 3 30 x 106 Polymer foams 3.6 x 104 4.8 x 103 3 x 103 Initial Foam Material Testing The Columbia disaster forced NASA and a ll those associated with the application of the foam insulation material to the ET to see it in an entirely diffe rent light. Numerous tests were conducted at cryogenic, room, and elevated temperatures to ascertain the properties of this material. In addition to st andard torsion and tension tests to obtain the elastic moduli, various tests we re performed to evaluate the failure characteristics, such as the plane strain fracture toughness, KIc. While these tests are an important first step in understanding how the foam behaves on a funda mental level, NASA also took extra steps to examine how the foam sheds from the tank through a specialized test covered in a later section. Determination of Elastic Constant s and Material Classification Some of the initial work toward understa nding the mechanical properties of foams took places in the 1960s and 1970s. Gent and Thomas (1959), and Patel and Finnie (1970) among others were among the first to expl ore the properties of cellular materials. The model developed by Gibson and Ashby (1982) however, is perhaps the most widely accepted way to analyze the propert ies of cellular solids. Thei r analysis is rooted in a unit cell comprised of a network of interconnected struts (Figur e 21) The equations that arise from this mathematical treatment are in terms of the strut dimensions, which in turn
PAGE 21
8 Figure 21: Idealized foam microstructure allow some material properties to be related to the relative density The relative density is the density of cellular material, *, divided by the density of the solid from which the cell walls are made, s. Generally, foams can be divided into two main categories: those with open or closed cells. By closed we mean material s where each cell is partitioned from its neighbor by membranelike faces. Within open celled structures the individual members, struts say, are individually connected. With Gibson and Ashbys (1982) model, many of the equations that determine various properties, such as the elastic moduli or the fracture toughness are derived considering the c onfiguration of the cell walls. Understanding the mechanical and metallur gical aspects of these struts at a microscopic level is very important. After all, it is these struts that transmit the applied mechanical loads and transfer heat energy. But from a practical standpoint, models that contain these microscopic parameters are not very useful to the engineering community. The engineer may not have access to the equi pment and facilities necessary to measure
PAGE 22
9 things like cell wall thickness. Generally, the analyst, or engineer, may only be privy to a quantity like the relative dens ity. This parameter is one of the most important, and useful, concepts that help define and understa nd the properties of cellu lar materials. The relative density can be related to the dimens ions of the cell walls. For honeycombs we have 1*st C l (21) where t is the cell wall thickness and the l is the edgelength. For opencell foams 2 2*st C l (22) Finally, for the foams of the closedcell variety we have 3*st C l (23) The C terms are numerical constants that de pend on the cell shape, normally taking a value of unity. As mentioned beforehand, certain opence lled foams can be modeled as a network of beams. Using an idealized array of connected beams, the theory in many solid mechanics texts such as Timoshenko and G oodier (1982) is adequate to determine the deflections, strains, and stress es. It can be shown that the Youngs modulus for the opencelled foam is given by 1 4* s CEI E l (24) where E* is the modulus of the foam material. The relative density is related to the dimensions, t and l as shown in (22). The second moment of the area, I can also be related to the cell wall dimensions
PAGE 23
10 4 I t (25) Now (24) can be rewritten as 2 1**ssE C E (26) Using the same reasoning, we can write down the relationship for the shear moduli 2 2**ssC E (27) Lastly, we need *, the Poisson ratio, which is defined as the ratio between the lateral and axial strains. The strains are proportional to the bending deflecti on and their ratio is constant. Thus, depends only on the cell wall geometry. We can show this for isotropic foams where s, the shear modulus is 21sE (28) Using (28) along with (26) and (27) we have 1 3 2*1 2 C C C (29) The above analysis is for opencelled foam s. When examining closedcell foams, the derivations are more complex because most foams are made from a liquid and surface tension draws the material to the edges which might leave only thin membrane across the cell. For sufficiently thin membranes, the opencelled formulas can be applied. But if an appreciable quantity of foam is present in these membranes, the amount, or fraction, of this material contributes to the stiffness. There are three potential contributions that sum to the total stiffness of a closedcell mate rial. The aforementioned cell walls make a contribution when they are stretched. The s econd contribution is present due to a fluid,
PAGE 24
11 usually air, trapped within th e cells, and finally the cell wa lls, and beams, can also be bent when a load is applied. The contribution from the stretching is de rived by considering how the applied load bends and stretches the cell face. The cell edge bending is proportional to S 2 where S is the stiffness of the cell edge, and S EsI/l3, and is the displacement that arises from the applied load. The contribution from the stretching of the face is proportional to Es 2Vf. The Vf term is the volume of solid in the cell face and /l and Vf l2tf. The thickness of the edges a nd faces are denoted as te and tf, respectively. In the linearelastic regime, we can define the work done by the applied force as 2 2 31 2s s fEI FElt ll (210) Equation (210) can be rewri tten when we note that I t4 e and E* (F/l2)/( /l) which yields 4 4* '' f e st t E Ell (211) For many foams the edge and face thickness is related to the relative density by the following empirical relations 1/2 1/2* 1.41 0.93f s e st l t l (212) If we substitute (212) into (211) we obtain the ratio of the elastic moduli in terms of the relative density and volume fraction of material in the cell edges is
PAGE 25
12 2 11*** '1 s ssE CC E (213) where , , , C1, and C1are all constants of proportionality. The above equations take care of the contribution from the stretching of the cell walls. The next contribution comes from accounting for the fl uid trapped in the cell walls, E*g. From Gent and Thomas (1963) we have 012* * 1g sp E (214) where p0 is the initial gas pressure. Gibson a nd Ashby (1988) note that if this pressure happens to be equal to the atmospheric pressure this contribution is small. Lastly, shear modulus for the closedcell foam is shown below 2 22*** '1 s ssG CC E (215) Finally the last contribution to consider comes from the bending on the cell struts, and that analysis is identical to the ope ncelled formulation presented earlier. We note that the above relations are for foams where the cell walls are equiaxed. Most manmade foams are anisotropic. When this particular foam is sprayed down, foaming agents cause it to rise and the cells are stretched in th e rise direction. Cell shape, then, can significantly impact the material properties. The shape anisotropy ratio, R, is 1 12 2 1 13 3L R L L R L (216)
PAGE 26
13 where 1L, 2L, and, 3L are the lengths of principal cell dimensions and 1L is the largest of the three. Following a similar analysis for the bending stresses within the beams for the isotropic lattice Gibson and Ashby (1988) go ont o define a quantity called the Youngs modulus anisotropy ratio for opencelled foams as 2 3 3 1* 2 *(1(1/))E R ER (217) For closed cells an additional term 2 1 1(1/)R R appears in (217). For the anisotropic shear modulus, we have 31 12* 2 *1G GR (218) The Poisson ratio is once again the ratio of late ral and normal strains and, just like in the equiaxed case, is independent of the relative density. As such this elastic constant is dependent solely on the cell geometry. Un fortunately, for the anisotropic foams, employing the same type of dimensional analysis used to determine the elastic moduli for the Poisson ratio will not offer any insight in to its dependence on cell wall geometry and Gibson and Ashby (1988) do not di scuss this in their text. So there appear to be two ways to try and determine the properties of the foam insulation material. The properties can be es timated using the relative density and shape anisotropy ratio or they can be experimenta lly evaluated. The BX265 foam contains closed cells, so the equations presented earlie r for this particular cell wall geometry along with corrections necessary to account for the anisotropy could be utilized for evaluating various material properties in te rms of the relative density.
PAGE 27
14 Evaluating the allimportant relative density though, may not be as straightforward as it seems. While it is not difficult to measure the density of the foamed material (one simply weighs the sample and divide it by the measured volume), getting a handle on the unfoamed density for this particular kind of polyurethane foam is not so simple unless the exact same foam is available in the literature so mewhere. Also, when NASA was preparing various foam test panels (from wh ich the test specimens are machined) they determined the density could vary within the specimen. Also, there are several ways in which the insulation can be sprayed down. The foam can be applied via a mechanical process where a machine sprays it down, or it is laid down by hand using special equipment that delivers the insulation to the testing surface (normally an aluminum plate). Thus, since the density could change th roughout the test panel, and even vary depending on the process by which it is spra yed, NASA developed a test program where a vast database of information was to be established for this particular foam material. It was decided that standard test procedures for evaluating both elastic properties and fracture response were to be employed, in stead of estimating the properties based on relations that require the relative density to be known beforehand. Figures 22 and 23 display some of the various test specimens used in various tension and torsion tests. A su mmary of typical stressstrain da ta is shown in Figure 24. Within this chart the results of several test s performed with different foam orientations (denoted as local rise, spray, and axial) at a few temperatur es (RT denotes room temperature and LN2 is the temper ature of liquid nitrogen).
PAGE 28
15 Figure 22. Square specimen tensile te st used to determine the Poisson ratio As expected the colder the temperature, the more stiff the foam becomes and in general the foam fractures with little deformation. It should be pointed out that the foam is not a material that behaves exactly like a generic isotropic steel alloy. Most engineering Figure 23. Foam test specimens
PAGE 29
16 Figure 24. Summary of stressstrain data metals, when placed in a uniaxial tension te st, have a linear rela tionship between stress and strain up until the yield point. Cellular materials generally have a nonlinear relationship between stress and strain and the constitutive relations defined in this chapter assume that this material has a linear stre ss strain relationship because a model that accurately describes the response for fo am has not yet been developed. The aforementioned tension a nd torsion tests have led NASA to classify the foam as a transversely isotropic material. Th ese materials are sometimes called hexagonal materials and five independent elastic constants relate the stresses to the strains in the constitutive matrix. The 1122 plane in Figure 25 is a plane of isotropy. Source: personal communication, Doug Wells (MSFC) Summary of Typical Tensile Stress Strain Curves BX265 / Category 40 20 40 60 80 100 120 00.020.040.060.080.10.120.14 Strain (in/in)Stress (psi) 33 (Local Rise) @ RT Axial & Spray @ RT Spray @ LN2 Axial @ LN2 33 (Local Rise) @ LN2
PAGE 30
17 Figure 25. Coordinate system for th e transversely isotropic material Using the coordinate axes from Figure 24 we can define the constitutive relations for a transversely isotropic material. The Young s moduli in the 11 and 22 directions are defined as being equal as are the shear moduli in the 33, or rise, direction. The other shear modulus, G12, is determined through a relation with E11 and 12: G12 = 2(1+ 12)/E11. Table 21 lists the room temperature values of the elastic constants for the orientation shown in Figure 25. Since one of the objectives of this study is to analyze the effect of temperature, or thermal loads, on the Ksolution the material properties at vari ous temperatures are needed for a meaningful analysis. NASA ev aluated the foams elastic moduli, thermal conductivity, and coefficient of thermal expansion. 13 12 123 23 21 11 11 123 22 22 3132 33 33 123 23 23 13 13 23 12 12 13 121 000 1 000 1 000 1 00000 1 00000 1 00000 EEE EEE EEE G G G (219) where the five independent elastic constants are 33 22 11
PAGE 31
18 112233 2313 1331 1221 31; EEE GG EE (220) Table 21. Measured material property values from experiment E11 = 950 psi 12 = 0.45 G12 = 328 psi E22 = 950 psi 31 = 0.3 G23 = 231 psi E33 = 2400 p si 1 3 = 0.3 G31= 231 p si These values are plotted in Figures 26 through 29. The tests performed by NASA indicated that the Poisso n ratios did not vary substantiall y with temperature. As of this writing, only one value for k, the thermal conductiv ity, is available. It is possible for that value to vary with direction as well as temperat ure. It will be assume d, however, that k is isotropic, or has no direction dependency. Young's moduli0 1000 2000 3000 4000 5000 6000 5504503502501505050150250 temperature (F)psi E11 E22 E33 Figure 26. Youngs moduli vs. temperature Source: personal communication, Doug Wells (MSFC)
PAGE 32
19 Shear moduli0 100 200 300 400 500 600 5004003002001000100200300 temperature (F)psi G12 G13 G23 Figure 27. Shear moduli vs. temperature Secant coefficient of thermal expansion0 0.00002 0.00004 0.00006 0.00008 0.0001 0.00012 0.00014 0.00016 5004003002001000100200300 temperature (F)in/in/F as11 as22 as33 Figure 28. Coefficient of th ermal expansion vs. temperature Source: personal communication, Doug Wells (MSFC)
PAGE 33
20 Thermal conductivity0 0.0002 0.0004 0.0006 0.0008 0.001 0.0012 0.0014 0.0016 0.0018 5004003002001000100200300 temperature (F)BTU/hr ft F k Figure 29. Thermal conductivity vs. temperature Fracture Testing for the BX265 Foam Material While many foam materials do not exhibit linearelastic response before yielding, brittle foams normally exhibit linearelastic beha vior (while in tension) up until fracture takes place. Thus, we would like to use LEFM c oncepts to compute the neartip stresses. These analytical expressions assume that we are dealing with a continuum. So if the crack is very small, cell size could influen ce these computations. Brittle fracture of foams has been studied by Fowlkes (1974), McIntyre and Anderton (1978), Morgan et al. (1981), Maiti et al. (1984), and Huang and Gibson (1991) among others. Perhaps the most widely accepted theory for the br ittle facture of foam is the Maiti et al. (1984) study. They present a relation to com pute the toughness for brittle foams as 1/23/2 8*/cfssKICl (221)
PAGE 34
21 where C8 is a constant of proportionality and fs is the failure strength of the struts. Equation (221) can also be rewritten as 3/2 8'*/cCssKICKI (222) where KICs is the toughness of the solid material Equation (222) assumes that KICs is proportional to fsl1/2 and once more the desired property is related to the relative density. The concern, now, is at what length scale is such an equation valid. Huang and Gibson (1991) analyze the frac ture toughness of brittle honeycombs. The primary aim of their paper is to examin e the effect of short cracks and determine some sort of correction factor for (221). They go on to show that cell size does indeed impact the fracture toughness calculations. Using fractographic analysis, Brezny and Green (1991) study short cracks in oxidized carbon foams. Their results dovetail with Huang and Gibsons (1991) in that cell si ze can influence the fracture response. However, in both studies short cracks are used in the various experiments. By short, we mean crack lengths less than an or der of magnitude greate r than the average cell size. But in Huang and Gibsons (1991) paper, the wellknown relations based on continuum assumptions seem to be valid as long as the crack length is an order of magnitude larger than the cell size. Thus, if one is considering cracks much greater than the size of the cell la ttice, it seams reason able to model the foam as a continuum. To that end NASA conducted numerous singleedge bending SEN(B), compact tension C(T), and middle tension M(T) tests to examine how the foam fractures under load and used the regular, continuumbased relations to determine fracture properties, such as the plane strain toughness. Pictures of these test specimens are shown in Figure 210. Since this material is anisotropic, each time the fracture toughness is calculated,
PAGE 35
22 the materials orientation is supposed to be reported along with it. Figure 211 summarizes some of the tested geometries at various temperatures (l eft and right sides of the vertical lines), for a few orientations, denoted by the labels along the bottommost horizontal axis. Figure 210. From left to right: M(T), C(T), and SNE(B) fract ure test specimens While NASA utilized many geometries to test the fracture toughness, the C(T) test was the most extensively tested specimen used to estimate a toughness. The plane strain fracture toughness is usually ob tained with slow loading rate s (Barsom and Rolfe, 1999). By plane strain, we mean thick plates, or test specimens, with deep cracks. The foam is brittle and fracture is sudden with little or no stable crack growth. Wells (1966) suggested that the fracture behavior near th e crack tip can be characterized by crack tip opening displacement (CTOD). He showed that CTOD is analogous to the crack extension force (sometimes referred to as Gc in the literature) and CTOD could be related to the plain strain fracture toughnessKIc. If the material is brittle and the subsequent load vs crack mouth opening displacement (CMOD) curve (Figure 212) is linearelastic and the C(T) specimen meets
PAGE 36
23 the standards mandated by the ASTM, KIc can be calculated from this test method and this approach is that NASA used to calcula te the fracture toughness for the foam. In order to determine KIc, a preliminary fracture toughness, KQ, is determined by the following relationship (Anderson, 1991). Where PQ (Figure 213) is maximum load appl ied to the C(T) specimen before failure, B is the specimen thickness, W is th e width of the specimen measured from the centerline of the pin holes to the rear of the specimen, and f(a/W) is a dimensionless polynomial function. The ASTM E399 standard denotes which function to use depending on what specimen one is using. If the a/W ratio and B (thickness) are correct per ASTM procedure and as long as Pmax is 1.10PQ then KQ = KIc. Figure 211. Summary Kvalues for various fracture specimens BX265 / Geometry and Constraint Overview / CAT 40.0 5.0 10.0 15.0 20.0 25.0 30.0K (psi(in)0.5) SEN(B)(AS)M(T)(AS)C(T)(AS)M(T)(33S)C(T)(33S) RT 320F RT 320F RT 320F RT 320F RT 320F 0.1 < a/W < 0.5
PAGE 37
24 (/)Q QP KfaW BW (25) Figure 212. Clip gauge on C(T) specimen used to measure CMOD Additional Testing: Divot Test Specimens While NASA was conducting tests to evaluate the material properties at various temperatures, additional experiments were devised to examine how and why the foam becomes detached from the tank. Gibson and Ashby (1988) discuss the design of sandwich panels with foam cores. Thes e types of structural members are normally comprised of two, stiff, skins separated by a lightweight core. These parts are utilized heavily by the aircraft industry, particularly in applications where resistance to bending and buckling loads are important, such as heli copter rotor blades. One mode of failure for these types of structural members is d ecohesion of the adhesive bond between the skin and the core. CMOD Clip
PAGE 38
25 Figure 213. Typical load vs CM OD curve for a brittle material Often the epoxy bond is stronger than the core material itself. So if the interface between the skin and core is defectfree, delamination is usually not a concern. The situation changes if there are defects within th e interface, however. This type of analysis is complicated by the difference in elastic constants between skin, adhesives and core. The delamination described by Gibson and As hby can be likened to weakening the bond between the core and the outer skin. When the strength of the bond becomes too weak, the foam core and skin can peal apart. The failure mode that NASA is encountering, however, is more of an explos ive and sudden blowout of foam from the ET. Since the foam loss seems to be the great est near areas where the surface of the tank is somewhat uneven (from bolts and fittings), NASA believes that voids created when the foam is first being sprayed down are the primary reason why the foam is able to break off. These are indeed defects between the ET, or skin, and the foam, or core material. To model this phenomenon, NASA devised an experiment to examine how CMOD Load, P Pmax = PQ
PAGE 39
26 various voids within foam test panels coul d contribute to large pieces of foam being blown out under certain conditi ons. Rectangular test panels like the one shown in Figure 214 are used for the divot testing. Figure 214. Divot test panel In this experiment a cylindrical bore is mach ined into the panel. At the top of that bore, a sharp notch is created to simulate a crack, or defect, near this void. Liquid nitrogen is then pumped into the bored hole and a heat flux is applied across the other side of the foam panel. Twodimensional sc hematics of this process are shown in Figures 215 and 216. The applied flux warms the liquid nitrogen and eventually a phase change takes place. The experiment takes place in a thermovacuum chamber to simulate the environment the foam is exposed to during lift off. As the Shuttle ascends to space, the atmospheric pressure is dropping, but the gaseous N2 applies pressure on the walls of the void and the crack faces. With enough force, the flaw turns and propagates toward the surface.
PAGE 40
27 Figure 215. 2D view of divot test In Figure 217, we see the aftermath of one of these experiments. A frustum shaped piece of foam has been blown out and the crater left behind resembles a divot, similar to the removal of a small piece of turf after a golf shot has been played. It is here that predicting the crack turning angle could have a practical application. Figure 216. 2D view of test pane l after heat flux has been applied For a given void and flaw size, along with the applied tractions, it might be possible to calculate this angle, denoted as c in Figure 216, by calculating neartip stresses and applying a turning theory such as the maxi mum tangential stress criterion. The concept of the crack turning angle and methods to comput e it will be covered in a later chapter.
PAGE 41
28 Figure 217. Foam blow out from divot test Summary The pertinent theory (Gibson and Ashby 1988) regarding cellular solids is presented and various formulas are availabl e depending on what t ype (open or closecelled) of foam one is analyzing. Formulae for anisotropic foams are presented as well. The elastic constants, for both isotropic a nd anisotropic constant s can be estimated through the relations presented in the early s ections of this chapter or via experiment. While there is extensive literature available on foams, NASA decided to perform tests to evaluate the elastic constants and fracture properties at various temperatures. This enabled NASA to have a large collection of experimental data for several material orientations, over a wide range of geometries. Once the tests were conducted, determin ing the fracture properties can be done using the wellknown empirical re lations for standard test specimens. One of the main assumptions behind these equations is that th ey are to be applied to a continuum. The foam is a porous material and it has been shown that cell size can impact fracture properties if one is conducting analyses with short cracks. However, continuum
PAGE 42
29 assumptions are acceptable as long as the crack is an order of magnitude longer than the cell length. Lastly, additional and highly specialized te sts were performed to investigate how the foam is able to free itself from the tank. NASA believes that voids near fittings and bolts allow liquefied air to collect in these cav ities prior to liftoff. During liftoff, the air rushing over the outside surface of the insulation warms the liquid nitrogen which converts it into a gas. The pressure, while sm all, appears to be sufficient to drive a crack toward the exterior and this seems to accoun t for the foam loss outlined in chapter one. An idealized case of this process is examined through the divot test. These experiments entail machining oval or cylindric alshaped voids into a foam test panel and pumping in liquid nitrogen. The cryogen is heated when a flux is applied to the top part of the foam panel. At the top of the void, a sharp notch is inserted and when the nitrogen is turned into a gas, the pressures exerted on the void walls and crack faces are sufficient to drive a crack toward the surf ace. It is here that a dete rmination of the crack turning angle could be of some use to NASA and ther e a few theories preval ent in the literature on this topic. Most of them require the evaluation of neartip stresses and/or strains. Since this an anisotropic material, the gene ral equations for isotropic materials are not applicable in most cases. In the next chap ter a detailed discussion covering linear, elastic, anisotropic fracture mechanics is co vered and these concepts are applied in chapter four when crack turning theory is presented.
PAGE 43
30 CHAPTER 3 AN INTRODUCTION TO FRACTURE MECHANICS Material Definitions For local crack tip calculations there are ma terial definitions, or classifications, that need to be defined before moving on with discussions on how to compute neartip stresses for materials that are anisotropic, or have directiondepe ndent properties. In general, most texts on fracture mechanics focu s on materials that are isotropic. Here the constitutive matrix contains three elastic co nstants, two of which are independent, E and For isotropic materials, there are an infini te number of planes of material symmetry and the strains and stresses are rela ted to each other via equation (31). 1000 1000 1000 2100 210 21 x x xx yy yy z z zz yz yz x z xz x y xyEEE EE E E E symE (31) Many other types of materials can be char acterized by the number of planes of internal symmetry (Dieter 1976). Cubic materi als, for example, have three independent material parameters and nine pl anes of symmetry. Materials with three internal planes of symmetry are known as orthotropic materials and many engineering metals fall into this category, such as rolled aluminum. A special class of orthotropic materials has a single plane of symmetry that is also isotropic. These are known as transversely isotropic materials and the constitutive relations for this material were defined in the previous
PAGE 44
31 chapter. Monoclinic materials possess a si ngle plane of material symmetry, but the behavior in that plane is not isotropic. Being cognizant of the constitutive relations for the particular material in question (be it isotropic or fully anisotropic) is importa nt because many of the initial derivations of neartip stress fields for an isotropic materials tried to take advantage of material symmetry to decouple inplane, and outofplane displaceme nts, which, in turn, makes the mathematical formulati on a little less complicated. Isotropic Fracture Mechanics For isotropic materials, the neartip stress fields have been analyzed by Westergaard (1939), Sneddon (1946), Williams (1957), and Irwin (1960) among others. These solutions pertain to specific cracked conf igurations with applied farfield loads and are governed by a single parameter, K, the stress intensity factor (Anderson, 1991). (,)()higher order terms 2 ijijK rf r (31) Figure 31. Coordinate sy stem used in equation (31) The coordinate system for equation (31 ) is shown in Figure 31, where the fij term is a trigonometric function. We note the Greek lette r, theta, is normally used to define the angle with the horizontal in the polar coordina te system presented above, but theta will be
PAGE 45
32 used later on to define the material orientation. A different variable name is selected to avoid confusion. The key points about (3 1) are that the neartip stresses have no dependence on the elastic constants and the first term is proportional to 1/ r and, therefore, the stresses approach infinity as r is made smaller and smaller. Also, as r approaches zero the higher order terms in (3 1) get smaller or remain at some finite value. Thus, in many texts and papers the analysis is restricted to small values of r such that the higher order terms are neglected. There is typically a roman numeral written next to the stress intensity factor to denote the mode of loading and Figure 32 displays the three modes of deformation. Figure 32. From left to right: mode I ( opening), mode II (shearing), mode III (tearing) While the stresses are proportional to 1/ r, the displacements are related to requation (32). 22 x 22 y s zcos 1+2sinsin +1+2cos0 2222 u 1r u=sin +12coscos 12sin0 22 2222 u 00sin 2 KI KII KIII (32)
PAGE 46
33 where = 34 for plane strain and = (3)/(1+ ) for plane stress and s is the shear modulus. Linear Elastic Fracture Mechan ics for Anisotropic Materials In general, the equations in the previous section cannot be applied to this problem directly because the material properties are directiondependent. The analytical work geared at developing the neartip stress fields for elastic anisotropic materials has been examined by Sih et al. (1965), Bowie and Freese (1972), Barnett and Asaro (1972), Bogy (1972), Tupholme (1974), Kuo and Bogy (1974) Rathod (1979), Hoenig (1982), and Dhondt (2002) among other researchers. The important work of Sih et al. is one of the most cited references for determining the stress fields near the crack tip for anisot ropic materials. Following Westergaards use of stress functions and complex variables to determine neartip stresses for isotropic materials, Sih et al. use a similar approach for materials with directiondependent properties. Whenever stress f unctions are utilized to solve this kind of problem, there is normally a corresponding charac teristic equation (a deta iled explanation of the characteristic equation can be found in Appe ndix A). For the most general case, the order of the characteristic equation is six. In Sih et al. s paper the material in question is monoclinic and the crack is resting in the symm etry plane. The constitutive relations for monoclinic materials are s hown in equation (33). 11121316 222326 3336 4445 55 6600 00 00 0 0xx yy zz yzyz x zxz x yxySSSS SSS SS SS S symS (33)
PAGE 47
34 This configuration is advantageous because the inplane and outofplane displacements become decoupled, which in turn reduces the order of the characteristic equation to four. 22 122121 12211221 12 12211221 12 121 ReRe0 1111 ReRe0 1 Re 2xx yy xy xz yzQQQQ QQQQ r 12 121212 3 3 3111 Re0 00Re 1 00Re KI KII QQQQ KIII Q Q (34) 31323633 zzxxyyxySSSS (35) where the elastic constants in compliance form are defined as iijjS (36) 1222112211 1212 1222112211 1212 3 4534411 ReRe0 211 ReRe0 00Rex y zpQpQpQpQ u KI r uqQqQqQqQKII u KIII Q CC (37) The roots 1, 2, and 3 for this particular formulation come from the special case for monoclinic materials listed in Appendix A. 2 111612'''jjj p SSS (38) 122622'''jjjqSSS (39) Sih et al. enforce a plane strain condition and the reduced compliance matrix becomes
PAGE 48
35 33 33 ij ijijSS SS S (310) cossiniiQ (311) The C45 and C44 terms in (37) are the elastic constants which are in stiffness form iijjC (312) where 1 ijijCS (313) The aforementioned decoupling of the anti plane displacement is now evident in equation (37). The uz displacement does not depend on KI and KII. This is similar to the isotropic case in equation (32); the z component of displacement only depends on KIII. When analyzing the foam problem for NAS A, a crack, or defect, can be at any orientation relative to the insulation. In th at sense, it does not seem very practical to employ the Sih et al. solution unless we know the crack always lies in the symmetry plane. More general solutions, ones that do not require a decoupli ng of the antiplane displacements, will be used instead. The asymptotic solution derived by Hoenig is the one we will adopt to get the neartip stress fields. Materials with virtually a ny orientation, irrespective of the location of the crack front, can be modeled with this me thod. Hoenig uses the pioneering work of Lekhnitksii (1963) to derive the stresses and displacements given below by using complex variables and stress functions to deri ve a coupled set of differential equations whose solution entails determining the roots of the ensuing sixth order characteristic
PAGE 49
36 equation. We notice that K once again arises as a scalar multiplier just as it did for the isotropic solution. 21 3 1 1 3 1 1 3 1 1 3 1 11 Re 2cossin 1 Re 2cossin 1 Re 2cossin 1 Re 2cossin 1 Re 2cossiniijj x i i ijj y i i iiijj xz i i iijj yz i i iijj xy iNK r NK r NK r NK r NK r 3 1 i (314) 313234353633 zzxxyyyzxzxySSSSSS (315) 3 1 12 Recossiniijjllj jr umNK (316) Here an important distinction is made be tween the general case and Sih et als formulation. Notice in (316) how the mode IIII SIFs are to be summed over all three modes of displacement. Thus, all three components of displacement are dependent on KIKIII. 2 11116121514 2224 2212625 4244 3414645''''' '' ''' '' '''iiii iii ii iii iimSSSSS SS mSSS SS mSSS (317)
PAGE 50
37 The i are the roots that occur in conjugate pa irs and they arise from the sixth order polynomial, equation (318). The coordinate sy stem used in equations (314) and (316) is the same as the one used in (31). 2 4230 lll (318) Where 2 2554544 32 3151456254624 432 41116126626222 222 lSSS lSSSSSS lSSSSSS The N matrix in equation (314) is defined as 123 123111ijN (320) Where i is defined as l3(i)/ l2(i). There are two situations, also called degenerate cases, where this solution is invalid, however. If antiplane shear and pl ane strain displacements becomes decoupledfor monoclinic materials, sayHoenigs form ulation falls apart. He does present an alternate solution for monoclinic materials which turns out to be consistent with the Sih et al solution mentioned earlier. His solution also encounters difficulties if the crack is lying in a plane that is isotropic. Th e Nmatrix, equation (320), becomes singular because 1 and 2 = i. However, we have neartip solutions for the isotropic cases. These are shown in equati ons (31) and (32). Determining the Stress Intensity Factor K In the above sections, it is necessary to develop solutions for determining the stress intensity factor, K. Over the past 50 years, numerous methods have been developed to
PAGE 51
38 determine this parameter. For simple geometries, a handbook (Tad a et al, 1973) of Ksolutions could be consulted. Alaibadi a nd Rooke (1991) present an excellent summary of the literature on the subject of numerical ev aluations of SIFs. They divide many of the available methods into three main categor ies, denoted as stages in Table 31. Table 31. Various methods for determining the SIF Stage 1 Stage 2 Stage 3 Handbooks Superposition Collocating, or mapping Stress concentration Integral transform Stress distribution Body force method Greens function Method of lines Approximate weight function Finite element method Compounding Boundary element Alternating technique As one might have guessed, the procedures listed in the Stage 3 category are the most accurate. Since computers have become less expensive and more powerful in recent years, many of the current studies geared to ward Ksolutions incorporate the use of the finite element (FE) method or boundary elem ent (BE) techniques. Boundary element methods are particularly attractive because the dimensionality of the numerical problem is reduced, i.e. 2D problems involve discre tizing the lineboundary of the domain and for 3D problems, just the surface of the domain is discretized. This analysis will call on the FE method to help compute K, and there are a few ways of determining it with this method, such as: extrapolation of neartip stress and displacement fields, Rices Jintegral (1968), strain ener gy, and the virtual crack extension method. FRANC3D Next Generation Software The Cornell Fracture Group has developed software capable of creating crack meshes and computing an anisotropic Ksolu tion. This particular method of solution would fall under Stage 3 in Table 31, and in this case Rices Jcontour integral along
PAGE 52
39 with displacements output from FE software will be used to determine K along the crack front. This software is very versatile in that it is not restricted to analyzing any one kind of material; Ksolutions can be generated fo r isotropic and fully anisotropic materials alike. Since many FE models were built and the subsequent Ksolutions computed with this software, it would be prudent to discuss the pertinent theo ry (BanksSills et al, 2005) that governs this particular way of calcula ting K. Rices J inte gral, a conservation integral that measures the energy flux into th e crack tip, is equivalent to the strain energy release rate (denoted as G) for small scale yielding. Furthermore, it can be shown for a plane problem that only has KI and KII that 22()GJKIKII (321) Where = (12)/E for plane strain and = 1/E for plane stress. One can see that the Jintegral equation, as is cannot give us the SIFs on an individual basis. i iu JWdyTds x (322) with W = (1/2)ijij, also known as the stra in energy density, T = ijnj and ui is the displacement vector. The path, is defined in Figure 33 Figure 33. Path of the Jintegral
PAGE 53
40 The Jintegral was initially formulated for 2D problems, however since this material has properties that ar e direction dependent, our analys is requires the use of 3D FE models. The Jintegral, then, needs to be modified so it can be used in integrations that take place over volu mes. Following Li et al (1985), th e Jintegral can be rewritten as ds x q W x u Jj j i ij 1 1 (323) In this formulation, the crack is undergoing a virtual extension, and only a small portion of the crack front is being a dvanced. When using the formulation from the Li et al paper, the Jintegral is now evaluated in a cylindrical domain centered over the crack front. Figure 34. Cylindrical do main and virtual crack advance for a 3D crack front Following Figure 34 we can define increm ental area of virtual crack advance as max q LAaqsds (324) The definition of J in equation (322) requires to be very small. It is difficult to compute stresses and strains about a vanishingl y small path about the crack tip. Instead,
PAGE 54
41 area integrals (Moran and Shih, 1987) are us ed (for 2D problems) because they are easier to implement numerically. So, the q term in (323) is simply a mathematical tool that enables us to recast the Jintegr al into a slightly more usable form. This virtual extension method is the mode rnday approach to solve computational fracture mechanics problems. One of the main advantages of such a method is the ability to compute the energy release rates across th e crack front for 3D problems. Since the energy release rates can vary over the front that implies K can also change with respect to crack front position as well q t tA J ds s q ds s q s J J ) ( ) ( ) ( (325) where qt(s) is the crack tip perturbation. Rices Jintegral does not allo w us to look at the SIFs on an individual basis for a mixed mode crack problem. This hurdle wa s first surpassed by Ishikawa (1980) where he decomposes the stress, strain, and disp lacement fields into symmetric and antisymmetric parts. This allows him to use th e Jintegral to comput e the Ks individually. Ishikawas method, though, has not been exte nded for anisotropic materials. Yau et al. (1980) uses the Mintegral developed by Freund (1978) to determine KI and KII, individually, by using the idea that two separate and inde pendent equilibrium states (denoted by superscripts (1) and (2)) for a linea r elastic, isotropic material can be related to a third equilibrium state, denoted with no superscript. That is to say the stresses, strains, displacements, and SIFs from these sepa rate equilibrium states can be superposed (1)(2)(1)(2) (1)(2)(1)(2) (1)(2)(1)(2) ijijij ijijij iiiKIKIKI KIIKIIKII uuuKIIIKIIIKIII (326)
PAGE 55
42 The Jintegral can now take the form (1)(2)(1,2)JJJM (327) The M(1,2) term is known as the Minte gral initially derived by Freund and it can be recast in terms of the aforementioned equilibrium states (2)(1) (1,2)(1,2)(2)(1,2) 1, 11 ii ijijj juu q MW x xx (328) This integral describes the interaction between the two equilibrium states, and sometimes the Mintegral is called the interacti on integral in the literature. The first step in the procedure to determine the SIFs is to define the strain energy release rate using Irwins (1957) crack closure in tegral. Consider the crack in Figure 35. Now a compressive, or closure stress, is applied in such a fashion to clamp the crack faces down along the length r. The work it takes to perform this closure can be related to the energy release rate, G. Irwin goes on to use the force displacement curve and equate that to the work required to close th e crack. Substitutions are made for neartip displacements and stresses and eventually an integral expression is obtained that relates the energy release rate to the stress intensity factor. In a similar manner, Hoenig (via Equa tion 330) provides us an analytical expression for the energy release for general an isotropy in terms of the SIFs and material constants. Irwins work was initially developed for brittl e, isotropic materials; his expression for the energy release rate includ es the isotropic Poi sson ratio and the Young and shear moduli. 0 01 lim(,0)(,) jjGrurdr (329)
PAGE 56
43 Figure 35. Definition of Irw ins crack closure integral Now, if the stress and displacement fiel ds are substituted in (329), we have 111 2131 Im()Im()Im() 2iijjiijjiijjGJKImNKKIImNKKIIImNK (330) When we use the relations for SIFs in (326), we get (1)(2)1(1)(2)1 21 (1)(2)1 3()Im()()Im() 1 2 ()Im()iijjiijj iijjKIKImNKKIIKIImNK GJ KIIIKIIImNK (331) The terms in (327) have bars over them because a certain definition of J is used in (323). Let us define an alternate form of (327) where J is not normalized with respect to the extended area as (1)(2)(1,2)JJJM (332)
PAGE 57
44 Now we substitute (331) into (332) and for clarity the individual terms of that equation are listed below (1)(1)1(1)(1)1(1)(1)1(1) 213 (2)(2)1(2)(2)1(2)(1)1(2) 2131 Im()Im()Im() 2 1 Im()Im()Im() 2iijjiijjiijj iijjiijjiijjJKImNKKIImNKKIIImNK JKImNKKIImNKKIIImNK (333) (1,2)(1)1(2)(2)1(1) 1 22 2 (1)1(2)(2)1(1) 11 (1)1(2)(2)1(2) 33ImIm ImIm ImImiijjiijj iijjiijj iijjiijjMKImNKKImNK KIImNKKIImNK KIIImNKKIIImNK (334) It is now evident how the decomposition of displacements, strains, and stresses into separate equilibrium states is helpful. One can see that we have two definitions of the Mintegral, equations (328) and (334). If they are equated, the individual SIFs (KI(1), KII(1), KIII(1)) can be computed. But before that step, we must define auxiliary solutions for the second equilibrium state, de noted by superscripts (2a, 2b, 2c). The terms with the (1) superscript are going be computed numerically via FE software. The three auxiliary solutions are define d as: 2a, 2b, and 2c. For case 2a, KI(2a) = 1 and KII(2a) = KIII(2a) = 0. For case 2b, KII(2b) = 1 and KI(2b) = KIII(2b) = 0, and for case 2c, KIII (2c) = 1 and KI(2c) = KII(2c) = 0. Since the SIFs are being prescribed beforehand, the stresses, for those equilibrium states, can be computed from (314). Now the two relations for the Mintegral are equated (1)1(2) 2 (2)1(1) (2)(1) 2 (1)(2)(1,2) (1)1(2) 1 1 11 (2)1(1) 1 (1)1(2) 3 (2)1 3Im() Im() Im() 1 2 Im() Im() Im(iijj iijj ii ijijj iijj j q iijj iijj iijKImNK KImNK uu q Wds KIImNK xxx A KIImNK KIIImNK KIIImNK (2))j (335)
PAGE 58
45 From (335) we can assemble a system of equations that contai n the three unknown SIFs we are after 11111 (1)(1,2) 2111223123 11111(1)(1 2211123213 (1) 11111 23311332332ImImIm Im2ImIm ImIm2Ima iiiiiiiiii q iiiiiiiiii iiiiiiiiiimNmNmNmNmN K IMA mNmNmNmNmNKIIM KIII mNmNmNmNmN ,2) (1,2) b q c q A M A (336) Finally, we note that (336) i nvolved the extraction of stresse s near the crack tip. These stresses are taken from the Gauss point locati ons within the elements that encompass the crack front. Also, the integral in (336) is evaluated numerically using Gauss quadrature. Computing SIFs Using Displacement Correlation Displacement correlation (DC) is perhap s the oldest and simplest method to compute SIFs via FE displacements. Chan et al (1970) were among the first group of researchers to employ such a method to obtain SI Fs this way. In general, the correlation point is selected where the crack tip displacements are the larg est so the relative error in displacements is small. Another advantag e in using such a method is the SIFs are separated into the three modes defined in Figure 32. In it s most basic form, crack tip opening displacement (CTOD) can be used to determine the SIF by equating the displacements at the points in questi on to some analytical solution. From equation (32) with = and Figure 36 for a plan e strain problem dealing with an isotropic material, we have 2 2 1ybyasuu KI r (337) 2 2 1xbxasuu KII r (338)
PAGE 59
46 where s, again, is the shear modulus determined by E/2(1+ ) and = 34 for plane strain We note that uy and ux are the opening and sliding m odes of displacement. The above method is sometimes called a onepoint matching techni que since it involves extracting displacements from one point near the crack tip. The downside of such a Figure 36. Correlation poi nts typically used for displacement correlation method is that very refined meshes are needed for accurate results. Improved accuracy is obtained through the use of a special elem ent specifically designed to capture the singularity present at the crack tip. From equation (31), we see how the stresses have a 1/ r singularity at the crack tip. For accurate computations, we would like our FE models to also have this singularity. Quarter point elements are one way to do this. The first attempts at modeling cracks were done so using quadril ateral elements with the midside nodes moved to quarter point locations (Barsoum, 1976 ) When this is done the desired strain singularity is achieved. Howeve r, one downside of using this type of element is that the singularity is only present on the edges that contain the quarter point. If triangularshaped elements are used instead, the singul arity exists both along the edges and within the element (Anderson, 1991). In Figure 37, we see how an eight node quadrilateral element can be transformed to a triangle and then cert ain nodes are moved to the quart er point locations. Let us
PAGE 60
47 derive how the required strain singularity is obtained along th e edge of the quadrilateral element with the midside nodes moved to th e quarter point locations. Consider the bottom edge of the element on the right in Fi gure 38. Nodes one, five, and two lie along this boundary. In order to show how the stra in singularity comes about, we will need to list the shape functions for this particular element. Th e shape functions (Logan, 2002) for the isoparametric element at the corner nodes are 1 111 4iiiiiN (339) where i is the number of the shape functions at th e corresponding node (i.e. N1 is the shape function for node 1) and 1,1,1,1 (1,2,3,4) 1,1,1,1 (1,2,3,4)i ii i (340) We also have shape functions for the midside nodes 2 21 11 1,1 (5,7) 2 1 11 1,1 (6,8) 2iii iiiNti Ni (341) The global origin is placed at the corner of the element. We will n eed the shape functions at nodes one, two, and five. Setting = 1 we have 1 2 2 51 1 2 1 1 2 1 N N N (342)
PAGE 61
48 Figure 37. Triangular quarter point elements Figure 38. Isoparametric element and dege nerated element with midside nodes moved to quarter point locations The global x coordinate is related to th e isoparametric coordinate system via ii x Nx (343) where the repeated index is an implied summation. Setting x1 = 0, x2 = L, and x5 = L/4 we have 21 11 24 L xL (344)
PAGE 62
49 We can now solve for 21 x L (345) In the same way that we can relate the globa l x coordinate to the shape functions, we can also determine the displacements using the shape functions iiuNu (346) Again determining the displacements along th e edge containing nodes one, two, and five we write 2 12511 111 22 uuuu (347) Now we can substitute (345) into (347) to obtain 1251 12221224 2 xxxxxx uuuu LLLLLL (348) We note that the displacement shown in (348) has a x term present. This is consistent with equation (32) where the displacements vary with r. We need to differentiate equation (348) to obtain the strain in the x direction 12513411424 22xu uuu xLLL xLxLxL (349) Finally, it is shown in (349) how the strain disp lays the desired singularity (x1/2). One can see why the quarter point element ha s been the standard for modeling cracks for the past 30 years. No special programming, or internal tampering, with the FE code is required. Any FE program that carries quadr atic, or higher order, elements can support these special ones with the nodes moved to the quarter points.
PAGE 63
50 Let us now return to our discussion of using displacement correlation to extract SIFs from FE models. If quarter point elemen ts are used, equations (337) and (338) are modified because r = L/4 4 2 1ybyasuu KI L (350) 4 2 1xbxasuu KII L (351) The SIFs extraction via the DC method need not be restricted to just one set of nodes, or correlation points. One can also consider th e displacements over the whole element. The form of these equations depends on the shape functions of the elements. Ingraffea and Manu (1980) derive a DC extraction method usi ng 3D elements. This formulation is particularly useful to the present study si nce 3D simulations are used to model the fracture response of the foam material. Ingraffea and Manu use a special labeling convention shown in Figure 310 for a collapse d 20 node brick element (Figure 39), also known as a 15 node wedge. Figure 39. 20 node brick and 15 node wedge elements The shape functions for a 20 node brick el ement are listed below. The local ( , ) coordinate system is placed at the center of the brick element. For the corner nodes at positions i = 1,2 and i, i, i = 1
PAGE 64
51 111 2 8iii iiiiN (352) For the nodes at positions i = 17, 18, 19, 20 and i = 0, i = i = 2111 4ii iN (353) At i = 10, 12, 14, 16 and i = i = 0, i = 2111 4ii iN (354) Finally, for the nodes at positions i = 9, 11, 13, 15 and i = i = i = 0 2111 4ii iN (355) For anisotropic materials, Saouma and Siki otis (1986) extend Ingraffeas work for anisotropic materials. They use the displacement relations developed by Sih et al A similar procedure is followed in this study, but instead of using displacement equations from Sih et als paper, the relations in (31 6) are to be employed. The equations to determine the anisotropic Ks, using Hoenig s displacement equations, are as follows 12 KI KIIBA L KIII (358) where [B] is 1 1 1 2 1 3cossin []Recossin cossinjjlj jjlj jjljmN BmN mN (359)
PAGE 65
52 Figure 310. Collapsed 20node brick element ''''''''' 2 ''' ''''''''' 2 '''1 22224444 2 1 22 2 1 22224444 2 1 22 2 2BCEFBCEFDBCEFBcEF FCDFCD BCEFBCEFDBCEFBcEF FCDFCD Bvvvvvvvvvvvvvvvvv vvvvvv uuuuuuuuuuuuuuuuu uuuuuu w ''''''''' 2 '''1 2224444 2 1 22 2CEFBCEFDBCEFBcEF FCDFCDwwwwwwwwwwwwwwww wwwwww (357) Once again, the repeated indices are an implied summation. The terms in the B matrix only depend on the roots of the sixt h order characteristic equation. Effects of Temperature on the SIF Solution In chapter two it is pointed out that temp erature plays no small role in causing the foam to break off from the ET. For example, when Shuttle is resting on the launching
PAGE 66
53 pad, one side of the foam is exposed to the temperature of li quid hydrogen and oxygen; roughly 300oF. The other side of the foam is exposed to air, perhaps at 75 oF prior to liftoff. Over typical acreage sp rays, the foam is approximately three inches thick. This small thickness coupled with a large thermal gradient suggests that thermal strains could be significant. As shown with the divot tes ting in chapter two, even though there are no applied farfield mechanical loads, the therma l gradient seems to be sufficient enough to impart stresses capable of driving a crack to ward the surface of the ET resulting in a loss of material. To account for thermal strains, we simply a dd this term to our constitutive matrix. Following Hyer (1998), the stressstra in relations can be written as 1111 11 111213 22 2222 212223 33 313233 3333 44 23 23 55 13 13 66 12 12(,) 000 (,) 000 000 (,) 00 0Th ref Th ref Th refTT SSS TT SSS SSS TT S S symS (360) where the total strain is { 11 22 33 23 13 12}T, Th 11(T, Tref) is the free thermal strain in the 11direction, and Tref is the reference temperature. The mechanical strains are 111111 22 2222 33 3333 23 23 13 13 1212(,) (,) (,)mechTh ref mech Th ref mech Th ref mech mech mechTT TT TT (361) or, stated more concisely the total strain is total = mech + thermal (362)
PAGE 67
54 If the thermal expansion is linear with the te mperature change, we can rewrite (360) as 11111 111213 22222 212223 33333 313233 44 23 23 55 13 13 66 12 12000 000 000 00 0 T SSS T SSS T SSS S S symS (363) The i term is the coefficient of thermal expansion and T is the difference of temperature between the state of inte rest and the reference state. The above equations allow us to compute th ermal strains in the numerical models. Incorporating thermal and mechanical loads will involve two sepa rate steps because a thermal gradient implies that a temperature distribution must be obt ained before solving the model that contains both the thermal a nd mechanical boundary conditions. The first step is to solve the conduction problem usi ng the thermal conductivity coefficient listed in chapter two. The nodal temperature distributi on is saved for later use. In a subsequent run, now with mechanical loads, the nodal temperatures are carried over as body forces along with any applied mechanical loads. Summary Summing up, the analytical expressions from Hoenigs neartip solution is presented. His equations depend solely on the roots of the characteri stic equation and the stress intensity factor, K Hoenigs solution is not as well known as Sih et als paper and as such the vast majority of studies dealing with neartip solutions for cracked anisotropic bodies use the Sih et al formulae to compute the stresses. This is acceptable as long as one realizes those particular equations are valid for cert ain configurations where a decoupling of inplane and outo fplane displacements takes pla ce. Hoenigs solution is
PAGE 68
55 more general and his formulation of the energy release rate is the one chosen to compute the SIFs within the FRANC3D software. The DC technique presented also uses Hoenigs equations, but only the ones that pertain to displacements. These two methods, interaction integral and DC, are presented to show how the mode IIII SIFs can be computed. These fo rmulations are different in a few ways. Computing K with the Jintegral is inherent ly more accurate because only a moderate level of mesh refinement is needed for th e solution. However, the implementation of such a method is very involved. Alternativel y, the DC approach is much easier to apply but the accuracy of the solution does depend on the level of mesh refinement. Finally, the relations dealing with handling thermal strains are presented. These equations assume an expansion that is linear with temperature change. Handling thermal loads is straightforward by first solving the conduction problem to obtain the nodal temperature distribution. With that, those values are forwarded as body forces that are applied alongside the mechanical loads.
PAGE 69
56 CHAPTER 4 CRACK TURNING THEORY AND FINITE ELEMENT MODELING Crack Turning Theory There are three theories preval ent in the literature to pr edict incipient crack turning angles for isotropic materials: maximum energy release rate, minimum strain energy density, and maximum hoop stre ss. All of these theories are based on LEFM concepts and were initially deve loped from experiments involving brittle plates. The turning theory proposed by Hussain et al. (1974) seeks the direction where the energy release rate, G will be a maximum. The minimum strain energy density theory (Sih, 1974) postulates crack turni ng in the direction where the st rain energy density is at a minimum. The maximum hoop stress theo ry proposed by Erdogan and Sih (1963) predicts crack propagation in the direction where, defined in equation (41), is a maximum relative to the crack tip. Crack turn ing, in this study, denotes the incipient turning angle. We define the crack turning angle, c, below in Figure 41. For isotropic materials, the ne artip stresses (in the r, polar coordinate system) have been derived by Williams (1957). If the derivative with respect to theta of equation (41) is taken and set to zero, the critical angle, c, can be determined; equation (44). 213 coscossin 222 2 KIKII r (41) 213 cos1sinsin2tan 2222 2rKIKIIKII r (42)
PAGE 70
57 Figure 41. Definition of crack turning angle 1 cossin3cos1 2 22rKIKII r (43) 2 1118() 2tan 4()cKIIKI KIIKI (44) The maximum hoop stress theory is the easies t to implement and th at is perhaps the main reason why it is widely used for turni ng angle predictions. The minimum strain energy density criterion is also very popular and there is some debate (Gdoutos, 1984) as to which one is superior. Some might argue that it is not correct to use just a single component of stress to predict the incipient turning angle, whereas a quantity such as strain energy density seems to be more co mprehensive as all components of stress and strain are included. Ma ccagno and Knott (1989, 1992) study mixedmode facture behavior of PMMA and lightly tempered steel alloys and in both cases the crack turning angles are well approximated by the maximum tangential stress theory. There have been some topical investigations for crack tu rning within materials that have direction dependent prope rties. Buzcek and Herakovi ch (1985) pred ict the crack
PAGE 71
58 extension angle for orthotropic ma terials and some of the recent work in this area mirrors their ideas. They assume the tensile strength, T of the orthotropic material varies with direction, Since it is difficult to measure the strength of a material in all possible orientations, they assume that T = f( , XT, YT) where is the angle in a polar coordinate system, is the material orientation and XT and YT are the strengths of the material in the axial and fiber directions respectively (they are analyzing orthotropic composites). The equation for T, also denoted as a fracture resistance parameter, must be independent of if the material is isotropic. The lo west order function th at satisfies these requirements is an ellipse. This f unction is plotted in Figure 42 in (T, ) polar coordinates. Now T can be expressed as 22sincosTTTXY (45) Figure 42. Elliptical relationship for T The crack will turn, according to Buzek a nd Herakovich, in the direction where the ratio of the hoop stress to T is a maximu m. The hoop stress, sometimes called the tangential stress (or ), is defined in Figure 43. Consider a point P and a vector
PAGE 72
59 connecting that point to the crack tip. Th e hoop stress is normal to this vector that connects P to the crack tip. Figure 43. Definition of hoop stress Mixed mode crack propagation in anisotr opic materials is also investigated by Saouma et al. (1987). The neartip stre sses are obtained via Sih et al. (1965) and in their analysis the maximum tangential stress theory is used to predict the angle of propagation. Saouma et al. use a slightly different criterion for cr ack turning, but it is of similar flavor to what Buzek and Herakovich propose. Instead of using strength as a fracture resistance parameter, the plane strain fracture toughness is utilized. Boone et al. 1987) use FE models to study cr ack propagation in orthotropic materials. The idea of maximizing the ratio of the hoop stress to some strength parameter is once again utilized to pred ict the crack turning angle. Chen (1999) analyzes crack propagation in aircraft fuselage s. He, too, makes use of a ratio of the hoop stress to the plane strain fracture toughness. In this case, the material in question has isotropic elastic moduli, but the fracture properties are dire ction dependent. A similar analysis is conducted by Pettit (2000). He examines crack turning in rolled aluminum, which again, has isotropic stiffness properties but an isotropic fracture characteristics.
PAGE 73
60 Carloni et al. (2003) and Nobile and Carloni ( 2005) consider incipient crack turning for an orthotropic pl ate under biaxial loading. Th e former researchers utilize essentially the same fracture resistance pa rameter as Buzek and Herakovich, whereas the latter study uses a slightly m odified version of the Saouma et al. turning criterion. Carloni et al. and Nobile and Carloni also call on complex variables and stress functions in their formulations for neartip stresses. Both studies, however, make assumptions that reduce the order of the characte ristic equation that arises in such derivations (see chapter three), down to order four. While Pettits analysis is for a material with isotropic sti ffness properties, he generalizes the fracture resistance parameter to three dimensions. His analysis accounts for an arbitrary location of the crack with respect to the mate rial. Thus, we would like to use the method to obtain stress inte nsity factors given by BanksSills et al. (2005) in conjunction with Pettits formul ation for a fracture resistance parameter. This approach is a bit more generalized than the aforementioned studies because their work seems to be confined to specific configurations that us ually decouple the inp lane and antiplane displacements. Pettits method for determining the fracture resistance is formulat ed for orthotropic materials which possess cubic symmetry, i.e. three orthogonal planes of symmetry and within each plane lay two principal toughness values resulting in six principal toughness values: K12, K21, K23, K32, K13, and K31. The first number in that nomenclature is the vector normal to the crack plane and the s econd number corresponds to the direction of propagation. It is difficult to obtain a toughness for all possibl e orientations for materials that exhibit some form direction dependence. Essentially we are interpolating between
PAGE 74
61 the six principal toughness values for the orthotro pic material. This gives us an estimate of the toughness in any direc tion we desire. Consider Fi gure 44 and the material orientation as shown. What we see in that figure is block of foam with various M(T) specimens oriented within it. If the M(T) specimens are aligned with the principal axes of the large block of foam, th ere are six possible orientati ons and with each orientation there is a corresponding toughness. Pettits interpolation equation, then, requi res six unique toughness values for the orientations shown in Figure 44. In our case, only thr ee toughness values are needed, however. In Figure 44, consider th e M(T) specimens that have the K31 and K32 labels. Here the load is applied in the 3directi on and the crack is propa gating in the 1 and 2directions, respectively. Since the 3direction is the rise di rection and the 12 plane is a plane of isotropy, it is re asonable to assume that K31 = K32. That is to say, the crack is resting in the plane of the knit line for thos e two cases. For the other four cases, the fracture plane is normal to the knit line plane. Thus, K21 = K12 and K13 = K23. So the total number of required toughness values ha s been reduced to three and NASA did conduct tests, at room temp erature, to obtain them: K32 = 22.4 psi in, K12 = 17.4 psi in, and K23 = 19.5 psi in. Pettits equation that will determine the e ffective K needs severa l variables besides the six principal toughness values. Remember the goal here is to try and predict the direction of propagation. We can convenientl y describe this direction with a vector, denoted as a. That vector lies in a plane tangent to the fractu re surface, whose normal vector will be defined as n. These two vectors are shown in Figure 45.
PAGE 75
62 Figure 44. Orthot ropic toughness values Figure 45. Definition of the a and n vectors To develop the equation for the effective K, we first define the trace angles of the vector a makes with the principal planes. These are 3 12 123 231tan tan tan a aa aaa (46) The subscript in the theta term denotes the axis normal to the principal plane. Using equation (45) we can write a 2D in terpolation function of the form
PAGE 76
63 22cossinkkkikkjkKKK (47) We also observe the following trigonometric identities 2 21 22 2 21 22costan sintan bc cbc bb cbc (48) and the components of unit vector a must satisfy the relation 222 1231 aaa (49) Pettit goes on to define equa tions that can be interpreted as the fracture resistance components of a in the principal planes. Consider Figure 44 once more where the 1axis is the loading direction. To estimate an effective K on this plane, we will interpolate between two toughness values, K12 and K13. Using (47) through (49) we can write K1 as 22 1122133 2 11 () 1 KaKaKa a (410) In a similar fashion we can determine K2 and K3 22 2233211 2 2 22 3311322 2 31 () 1 1 () 1 KaKaKa a KaKaKa a (411) Finally, Pettit sums the relations in (410) th rough (411). The effective K, denoted as Kp, can now be written as 2 22222 1 112233122133 2 1 2 2 2222 3 2 233211311322 22 23(a,n) 1 11pn KKnKnKnKaKa a n n K aKaKaKa aa (412)
PAGE 77
64 FE Modeling of the M(T) Specimen While the importance of crack turning was made evident by th e divot test, this experiment is highly specialized and is not recognized as a standard fracture test specimen. A more conventional way of anal yzing crack turning can be performed with the middle tension, M(T), test specimen. Material orientation, mode mixity, and boundary conditions can all be carefully controlled with the M(T) specimen. The idea behind creating FE models of M(T) specimens is twofold. One of the objectives of this study is to examine the e ffect of mode mixity on the Ksolution, or stated differently, as the material orientat ion is changed how does this impact KIKIII along the crack front? Also, several foam M(T) specimens were fabricated and the resulting crack turning angles measured. Nu merical predictions, ba sed on Ks extracted from the FE models are used to make pr edictions via the maximum hoop stress theory discussed earlier. M(T) FE models are built using ANSYS FE software. Crack meshing is handled by FRANC3D, a FE software package devel oped by the Cornell Fracture Group. Figure 46 shows a schematic of the M(T) specime n and Figure 47 displays one of the FE models used in this analysis. The FE mode l has a unit applied pre ssure at the top and bottom faces and is simply supported. The model is comprised of entirely of SOLID92 and SOLID95 elements (ANSYS, 1999) which are used to make 3D models and have the capability to handle anisotr opic material properties. The model is twelve inches high, fiveandahalf inches wide, and an inchandha lf thick. The boxed in area in Figure 47 highlights the densely meshed region containing the crack.
PAGE 78
65 Figure 46. Schematic of M(T) specimen Figure 47. M(T) m odel created in ANSYS Since the foam has direction dependent prope rties, this needs to be accounted for in the FE models. This is best done by using di rection cosines. Using Figure 48, let us look at simple transformation of axes to define a new material orient ation. Consider two sets of axes: x, y, z and x, y, z. Initially the both sets of axes are aligned and then a rotation, anticlockwise (of angle zeta) about the x axis, takes place. Table 41 lists the 12 5.5 1.5
PAGE 79
66 direction cosines, i.e. the co sines of the angles between the primed and unprimed axes. Here 1 is the cosine of the angl e between the x and x axis, 2 is the cosine of the angle between the y and x axis, and so on. We w ill examine several material orientations with different initial crack inclinations. In this way, mode mixity is introduced into the problem and we also examine the effect of material orientation on the predicted propagation path. Figure 48. Sample transformation Table 41. Direction cosines x y z x 1 1 1 y 2 2 2 z 3 3 3 As mentioned in an earlier section, for typi cal acreage sprays there is little offset between the material and substrate coordinate systems except when the foam is applied near fittings and/or bolts where the foam ca n to rise at angles as high as 30 degrees relative to the ET. Twelve M(T) fracture te st samples with varyi ng material orientation and crack inclination were fabr icated and tested. A broken M(T) specimen is shown in Figure 49 and a comparison between the numer ical prediction and the measured turning
PAGE 80
67 angle, as well as how the foam M(T) specimens are fabricated, will be made in chapter five. Figure 49. Fractured M(T) specimen: the dotte d line is added to show the location of the knit line
PAGE 81
68 CHAPTER 5 RESULTS AND DISCUSSION Effect of Material Orientation on the SIF Solution On most sections of the ET, the foam is sprayed on in such a fashion that it rises normal to the surface. In that sense, there is very little offset between the tank (substrate) coordinate system and the foam (material) coordinate system. Also, when modeling these specimens, we should note that the curvatur e of the ET is rather large; some 28 ft in diameter and our specimens tend to include sm all defects (usually two inches). Thus, it seems reasonable to assume that curvature e ffects are minimal. But when the spraying process takes place near a bolt or fitting point the foam must be applied around, and over top of, these parts. In this situation it is possible to have the foam rising at different angles relative to the surface. The anisotropic nature of the foam require s 3D FE models to be employed so the SIFs can be evaluated along the crack front. To that end, the effect of material orientation on the Ksolution will be examined within models that have both straight and diagonal cracks. Before moving on with this discussion, let us define a few important terms that we will use extensively throughout th is chapter. The crack inclination angle, is the angle the crack makes with the horizon tal as shown in Figure 51. M(T) models use thoughcracks, or flaws that are placed through the thickness of the specimen, or model, in question. The crack front distance defined as the length of the crack in the thickness direction, shown in Figure 52.
PAGE 82
69 Figure 51. Definition of crack inclination angle To analyze the effect of material orie ntation of the SIF solution, 17 material orientation cases are defined; see Figures 54 and 55. M(T) models are constructed using ANSYS FE software and the flaw is in serted using FRANC3D Next Generation. Figure 52. Definition of crack front distance
PAGE 83
70 A unit pressure is applied to the simply s upported M(T) models and in all cases, the crack length is two inches. A picture of the ANSYS model (with dimensions) is shown below. The material orientat ions are also input into the FE models along with the loads 12 5.5 1.5 12 5.5 1.5 Figure 53. M(T) model built in ANSYS Figure 54. Definition of material orientation z y x 3 2 1 445o 8 5 6 7 10 11 12 13 9 14 15 16 30o
PAGE 84
71 Figure 55. Top view of the cones shown in Figure 54 and boundary conditions. There are many differ ent orientations at which the foam can rise relative to the tank. But when the foam is sprayed down, at some locations, the rise direction relative to the tank can be as high as thirty degrees. An example of this happening is shown in Figure 55. Figure 56. Hypothetical material orientation relative to the ET One can see how the foam can rise at angl e when it sprayed out the bolt shown in the above in Figure 56. As such we would like to have a reasonable t olerance that the rise 9 13 7 6 8 5 4 2 1 3 0 10 11 12 14 15 16 45
PAGE 85
72 direction can fall into. This w ould give rise to a range, or set, of possible orientations that could occur when the foam is being applied. Si nce the foam can rise at thirty degrees, it might be interesting to examine the possibility of the foam rising at another, larger, angle; 45 degrees, for example. Thus we now have two possible cases to examine; a rise angle that is typical near fitting points and bolts, and one that is slightly larger. The material orientations used in this study can be visualized by setting up two concentric cones. The points on the rims of the cones denote where the zaxis (local rise direction) will be for each material case. Let us consider an example. Point one in Figure 54 lies on the inner, 30 degree, cone. To define the material axes for this case, a single transformation of 30 degr ees, denoted by zeta, in Figure 57 is performed. This can be likened to the knit lines being orient ed within the M(T) specimen shown in Figure 57. Table 51 lists the dire ction cosines, i.e. the cosine s of the angles, between the primed and unprimed axes. Here 1 is the cosine of the angl e between the x and x axis, 2 is the cosine of the angle betw een the y and x axis, and so on. Figure 57. Definition of case one mate rial orientation for the M(T) FE model
PAGE 86
73 Table 51. Direction cosines for the case one material orientation x y z x 1 = 1 1 = 0 1 = 0 y 2 = 0 2 = 0.866 2 = 0.5 z 3 = 0 3 = 0.5 3 = 0.866 Figure 58 is a plot of the mode I SIF against normalized crack front distance for cases 14. The mode I SIF for the case zero orientation, the case where the material and substrate axes are coincident, is also plotted to make a compar ison. Figure 58 pertains to a straight crack that is para llel to the horizontal (no in clination). Even though this material is anisotropic, for this geometry, KII and KIII remain small compared to KI. This is not always the case, however, particular ly when the crack is inclined. To examine the effect of mode mixity, th e crack in the M(T) model is inclined by 30 degrees. All cases are rerun and in Figures 59 through 511, the mode IIII SIF is plotted versus the normalized crack front distance for cases one through four. As expected, when the crack is inclined KI decreases, a nd KIIKIII increase. This is consistent with isotropic assumptions. For the configur ation shown in Figure 61, th e KI and KII are related to or the crack inclination angle (Anderson, 1991) 2cos sincos KIa KIIa (51) We can see that KII is zero when = 0 and that KI will decrease for increasing values of
PAGE 87
74 Figure 58. Mode I SIF for cases 14; 0 degree crack inclination Figure 59. Mode I SIF for cases 14; 30 degree cr ack inclination
PAGE 88
75 Figure 510. Mode II SIF for cases 14; 30 degree crack inclination Figure 511. Mode III SIF for cases 14; 30 degree crack inclination
PAGE 89
76 The mode I SIF for = 0 is listed in tables 52 and 53. Here we can see the variation of the mode I SIF as the material orientation is changed. For these cases, the ones that resulted in the largest KI are: 11 and 15 (2.27 psi in). When you compare that number to the KI for case zero, the orientation for t ypical acreage sprays, the largest KI is just 14% higher. This is an interes ting result in that even if we l ean the zaxis by quite a bit, the resulting mode I SIF is not much great er than the mode I SIF for the case 0 orientation. However, it is po ssible to orient this material in an infinite number of ways and it could be possible to obtain a KI value fo r an arbitrary orientation that is greater than the largest value determined in this section (2.27 psi in). Table 52. KI for cases 08, 0 degree inclination Case 0 1 2 3 4 5 6 7 8 KI (psi in) 1.99 2.01 2.14 2.252.142.012.142.25 2.14 Table 53. KI for cases 916, 0 degree crack inclination Case 9 10 11 12 13 14 15 16 KI (p si in ) 2.07 2.192.272.192.072.192.272.19 Crack Turning Predictions We can numerically predict the crack turn ing angles using the maximum tangential stress criteria since we have SIFs from FE models and the neartip stress solutions from Hoenig (1982). As discussed in chapter four we will use a ratio of the tangential stress, to Kp for turning angle estimations. Let us denote the ratio of /Kp as R and where R is a maximum, that is the directi on of predicted propagation. To make a comparison, NASA has donated several broken M(T) specimens with various material orientations. The turning angles within these specimens will be measured and compared to their numerical counterparts.
PAGE 90
77 M(T) Specimen Fabrication and Determ ination of Local Knit Line Axes Fabrication of the M(T) specimens involve s spraying the foam on a metal plate, it rises and eventually a rind forms. When the foam has settled, the next layer is applied until the desired thickness is achie ved. From here, the foam is cut, or machined, from this parent block so the knit lines are running in the desired directions. One of the foam samples is shown in Figures 512. Three M(T) foam test specimens are analyzed in this section. These samples have various orientations a nd crack inclinations ( values). The first step in this analysis is to determine what the material orientation is for each specimen and we do this by measuring Figure 512. From left to right: front, left, right and rear sides of the M(T) specimen the angles of the knit lines on the specimens surface. Consider Figure 513. We need to determine the orientation of the x, y and z vectors in order to define the material orientation. This information is used to define the element c oordinate system in ANSYS. The intersection of two planes is determ ined by taking the cr ossproduct of their normal vectors. If one wants the trace angl e on a particular face of the M(T) specimen, the dot product can be used. In our case, th e normal to the knit line plane is not known beforehand. Instead we will measure two angles, and and use them to define the knit
PAGE 91
78 line plane shown in Figure 513. What we are interested in, is the trace, or intersection, the knit line plane makes on the x and y = 0 faces. Along the x = 0 face, we need the angle of the knit line measured from the horiz ontal, or yaxis. In this case the knit lines on the x = 0 face are relatively straight; is approximately zero. On the y = 0 face the angle, is also measured relative to the horizontal, or the xaxis. Thus, the traces on the y and x = 0 faces define the x and y material axes respectively. Using the cross product, the axis normal to knit line plane is determined, i.e. z = x y. Figure 513. Determination of knit line plane Local Crack Tip Computations For crack tip computations i nvolving stresses, there are a few subtleties that should be addressed in this section because the materi al classifications discussed in chapter three are in a global sense. When looking at a foam M(T) specimen, one would say it is indeed a transversely isotropic material. However, the postprocessing of neartip
PAGE 92
79 stresses requires material properties rotated into a local crack front coordinate system. Let us consider an example. First, assume we have the case one orientation; a material rotation ( = 30o) takes place about the global x axis shown in Figure 513. If the inclined crack is oriented in such a fashion so that it rests in the plane of symmetry, the outofplane and inplane displacements become decoupled and this is a socalled degenerate case described by Hoenig. His fo rmulation cannot be applied for this special case. As described in chapter three, when this situation arises (crack lying in a symmetry plane) the usual isotropic neartip formulas are utilized. Figure 514. M(T) model with inclined crack Now set = 15 degrees. The crack is no longer resting in the symmetry plane. Let us denote this as a general case and a rotati on is needed to set the properties in the local crack front system. The way the constitutiv e matrices are transformed depends on how they are defined. If they are in compliance form
PAGE 93
80 S (52) then the transformation formulas (Ti ng, 1996) for the elastic constants are 'T SSSQSQ (52) where S is the matrix of elastic constants in the rotated system and 222 111111111 222 222222222 222 333333333 232323232323232223 313131313131313131 121212121212121212222 222 222SQ (53) The order of the terms in and (52) and (5 3) is based on the convention commonly used by the composites and elasticity community wher e the stresses and strains are ordered as { x, y, z, yz, xz, xy}T. If the constitutive matrix is in stiffness form C (54) there is an alternate definition of the [QS] matrix 222 111111111 222 222222222 222 333333333 232323232323232223 313131313131313131 121212121212121212222 222 222CQ (55) The transformed elastic constants in sti ffness form are therefore determined by '[][][]T CCCQCQ (56) The two [Q] matrices are related 1T SCQQ (57)
PAGE 94
81 Equations (52) or (56) allow us to transform the material properties from one Cartesian coordinate system to another. The only info rmation we need for such a transformation is the direction cosines of the angles between th e coordinate axes where we start, or the system where the material properties are initia lly defined, and the system where we want to set the properties. The incipient turning angles analyzed in this program use the maximum tangential stress theory. In order to apply it, the first thing we need is the hoop, or component shown in Figure 515. The neartip stress es derived by Hoenig are in a (x, y, z) Cartesian coordinate system. To obtain the hoop stress, we transform axes, in this case a rotation anticlockwise about the zaxis by degrees. This is done using standard transformation equations (Boresi et al 1993). Figure 515. Cartesian and cylindrical stresses Calculating is done via 222 222222222222 x xyyzzyzzxxy (54) Here the 2, 2, 2 terms are the direction cosines of the angles between the old and new axes
PAGE 95
82 Now we can calculate the hoop stress at any r and of our choosing. The next step is to divide the hoop stress by the effective K, Kp. We do this because not only are the elastic constants direction dependent, but the fr acture properties are as well. The crack is predicted to turn in the direction where R is a maximum. The crack turning angle, defined as c is shown in Figure 516. Here we have a closeup view of a M(T) specimen containing an inclined crack. The inclined crack offset by angle, and the direction of incipient turning is denoted by the dashed line. Figure 516. Crack turning angle Let us pause briefly and quickly summarize the procedure needed to compute neartip hoop stresses and the pr edicted turning angle Enter material orientation (denoted by cas e number in Figure 52) into a FE model coordinate system
PAGE 96
83 Run the model and extract the mode IIII SIFs via FRANC3D Rotate material properties into a local crack front coordinate system Compute neartip stresses using Hoenigs formulas Transform yy to via the transformation equati ons found in Boresi et al Divide yy by Kp and determine where this ratio is a maximum to predict the turning angle, c Table 54. Measured and predicte d turning angles (in degrees) c (measured) c (predicted) c (isotropic prediction) specimen A 27 32 43 specimen B 9.7 21 27 specimen C 12 19 27 For the three test specimens, one of them had a crack inclination of 30 degrees (specimen A). The other two specimens (B a nd C) had crack inclin ations of 10 degrees. All three specimens have different material or ientations. Table 54 lists the measured and predicted incipient turn ing angles. The SIFs can vary al ong the crack front and since the free surface can influence the numerical Ksolu tion, the measured and predicted angles are computed at the midplane. For comparis ons sake, isotropic predictions for the same M(T) geometries are also tabulated. The isotropic angle is the critical angle, c, from equation (44). The experimental tests were carried out at 75oF and as such, the elastic properties at that temperature are us ed to predict the turning angle. The turning angle was measur ed using a Brown and Shar pe coordinate measuring machine. These devices are used to inspect various geometries such as cylinder heads or gear boxes. This instrument uses a probe that physically touches the work piece and the points are recorded digitally. In this case, we wish to use the points to define the fracture planes. The angle between these planes is th e crack turning angle. Within each specimen
PAGE 97
84 are four fracture surfaces. The fracture angle was measured at a location near the middle and an average of the four values is taken and listed in Table 54. Specimen A had the best agreement with experiment whereas the other two predictions fall somewhere between the isot ropic prediction and actual values. Even though both the direction dependence of both th e fracture toughness a nd neartip stresses are accounted for in these predictions, the results presented here still vary by large amount. One possible source for the difference is how the material orientation is defined and modeled. When tasked to predict the turning angle, NASA fu rnished us with the three test specimens but the orientations were not specified beforehand. These particular specimens were initially fabri cated to analyze mode mixity, not the effects of material orientation. As such, the or ientations of the knit lines vary throughout the specimens. When defining the knit line plane, the knit li nes (denoted with black dots in Figure 512) closest to the crack were used to define the orient ation, even though the material orientation can vary over the specimen. Based on the results seen here, it seems that perhaps the crack is seeking a pure mode I condition, or that the crack is seeki ng a direction that is normal to the highest principal stress. In such cases, if a crack, for example, is orientated at 30 degrees relative to the horizontal, the crack will turn 30 degr ees to it is running straight across or in a direction normal to the maximum principal stress. However, only three specimens are analyzed in this section and to some extent, the material orientation can be more carefully controlled. Obviously more samples would have to be an alyzed to determine if the cracks favor a direction that is norm al to the highest principal stress.
PAGE 98
85 Influence of Thermal Loads on the SIF Solution FRANC3D is very versatile in that it can generate different types of crack meshes and also compute a Ksolution for both isotropi c and fully anisotropic materials but at the time of this writing, it is unable to handle proble ms with a thermal gradient. This is main reason why the DC technique outlined in chapter three is utilized. To examine the effect of a thermal gradient throughout the M(T) specimen we will examine the case two orientation with a straight and inclined crack. This varies the amount of mode mixity and two different sets of temperatures are appl ied to the top and bottom faces of the M(T) specimen shown in Figure 517. In the figures that follow the primary ai m is to see how the SIFs change when temperatures are applied to the M(T) speci men. Since the Ks are evaluated using a different method, it is useful to compare th e results of these two approaches when the same set of boundary conditions is applied. We have results tabulated for several cases at room temperature. Selecting one of them, case 1 for example, we can compare the SIF results for a straight crack. The midpl ane KI at room temperature is 2.01 psi in. Using the DC method we obtain a value of 1.8 psi in, a difference of appr oximately 10%. The quality and density of the mesh does impact the alternate solution presented here. In BanksSills et al. (2005) study a similar routine is used to check the Mintegral computations. Recall that the interaction, or M, integral is used to compute the anisotropic Ksolution in FR ANC3D. Excellent agreement is shown between their various schemes, but much denser meshes are utilized within their FE models. In this study it is possible to refine th e mesh, perhaps by changing the location of the quarter point elements by some degree, but this may not be feasible because there is a
PAGE 99
86 node and element limitation within this FE software. BanksSills et al. use 2D models and that allows one to change the mesh de nsity extensively without adding a lot of computational cost in terms of phys ical model size or run time. Consider Figure 517 with = 0 degrees (straight crack). Here T1 = 0oF and T2 = 100oF. The same unit pressure is also applied. All SIF data from FRANC3D is at room temperature (75oF). The DC method used to compute th e SIFs is able to determine these parameters at any temperature. To examine the influence of thermal loads on the Ksolution, the mode IIII SIF for case two at both room temperature and with a thermal gradient is plotted in Figure 518. From Table 52, the room temperature midplane value for KI is approximately 2.14 psi in. With the prescribed thermal boundary conditions, the midplane KI drops to 1.44 psi in; a 33% reduction. Also, in a few of the figures, some odd ac tivity takes place near the free surface. When looking at Figure 518, and the line that pertains the to KI, for the DC solution, the third point drops down to approximately 0.946 psi in. This happens again in the rest of the figures that pertain to th e DC solution. One possible ex planation is that the free surface along with the boundary conditions is influencing the displacement results at that point. Sometimes quadratic elements are prone to spurious oscillations. Within Dhondts (2002) paper on computing mixedmode anis otropic SIFs, he pe rforms a smoothing operation via an averaging technique to remove the unwanted oscillations. Likewise Man et al (1993) encounter a si milar problem when using boundary elements to model a contact problem. Their solu tion involves removing the hi gher order elements in the
PAGE 100
87 problem area, and replacing them with linea r elements instead. No such averaging technique or element modi fication is employed here. There is a marked change in KI, and the other SIFs (KII and KIII) also had fairly substantial changes for this particular orienta tion and geometry. In this case KIII remains small across the crack front. But KII increases by quite a bit, just like KI does. KII at the midplane is roughly 0.04 psi in at room temperature. W ith the addition of temperature to this model, KII rises to 0.315 psi in. Granted for this geometry and set of loads, KII and KIII still remain small when compared to KI. Let us now consider a case where we have a larger gradient applied to the M(T) model in Figure 517. The next set of results for this material orientation are for T1 = 200oF and T2 = 100oF. In Figure 519 we see how the larger ther mal gradient impacts the mode I SIF. With these boundary conditions, KI, at the midplane, is 1.78 psi in. But this time, KII and KIII both show a fairly large increase after th is gradient is applied. Here KII at the midplane reaches 0.122 psi in and KIII = 0.355 psi in. Figure 517. Applicatio n of thermal loads
PAGE 101
88 KIKIII vs normalized crack front distance = 0 degrees0.5 0 0.5 1 1.5 2 2.5 00.20.40.60.81 normalized crack front distance psi in KI (DC) KI (F3D) KII (DC) KII (F3D) KIII (DC) KIII (F3D) Figure 518. Room temperature KIKIII (F3D) vs. KIKIII (DC) with T1 = 0oF and T2 = 100oF KIKIII vs normalized crack front distance = 0 degrees1.5 1 0.5 0 0.5 1 1.5 2 2.5 00.20.40.60.81 normalized crack front distancepsi in KI (DC) KI (F3D) KII (DC) KII (F3D) KIII (DC) KIII (F3D) Figure 519. Room temperature KI (F3D) vs. KI (DC) with T1 = 200oF and T2 = 100oF
PAGE 102
89 We see in the above simulations that a dding a thermal gradient can impact the mode IIII SIFs. In the first section of th is chapter, we looked at SIFs for various orientations for straight a nd diagonal cracks alike. When the crack is straight ( = 0), at room temperature, the KII and KIII are sma ll compared to KI. Again using the same geometry, when thermal loads are applied to the top and bottom faces of the M(T) model, KII and KIII are still small co mpared to KI, but they do exhibit a change from the room temperature reference state. As we have s een with the room temperature results, KII and KIII become more significant when the crack is inclined. This introduces a degree of load asymmetry, or mode mixity. To furthe r explore the influence of mode mixity along with the thermal gradient, let us consider the case two orientati on but with a 30 degree inclined crack. KIKIII vs normalized crack front distance = 30 degrees0.7 0.2 0.3 0.8 1.3 1.8 00.20.40.60.81 normalized crack front distancepsi in KI (DC) KI (F3D) KII (DC) KII (F3D) KIII (DC) KIII (F3D) Figure 520. = 30 degrees: room temperature KI (F3D) vs. KI (DC) with T1 = 0oF and T2 = 100oF
PAGE 103
90 Taking a look at KI at the midplane in Figure 520, KI goes from 1.37 psi in at room temperature to 0.978 psi in when the thermal gradient is present. KII and KIII are 0.545 and 0.188 psi in at the midplane, respectively. Our final set of plots is again for case two with a thirty degr ee inclined crack, but now T1 = 200oF and T2 = 100oF. KI at the midplane is approximately 1.74 psi in, and KII and KIII are 0.455 psi in 0.0394 psi in, respectively. KIKIII vs normalized crack front distance = 30 degrees1.75 1.25 0.75 0.25 0.25 0.75 1.25 1.75 2.25 00.20.40.60.81normalized crack front distance psi in KI (DC) KI (F3D) KII (DC) KII (F3D) KIII (DC) KIII (F3D) Figure 521. = 30 degrees: room temperature KI (F3D) vs. KI (DC) with T1 = 200oF and T2 = 100oF A final summary of the midplane values shown in the tables below. The abbreviations are as follows RT = room temperature, TG1 (T1 = 0oF and T2 = 100oF), TG2 (T1 = 200oF and T2 = 100oF). The values are taken at the midplane. The units for KIKIII are psi in.
PAGE 104
91 Table 55. Summary of midplane KIKII KI (RT) KI (TG1) KI (TG2) KII (RT) KII (TG1) KII (TG2) Case 2 ( = 0) 2.13 1.44 1.78 0.038 0.315 0.122 Case 2 ( = 30) 1.37 0.978 1.74 1.39 0.545 0.455 There are many more ways to apply boundary conditions to this problem. In the above analysis, since the specimen is 12 inches long, the temperature is different than room temperature, but still fairly uniform across the crack front. Table 56. Summary of midplane KIII KIII (RT) KIII (TG1) KIII (TG2) Case 2 ( = 0 ) 0.049 0.2580.289 Case 2 ( = 30 ) 0.242 0.2060.0394 The information presented in chapter tw o allows an estimation of the elastic constants over wide range of temperatures. Perhaps in some instances, the flaw can be six inches, or more, away from the surface of th e ET. In that case, given the foams very low thermal conductivity, it is po ssible to have a temperatur e over the crack front that does not vary substantially. But over t ypical acreage sprays the thickness of the insulation is roughly thre e inches. That is to say, ther e certainly are c onditions where the temperature across the crack front will not be uniform. To investigate this scenario let us reuse the same M(T) model as before but now the gradient will applied to th rough the thickness, or 1.5 inch dimension. Once again, the case two orientation along with a 30 degree inclined crack. The TG2 boundary conditions are applied and the resulting temp erature distribution al ong the crack front is plotted in Figure 522. Here the temperature varies appr eciably. To analyze this problem, the elastic constants are evaluated usi ng Figures 25 and 26. With those values the sixth order characteristic equation is solv ed at each point, and those roots are used in equation in (358) to obtain the corresponding Ks along the front.
PAGE 105
92 Crack front temperature (F) Case 2 orientation; = 30 degrees 250 200 150 100 50 0 50 100 150 00.20.40.60.81 normalized crack front distancetemperature (F) Crack front temperature (F) Figure 522. Temperature dist ribution along the crack front KIKIII vs normalized crack front distance = 30 degrees1.5 1 0.5 0 0.5 1 1.5 2 2.5 00.20.40.60.81 normalized crack front distancepsi in KI (DC) KII (DC) KIII (DC) Figure 523. KIKIII vs. normalized crack front distance; TG2 applied in thickness direction
PAGE 106
93 The results of this model are shown in Fi gure 523. Since the temperatures changes with respect to crack front distance, the elas tic constants also vary which impact the Ksolution. In some parts of the plot, KII is la rger than KI, but the mode I SIF dominates in the middle portion of the crack front. At the midplane, KI = 1.92 psi in. Comparing this to the midplane value in Figure 521, KI = 1.74 psi in; the value increased by approximately 10%. Summary The effect of material orientation on th e anisotropic Ksolution is examined. Seventeen material cases are examined and wh en the crack is straight we do not see a large difference in the mode I SIF even for or ientations that depart quite a bit from the norm. The SIFs are also computed for an incl ined crack. From the isotropic relations in (51) we see that KI should decrease while the remaining modes will increase as the gets larger. This trend is evid ent in Figures 59 through 511. Crack turning, at room temperature, was also estimated using both anisotropic and isotropic turning assumptions. In all cases the isotropic calculations predicted a higher turning angle than seen with in the three samples presente d here. Along with isotropic predictions, crack turning estimations that account for the direction dependence of the elastic constants and fracture toughness. Here the predictions are somewhat closer to the experimental results, but there is still a sizable difference. Perhaps one reason why this is so, is because when these specimens were fa bricated the material orientation was not carefully controlled. This makes the determ ination of the knit line orientation somewhat difficult and this can certainly influence the results because nearly all of the terms in Hoenigs equations depend on the elastic constants.
PAGE 107
94 Finally, the effect of temper ature is analyzed for the case two orientation. Several different sets of temperature boundary c onditions are applied along with the unit pressure. This analysis is conducted for both a st raight and inclined crack It is here that an important distinction can be made between the SIFs ev aluated at room temperature and those done so with the DC method with the applied te mperatures. This method, while simple to use, is not without its defici encies. Midplane results were compared for the case 1 orientation and a di fference of 10% between the mode I SIF is encountered. Thus, the numbers presented here should not be taken as absolutes and better agreement might be obtained by using more refined mode ls. Within the realm of the material orientations examined in this chapter, the m ode I SIF did not vary a great deal, but when there is a thermal gradient present, significan t changes in KI are encountered. The mode II and III SIFs also saw marked changes wh en the crack is no longer horizontal. In particular KIII, which is relatively small comp ared to KI and KII, for all of the room temperature tests showed larg e increases particularly when the second set of thermal boundary conditions are applied. In our last example, the thermal gradient is applied through the thickness where this dimension is much smaller than the height. Such a large gradient over a small thickness (1.5 inches) should translate in to a thermal gradient that va ries over the crack front. In this case the gradient does indeed change substantially and, unlike the prior thermal analysis, one can no longer assume a uniform temperature and prope rties along the crack tip. Thus, at each point the constants are re evaluated and the resulting Ks are plotted against the normalized crack front distance. In this scenario, KI is the largest at the midplane.
PAGE 108
95 CHAPTER 6 CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE WORK Conclusions Initial investigations of analytical solu tions for cracked anisotropic elastic solids have existed since the 1960s and the overwhelm ing majority of literature on this subject assumes the crack lies in a sy mmetry plane where antiplane decoupling of displacements takes place. Relatively few studies take on the issue of general anisotropy where no assumptions are made to simplify the mathema tical treatment of such a complex problem. In this study, two such methods are presente d to provide a solution to obtain neartip stresses for a cracked anisotropic elastic body. The first of which uses Hoenigs (1982) derivation of an energy release rate for genera l anisotropy. This is the preferred choice for the computation of SIFs because of the inherent accuracy in this method. As such, the Cornell Fracture Group developed a software package, FRANC3D, which is able to not only create crack meshes but also compute an anisotropic Ksolu tion for a variety of material orientations: isotropy, m onoclinic, and general anisotropy. Alternatively, using the DC method, another approach to determining the anisotropic Ksolution is deve loped to not only provide anot her general approach as far as material orientation goes but also one that can pr ovide a solution for problems containing both mechanical and thermal load s. FRANC3D does provi de a very versatile solution. But only structures subjected to isot hermal conditions can be analyzed with this software. However, the problem that NASA is encountering likely invo lves the effects of a thermal gradient within the foam insulation. It was intended that these methods would
PAGE 109
96 provide NASA with a more complete and funda mental understanding with regard to the fracture response of this material. The va rious conclusions are summarized below: 1. Numerical studies conducted at room temperature ove r a range of material orientations showed us that the mode I SIF did not depart substantially from the solution that arises from the orientation (case 0) which is present over the vast majority of the tank. The maximum diffe rence in the mode I SIFs between the prescribed orientations and the cas e 0 scenario was shown to be 14%. 2. The specialized divot tes ting conducted by NASA gave so me insight as to how subsurface crack, perhaps emanating from sharp void, could propagate toward the surface. It is here that crack turning an alysis could be useful. While the physics, loading, and geometry are quite complicated there are several th eories available in the literature that are app licable to this problem sin ce it is one involving brittle fracture. The fracture toughness and elastic co nstants are direction dependent in this type of foam material. As such, crack turning predictions that account for both of these characteristics were performed on th ree M(T) specimens provided by NASA. One of the three analytical predictions had reasonable agreement and while the remaining two had a large difference relativ e to the experimental value, they are closer than what isotropic assumptions would predict. 3. Temperature effects play an important role in this study. This is among the main reasons why NASA chose to initiate an exha ustive series of tests over range of temperatures. A large body of literature is available on foams, covering their applications and ways to estimate elastic a nd fracture response, but few studies deal with the harsh environment that this foam is subjected to. Since the range of temperatures involved in this problem is quite large, the effects of thermal gradients on the anisotropic Ksolution were studied. Several variables come into play with this type of analysis, namely the am ount of mode mixity, specimen geometry, material orientation, and boundary conditi ons (be it thermal, mechanical, or a combination of both). Many of the above variables were manipulated and the DC method developed was able to produce an anisotropic SIF soluti on that showed how the individual Ks varied ove r the crack front. It was shown that the SIFs could change dramatically depending on the level of mode mixity and the manner in which the thermal boundary conditions are applied. The results obtained in this study represent a quantitative basis for determini ng the response of this material for a given set of boundary conditions using the computational methods outlined in earlier chapters. Recommendations for Future Work One of the primary aims of this st udy was to help NASA acquire a more fundamental understanding of the fracture res ponse of the BX265 foam insulation. To some extent that was achieved as made evident by the conclusions made in the previous
PAGE 110
97 section. Nevertheless, much remains to further our understanding of such a complex material. What follows is a list of recommended work for future study: 1. Various models exist in the literature to evaluate the material properties and fracture response for foam materials. Such estima tions come from considering the geometry of the cells that comprise the cellular mate rial. While the density does vary over any given test specimen, or volume of foam, pe rhaps a wide range of measurements could be taken so an average can be made. This enables us to compute the relative density and this dimensionless parameter is quite us eful in estimating the elastic properties, even for anisotropic materials. 2. The crack turning investigation presented he re is far from comprehensive. A more complete testing program invol ving many different orientati ons and crack inclinations should be developed so better criteria can be developed for more accurate predictions. Within the three test samples studied, it a ppears that the crack is seeking a path normal to the maximum principal stress. Bu t with only three specimens on hand, it is difficult to generalize results. 3. While this study made use of simple and st andardized geometries, the problem NASA is faced with is one that involves a rather complicated geometry. Divot models can be built using FE software and crack turning analysis similar to what was done in this program can be performed. This should provide NASA with a very comprehensive representation of the actual problem since anisotropy and temperature effects can be accounted for. Many different geometries and loads could be applied and the resulting turning angles coul d be compared to experime nt. This could lead to additional refinements being made to the crack turning theories for anisotropic materials. 4. While the analysis in this study is confined to foam, there is no reason why it cannot be extended to other engineering materials th at exhibit direction de pendent properties. For example, perhaps crack turning calcula tions can be performed on various types of aluminum, or other metals where the crysta l orientation influences the behavior. Additional work can be done to modify the various theories, crack turning for example, to account for a plastic zone in front of the crack tip. Plasticity is neglected in this study because the foam is so brit tle, but many engineering materials are quite a bit more ductile which means that a process, or plastic zone, needs to be accounted for.
PAGE 111
98 APPENDIX GENERAL SOLUTION FOR THE PLANE PROBLEM OF A LINEAR ELASTIC ANISOTROPIC BODY Following Lehkhnitskii (1963), Sih et al (1965), Ting (1996) and Cherepanov (1979) the characteristic equation described in chapter three will be derived here. The ensuing discussion is for a plane probl em, meaning the stresses, strains, and displacements have no dependence on the z coor dinate. We start by listing the equations of equilibrium 0 0 0xyyxyyz xxzxyyxxy (A1) Since this a plane problem, the reduced compliance matrix becomes 33 33'ij ijijSS SS S The strains and displacements are related to each other 0 xyz xyxzyzuvw xyz uvww yxxy (A2) Finally we have the compatibility equations 22 2 22 yxyyz xxz y xxyyx (A3) Two arbitrary stress functions (U and F) are introduced as a means to solve the equilibrium equations
PAGE 112
99 22 22 2 xy xyxz yzUU y x UF x yy F x (A4) If we substitute (A4) into the compatibility relations we have a system of equations in terms of U and F 43 320 0LULF LULF (A5) Where 44444 4222612661611 432234 3333 3242546145615 3223 222 2554544 22'2'2''2'' '''''' '2''LSSSSSS x xyxyxyy LSSSSSS xxyxyy LSSS xxyy (A6) Through equation (A5) we can eliminate one of the stress functions, F for example, and it can be shown that 2 4230LLLU (A7) We note that the stress functions are in terms of the x and y coordinates and the 2 423LLLU operator is of order of order six. The solution to this problem will use a composite variable, z such that (,)()UxyFz (A8) zxy (A9)
PAGE 113
100 Before moving on with this discussion, let us define the term characteristic equation with a simple example. Consider, for a moment, a linear seco nd order homogeneous equation such as '''0yayby (A10) where a and b are constants, one method (Kre yszig, 1962) used to solve this type of differential equation is to assume the solution takes the form x y e (A11) with first and second derivatives 2' '' x x y eye (A12) If we substitute the derivates into equations (A12) we get 20ab (A13) This equation is called the characteristic, or auxiliary, equation w hose solution is readily determined by the quadratic formula. Here we ll will obtain two roots th at can be real and distinct, complex conjugates, or double real. The above pro cedure is not just valid for second order linear differential equations, but ca n be applied to another equation of order n. Continuing on, we now substitute (A 8) into (A7) which results in 2 243lll (A14) where 2 2554545 32 3151456254624 432 4111612662622'2'' '''''' '2'2''2'' lSSS lSSSSSS lSSSSSS (A15)
PAGE 114
101 The coefficients in (A14) are real a nd the corresponding roots are going to be complex and occur in conjugate pairs. The general case of the plane problem makes no assumptions that decouple the outofplane and inplane displacements. This decoupling can take place, depending on the type of materi al one is analyzing. Sih et al consider monoclinic materials. In this special cas e, Cherepanov derives an alternate solution where a fourth order character istic equation is developed 432 111312332322'2'2''2''0 SSSSSS (A16) 2 554544'2''0SSS (A17) The roots 1, 2 and their conjugates are determin ed from the first equation and 3 and its conjugate is obtained by the la tter. Finally, the neartip solution for such a material is 22 122121 12211221 12 12211221 11 ReRe 1111 ReRe 1 Re 2xx yy xy zx yzKIKII QQQQ KIKII QQQQ KI r 212 12121212 3 3 3111 Re Re 1 Re KII QQQQ KIII Q KIII Q (A18) 31323633zzxxyyxySSSA (A19) 2 111612'''jjj p SSS (A20) 122622'''jjjqSSS (A21) cossiniiQ (A22)
PAGE 115
102 1222112211 1212 1 21222112211 1212 3 3 4534411 ReRe 211 ReRe Re KIpQpQKIIpQpQ u r uKIqQqQKIIqQqQ u Q KIII CC (A23) where C45 and C44 are the elastic constants in stiffness form C (A24) where 1 ijijCS (A25)
PAGE 116
103 LIST OF REFERENCES Aliabadi and Rooke D.P. Numerical Fracture Mechanics, Boston: Computational Mechanics Publications, 1991, pp. 45. Anderson T.L. Fracture Mechanics, Boca Raton, FL: CRC Press, 1991, pp. 70, 120121, 430431. ANSYS Elements Reference, ANSYS Rel ease 5.6; ANSYS, Inc. November 1999. BanksSills L., Hershkovitz I., Wawrzynek P., Eliasi R., and Ingraffea T., Methods for Calculating Stress Intensity Factors in An isotropic Materials: Part 1Z = 0 is a Symmetric Plane, Engineering Fracture Mechanics, 2005, Vol. 72, pp. 23282358. Barnett D.M., and Asaro R.J., The Fracture Mechanics of SlitLike Cracks in Anisotropic Elastic Media, Journal of th e Mechanics and Physics of Solids, 1972, Vol. 20, pp. 353366. Barsom, and Rolfe S. Fracture and Fatigue Control in Structures: Applications of Fracture Mechanics, Woburn, MA: ButterworthHei nemann, 1999, pp. 55, 77, 8081. Barsoum R.S., On the Use of Isoparametric Finite Elements in Linear Fracture Mechanics, International Journal for Nu merical Methods in Engineering, 1976, Vol. 10, pp. 2537. Bowie O.L., and Freese C.E., Central Crack in Plane Orthotropic Rectangular Sheet, Engineering Fracture Mechanics, 1972, Vol. 8, pp. 4958. Bogy D.B., The Plane Solution for Anisot ropic Elastic Wedges Under Normal and Shear Traction, Journal of Applied Mechanics, 1972,Vol. 39, pp. 11031109. Boone T.J., Wawrzynek P.A., and Ingraffea A.R., Finite Element Modeling of Fracture Propagation in Orthotropic Ma terials, Engineering Fract ure Mechanics, 1987, Vol. 26, pp. 185201. Boresi A.P., Schmidt R.J., Sidebottom O.M. Advanced Mechanics of Materials, New York: John Wiley and Sons, Inc., 1993, pp. 36. Brezny R., and Green D.J., Factors Controlli ng Fracture Resistance of Brittle Cellular Materials, Journal of the American Ce ramics Society, 1991, Vol. 78, 10611065.
PAGE 117
104 Buczek M.B., and Herakovich C.T., A Normal Stress Criterion for Crack Extension Direction in Orthotropic Co mposite Materials, Journal of Composite Materials, 1985, Vol. 19, pp. 544553. Carloni C., Piva A., and Viola E., An Altern ative Complex Variable Formulation for an Inclined Crack in an Orthotropic Medium , Engineering Fracture Mechanics, 2003, Vol. 70, pp. 20332058. Chan S.K., Tuba I.S., and Wilson W.K., O n the Finite Element Method in Linear Fracture Mechanics, Engin eering Fracture Mechanics, 1970, Vol. 2, pp. 117. Cherepanov G.P. Mechanics of Brittle Fracture, New York: McGrawHill, 1979, pp. 7178. Chen C.S., Crack Growth Simulation and Re sidual Strength Predic tion in Thin Shell Structures, Ph. D Thesis, Cornell University, 1999. Dhondt G., MixedMode KCalculations fo r Anisotropic Materials, Engineering Fracture Mechanics, 2002, Vol. 69, pp. 909922. Dieter G.E. Mechanical Metallurgy, New York: McGrawHi ll Inc., 1976, pp. 5564. Erdogan F., and Sih G.C., On the Extensi on of Plates Under Plane Loading and Transverse shear, Journal of Basic E ngineering, 1963, Vol. 85, pp. 519527. Fowlkes C.W., Fracture Toughness Tests of a Rigid Polyurethane Foam, International Journal of Fracture, 1974, Vol. 10, pp. 99108. Freund L.B., Stress Intensity Factor Calcul ations Based on a Conservation Integral, International Journal of Solids and Structures 1978, Vol. 14, pp. 241250. Gent A.N., and Thomas A.G., The Deforma tion of Foamed Elastic Materials, Journal of Applied Polymer Science, 1959, Vol. 1, pp. 107113. Gent A.N., and Thomas A.G., Mechanics of Foamed Elastic Materials, Rubber and Chemistry, 1963, Vol. 36, pp. 597610. Gibson L.J., and Ashby M., The Mechanics of ThreeDimensional Cellular Solids, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 1982, Vol. A 382, pp. 4359. Gibson L.J., and Ashby, M. Cellular Solids, Elmsford, NY: Pergamon Press Inc., 1988, pp. 7681, 241256. Gdoutos E.E. Problems of Mixed Mode Crack Propagation, The Hague, The Netherlands: Martinus Nijho ff Publishers, 1984, pp. 2.
PAGE 118
105 Hoenig A., Near Tip Behaviour of a Crack in a Plane Anisotropic Elastic Body, Engineering Fracture Mechanics, 1982, Vol. 16, pp. 393403. Huang J.S., and Gibson L.J., Fracture Toughness of Brittle Honeycombs, Acta Metallurgica Materialia, 1991, Vol. 39, pp. 16171626. Hussain M.A., Pu S.L., and Underwood J.H., Strain Energy Release Rate for a Crack Under Combined Mode I and Mode II, Frac ture Analysis ASTM STP, 1974, Vol. 560, pp. 228. Hyer M.W. Stress Analysis of FiberRein forced Composite Materials, New York: McGrawHill Inc., 1998, pp. 6975. Ingraffea A.R., and Manu C., Stress Intensity Factor Computation in Three Dimensions with QuarterPoint Elements, International Journal for Numerical Methods in Engineering, 1980, Vol. 15, 14271445. Irwin G.R., Analysis of Stresses and Strains near the End of a Crack Traversing a Plate, Journal of Applied Mechanic s, 1957, Vol. 24, pp. 361364. Irwin G.R., Fracture Mechanics, In: G oodier J.N., and Hoff N.J., (editors), Structural Mechanics: Proceedings of First Naval Symposium, Pergamon Press, 1960, pp.557. Ishikawa H., A Finite Element Analysis of Stress Intensity Factors for Combined Tensile and Shear Loading by Only a Vi rtual Crack Extensi on, International Journal of Fracture, 1980, Vol. 16, pp. 243246. Ko W.L., Deformations of Foamed Elasto mers, Journal of Cellular Plastics, 1965, Vol.1, pp. 4550. Kuo M.G., and Bogy D.B., Plane Solutions for the Displacement and TractionDisplacement Problem for Anisotropic Elastic Wedges, Journal of Applied Mechanics, 1974, Vol. 41, pp. 197203. Kreyszig E. Advanced Engineering Mathematics, New York: John Wiley and Sons Inc., 1962, pp. 105105, 136. Lekhnitskii S.G. Theory of Elasticity of an Anisotropic Body, San Francisco: HoldenDay, Inc., 1963, pp. 3, 104110, 117119, 122123. Li F.Z., Shih C.F., and Needleman A., A Comparison of Met hods for Calculating Energy Release Rates, Engineering Fract ure Mechanics, 1985, Vol. 21, pp. 405421. Logan D.L. A First Course in the Finite Element Method, Pacific Grove, CA: Wadsworth Group, 2002, pp. 413414, 433434.
PAGE 119
106 Maccagno T.M., and Knott J.F., The Fracture Behaviour of PMMA in Mixed Modes I and II, Engineering Fracture Mechanics, 1989, Vol. 34, pp. 6586. Maccagno, T.M., and Knott J.F, The Low Temp erature Brittle Fracture Behaviour of Steel in Mixed Modes I and II, Engineer ing Fracture Mechanics, 1992, Vol. 38, pp. 111128. Maiti, S.K., Ashby M.F., and Gibson L.J., Fracture Toughness of Brittle Cellular Solids, Scripta Mettalurgica, 1984, Vol. 18, pp. 213217. Man K.W., Aliabadi M.H., and Rooke D.P ., BEM Frictional Contact Analysis: Modeling Considerations, Engineering Analysis with Boundary Elements, 1993, Vol. 11, pp. 7785. McIntyre A., and Anderton G.E., Fracture Properties of a Rigid Polyurethane Foam Over a Range of Densities, Polymer, 1979, Vol. 20, pp. 247253. Moran B., and Shih C.F., A General Treat ment of CrackTip Contour Integrals, International Journal of Fracture, 1987, Vol. 35, pp. 295310. Morgan J.S., Wood J.L., and Bradt R.C., Ce ll Size Effects on the Strength of Foamed Glass, Materials Science and Engi neering, 1981, Vol. 47, pp. 3742. Nobile L., and Carloni C., Fr acture Analysis for Orthotropic Cracked Plates, Journal of Composite Structures, 2005, Vol. 68, pp. 285293. Patel M.R., and Finnie I., Str uctural Features and Mechan ical Properties of Rigid Cellular Plastics, Journal of Materials, 1970, Vol.5, pp. 909932. Pettit R., Crack Turning and Arrest in Integr ally Stiffened Aircraft Structures, Ph. D Thesis, Cornell University, 2000. Rice J.R., A Path Independent Integral and the Approximate Analysis of Strain Concentration by Notches and Cracks, J ournal of Applied Mechanics, 1968, Vol. 35, pp. 379386. Rathod H.T., A Study of Asymmetrically Load ed Griffith Crack in a Plane Anisotropic Elastic Body, Engineering Fracture M echanics, 1979, Vol. 11, pp. 8793. Saouma V.E., and Sikiotis E.S., Stress Inte nsity Factors in Anis otropic Bodies Using Singular Isoparametric Elements, Engineer ing Fracture Mechanics, 1986, Vol. 25, 00. 115121. Sih G.C., Paris P.C., and Irwin G.R., On Cr acks in Rectilinearly Anisotropic Bodies, Engineering Fracture Mechanics, 1965, Vo1. 1, pp. 189203. Sih G.C., Strain Energy Density Factor Applied to Mixed Mode Crack Problems, International Journal of Fracture, 1974, Vol. 10, pp. 305321.
PAGE 120
107 Sneddon IN., The Distribution of Stress in the Neighbourhood of a Cr ack in an Elastic Solid, Proceedings, Royal Society of London, 1946, Vol. A187, pp. 229260. Tada, H., Paris, P., and Irwin, G. The Stress Analysis of Cracks Handbook, Hellerton, PA: Del Research Corp., 1973. Ting T.C. Anisotropic Elasticity: Theory and Applications, Oxford: Oxford University Press, 1996, pp. 5455, 118133. Timoshenko S.P., and Goodier J.N. Theory of Elasticity, Tokyo: McGrawHill, 1982, pp. 354380. Tupholme G.E., A Study of Cracks in Orthotro pic Crystals Using Di slocation Layers, Journal of Engineering Mathematics, 1974, Vol. 8, pp. 5769. Wells A., Unstable Crack Propagation in MetalsCleavage and Fracture, Cranfield Crack Propagation Symposium, 1961, Vol. 1, pp. 210. Westergaard H.M., Bearing Pressures on Cracks, Journal of Applied Mechancis, 1939, Vol. 6, pp. 4953. Williams M.L., On the Stress Distribution at th e Base of a Stationary Crack, Journal of Applied Mechanics, 1957, Vol. 24, pp. 109114. Yau J.F., Wang S.S., and Corten H.T., A MixedMode Crack Anal ysis of Isotropic Solids Using Conservation Laws of Elasti city, Transactions of the ASME, 1980, Vol. 47, pp. 335341.
PAGE 121
108 BIOGRAPHICAL SKETCH Erik Knudsen was born in West Chester, PA, on February 18th, 1979. At the age of four his family moved to Naples, FL, wher e he spent the remai nder of his childhood and graduated from Lely High School in 1997. He earned a B.S. in mechanical engineering from Villanova in the Spring of 2001 and returned to Florida for his graduate work. Under the tutelage of Dr. Nagaraj Arakere he earned a M.S. in mechanical engineering in the Fall of 2003 working on a NASA related pr oject. Upon completion of the doctorate from the University of Florida, he would lik e to become a fulltime teacher at the college level or perhaps retain a positi on in the aerospace industry.

