Citation
Anisotropic Fracture Analysis of the BX-265 Foam Insulation Material Under Mixed-Mode Loading

Material Information

Title:
Anisotropic Fracture Analysis of the BX-265 Foam Insulation Material Under Mixed-Mode 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 non-profit 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:

knudsen_e ( .pdf )

knudsen_e_Page_002.txt

knudsen_e_Page_113.txt

knudsen_e_Page_020.txt

knudsen_e_Page_053.txt

knudsen_e_Page_006.txt

knudsen_e_Page_073.txt

knudsen_e_Page_052.txt

knudsen_e_Page_057.txt

knudsen_e_Page_007.txt

knudsen_e_Page_049.txt

knudsen_e_Page_094.txt

knudsen_e_Page_011.txt

knudsen_e_Page_116.txt

knudsen_e_Page_024.txt

knudsen_e_Page_085.txt

knudsen_e_Page_025.txt

knudsen_e_Page_030.txt

knudsen_e_Page_107.txt

knudsen_e_Page_056.txt

knudsen_e_Page_096.txt

knudsen_e_Page_037.txt

knudsen_e_Page_069.txt

knudsen_e_Page_086.txt

knudsen_e_Page_058.txt

knudsen_e_Page_108.txt

knudsen_e_Page_012.txt

knudsen_e_Page_106.txt

knudsen_e_Page_063.txt

knudsen_e_Page_031.txt

knudsen_e_Page_046.txt

knudsen_e_Page_054.txt

knudsen_e_Page_064.txt

knudsen_e_Page_028.txt

knudsen_e_Page_060.txt

knudsen_e_Page_016.txt

knudsen_e_Page_081.txt

knudsen_e_Page_103.txt

knudsen_e_Page_040.txt

knudsen_e_Page_117.txt

knudsen_e_Page_042.txt

knudsen_e_Page_022.txt

knudsen_e_Page_051.txt

knudsen_e_Page_070.txt

knudsen_e_Page_023.txt

knudsen_e_Page_092.txt

knudsen_e_Page_039.txt

knudsen_e_Page_013.txt

knudsen_e_Page_048.txt

knudsen_e_Page_034.txt

knudsen_e_Page_038.txt

knudsen_e_Page_072.txt

knudsen_e_Page_047.txt

knudsen_e_Page_102.txt

knudsen_e_Page_014.txt

knudsen_e_Page_111.txt

knudsen_e_Page_077.txt

knudsen_e_Page_065.txt

knudsen_e_Page_114.txt

knudsen_e_Page_101.txt

knudsen_e_Page_066.txt

knudsen_e_Page_008.txt

knudsen_e_Page_087.txt

knudsen_e_Page_119.txt

knudsen_e_Page_043.txt

knudsen_e_Page_115.txt

knudsen_e_Page_062.txt

knudsen_e_Page_041.txt

knudsen_e_Page_120.txt

knudsen_e_Page_105.txt

knudsen_e_Page_010.txt

knudsen_e_Page_074.txt

knudsen_e_Page_004.txt

knudsen_e_Page_050.txt

knudsen_e_Page_071.txt

knudsen_e_Page_015.txt

knudsen_e_Page_082.txt

knudsen_e_Page_005.txt

knudsen_e_Page_095.txt

knudsen_e_Page_059.txt

knudsen_e_Page_001.txt

knudsen_e_Page_075.txt

knudsen_e_Page_003.txt

knudsen_e_Page_110.txt

knudsen_e_Page_079.txt

knudsen_e_Page_027.txt

knudsen_e_Page_100.txt

knudsen_e_Page_088.txt

knudsen_e_Page_035.txt

knudsen_e_Page_033.txt

knudsen_e_pdf.txt

knudsen_e_Page_026.txt

knudsen_e_Page_029.txt

knudsen_e_Page_098.txt

knudsen_e_Page_061.txt

knudsen_e_Page_044.txt

knudsen_e_Page_018.txt

knudsen_e_Page_099.txt

knudsen_e_Page_078.txt

knudsen_e_Page_093.txt

knudsen_e_Page_084.txt

knudsen_e_Page_045.txt

knudsen_e_Page_036.txt

knudsen_e_Page_121.txt

knudsen_e_Page_080.txt

knudsen_e_Page_032.txt

knudsen_e_Page_112.txt

knudsen_e_Page_089.txt

knudsen_e_Page_097.txt

knudsen_e_Page_091.txt

knudsen_e_Page_009.txt

knudsen_e_Page_090.txt

knudsen_e_Page_076.txt

knudsen_e_Page_019.txt

knudsen_e_Page_055.txt

knudsen_e_Page_083.txt

knudsen_e_Page_118.txt

knudsen_e_Page_067.txt

knudsen_e_Page_068.txt

knudsen_e_Page_021.txt

knudsen_e_Page_104.txt

knudsen_e_Page_109.txt

knudsen_e_Page_017.txt


Full Text












ANISOTROPIC FRACTURE ANALYSIS OF THE BX-265 FOAM INSULATION
MATERIAL UNDER MIXED-MODE 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 BX-265 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 Factor-K............... ..... ............... 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

2-1 Measured material property values from experiment ........................ ...........18

3-1 Various m ethods for determ ining the SIF ..................................... .................38

4-1 Direction cosines .......... ....... ........ .................... ................... 66

5-1 Direction cosines for the case one material orientation .......................................73

5-2 KI for cases 0-8, 0 degree inclination ........................................... ............... 76

5-3 KI for cases 9-16, 0 degree crack inclination ................ .................................. 76

5-4 Measured and predicted turning angles (in degrees)................... ... .............83

5-5 Sum m ary of m id-plane K I-K II......... ............... ........................... .. ...................91

5-6 Sum m ary of m id-plane K ill......... ............... .................................. ............... 91
















LIST OF FIGURES

Figure page

1-1 Foam loss during Shuttle ascension ........................................ ........ ............... 2

1-2 Spray-on foam insulation (SO FI)................................................................... ......3

2-1 Idealized foam microstructure...................................................... ........................ 8

2-2 Square specimen tensile test used to determine the Poisson ratio..........................15

2-3 F oam test specim ens........... ............................................................ ...... .... ... ... 15

2-4 Summary of stress-strain data .............................................. 16

2-5 Coordinate system for the transversely isotropic material............................... 17

2-6 Young's moduli vs. temperature .................... ........................... 18

2-7 Shear m oduli vs. tem perature...................... ................................. ............... 19

2-8 Coefficient of thermal expansion vs. temperature ...........................................19

2-9 Thermal conductivity vs. temperature....... .......... ................... ..............20

2-10 From left to right: M(T), C(T), and SNE(B) fracture test specimens ....................22

2-11 Summary K-values for various fracture specimens ...........................................23

2-12 Clip gauge on C(T) specimen used to measure CMOD .....................................24

2-13 Typical load vs CMOD curve for a brittle material ...........................................25

2 -14 D iv ot test p an el............ .... ........................................................................... .. ...... .. 2 6

2-15 2-D view of divot test................................................................. ............... 27

2-16 2-D view of test panel after heat flux has been applied .......................................27

2-17 Foam blow out from divot test ............................... ............... 28

3-1 Coordinate system used in equation (3-1)................................... ............... 31









3-2 From left to right: mode I (opening), mode II (shearing), mode III (tearing) .........32

3-3 Path of the J-integral .............................................. .... .... .. ........ .... 39

3-5 D definition of Irwin's crack closure integral ............................ .............................43

3-6 Correlation points typically used for displacement correlation ............................46

3-7 Triangular quarter point elem ents ........................................ ........................ 48

3-8 Isoparametric element and degenerated element with mid-side nodes moved to
quarter point locations .............................................................................48

4-1 D definition of crack turning angle ........................................ ......................... 57

4-2 Elliptical relationship for T ............................................................................58

4-3 D definition of hoop stress ................................... .... ............................ ...............59

4-4 Orthotropic toughness values ........................................................ ............. 62

4-5 Definition of the a and n vectors ................................................... ............... 62

4-7 M (T) m odel created in AN SY S ........................................ .......................... 65

4-8 Sam ple transform action ................................................. ................................ 66

4-9 Fractured M(T) specimen: the dotted line is added to show the location of the
k n it lin e ............................................................................ 6 7

5-1 Definition of crack inclination angle.................................................................. 69

5-2 D definition of crack front distance....................................... .......................... 69

5-3 M (T) m odel built in A N SY S ........................................................................ .. .... 70

5-4 D definition of m material orientation..........................................................................70

5-5 Top view of the cones shown in Figure 5-4......................................................71

5-6 Hypothetical material orientation relative to the ET..............................................71

5-7 Definition of case one material orientation for the M(T) FE model ......................72

5-8 Mode I SIF for cases 1-4; 0 degree crack inclination ...........................................74

5-9 Mode I SIF for cases 1-4; 30 degree crack inclination..................... ...............74

5-10 Mode II SIF for cases 1-4; 30 degree crack inclination .......................................75









5-11 Mode III SIF for cases 1-4; 30 degree crack inclination.............................75

5-12 From left to right: front, left, right, and rear sides of the M(T) specimen..............77

5-13 D eterm nation of knit line plane..............................................................................78

5-14 M (T) m odel with inclined crack ..................................................... ............. 79

5-15 Cartesian and cylindrical stresses................................ ................... ...... ........ 81

5-16 Crack turning angle .................. ............................... ..... .. ............ 82

5-17 Application of therm al loads .............................................................................87

5-18 Room temperature KI-KIII (F3D) vs. KI-KIII (DC) with T1 = 0F and T2 =
10 0 F ............................................................................... 8 8

5-19 Room temperature KI (F3D) vs. KI (DC) with T1 = -2000F and T2 = 100F..........88

5-20 (p = 30 degrees: room temperature KI (F3D) vs. KI (DC) with T1 = 0F and T2 =
10 0 F ............................................................................... 8 9

5-21 (p = 30 degrees: room temperature KI (F3D) vs. KI (DC) with T1 = -2000F and
T 2 = 10 0 F ........................................................................ 9 0

5-22 Temperature distribution along the crack front..................................................... 92

5-23 KI-KIII 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 BX-265 FOAM INSULATION
MATERIAL UNDER MIXED-MODE 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 closed-cell

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 spray-on 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 re-entry. 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 non-uniform 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 mixed-mode 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 BX-265 foam insulation material.














CHAPTER 1
BACKGROUND

Problem Motivation

On February 1st, 2003, the Space Shuttle Columbia suffered a catastrophic failure

during re-entry. NASA has conducted an exhaustive investigation of the failure and the

consensus now is the breakup was caused by a segment of BX-265 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 re-entry, 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 1-1).


























Figurel-1. 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

1-1 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 1-2).

Pictures of the Shuttle during liftoff, such as those shown in Figure 1-1, 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 1-2. Spray-on 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, crack-like, defects for this sort of material is

through linear elastic fracture mechanics (LEFM). For isotropic materials, the near-tip

stresses are completely characterized by K-the stress intensity factor (SIF). The

analytical K-solutions for anisotropic materials take a similar form, but have additional

terms that come from the various elastic moduli. The formulations of near-tip stresses for

isotropic materials are well established and are now incorporated in many finite element

(FE) software packages. Some K-solutions for certain kinds of direction-dependent

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 single-crystal 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 K-solutions are not readily

available, the Cornell Fracture Group has developed a software package, FRANC3D, that

simplifies the meshing process necessary for FE modeling of crack-like defects in elastic

solids and also computes a K-solution 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 K-solutions 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 K-solution is
examined

2. The impact of mode mixity on the K-solution 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
K-solution 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 two-dimensional arrays of

polygons that fill a planar area and 2) foams, the three-dimensional 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 2-1 (Gibson and Ashby, 1988) lists

several properties of polymer foams and true solids (solid metals and ceramics).

Table 2-1. Properties of solids and polymer foams

Density (lb/in3) Conductivity (BTU/hr in F) Young's modulus (psi)
Solids 3.6 x 10-2- 3.6 x 10- 4.8 x 10 48 3 30 x 106
Polymer 3.6 x 10-4 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 2-1) The equations that

arise from this mathematical treatment are in terms of the strut dimensions, which in turn




























Figure 2-1: 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 membrane-like faces. Within open-celled 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 (2-1)
P, I

where t is the cell wall thickness and the 1 is the edge-length. For open-cell foams

-= C2 (2-2)
P, 1

Finally, for the foams of the closed-cell variety we have

p* t
=C3 (2-3)
P, I

The C terms are numerical constants that depend on the cell shape, normally taking a

value of unity.

As mentioned beforehand, certain open-celled 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 (2-4)
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 (2-2). The second moment of the area, I, can also be

related to the cell wall dimensions









I oc t4 (2-5)

Now (2-4) can be re-written as

E C,P (2-6)


Using the same reasoning, we can write down the relationship for the shear moduli

P =C2 (2-7)


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- ) (2-8)
2(1+ v)

Using (2-8) along with (2-6) and (2-7) we have

v* C= 1 = C (2-9)
2C,

The above analysis is for open-celled foams. When examining closed-cell 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 open-celled 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 closed-cell 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 (2-10)
2 /f

Equation (2-10) can be re-written when we note that I t4e and E* ocx (F/12)/(6/1) which

yields


=a' +/f' (2-11)
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,
(2-12)

t 0.9301/2 0P*
I P

If we substitute (2-12) into (2-11) 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-: ), (2-13)
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*) (2-14)



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 closed-cell foam is shown below

G=C22 p +2(1-) (2-15)
E P, P,

Finally the last contribution to consider comes from the bending on the cell struts, and

that analysis is identical to the open-celled formulation presented earlier.

We note that the above relations are for foams where the cell walls are equiaxed.

Most man-made 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
(2-16)

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 open-celled foams as

E *3 2R2
3- 2(2-17)
E*1 (1+(1/R)3)

2R
For closed cells an additional term (1- +) appears in (2-17). For the
1+ (1 / R)

anisotropic shear modulus, we have

3 (2-18)
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 BX-265 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 all-important 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

un-foamed 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 2-2 and 2-3 display some of the various test specimens used in various

tension and torsion tests. A summary of typical stress-strain data is shown in Figure 2-4.

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 2-2. 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 2-3. Foam test specimens







16



Summary of Typical Tensile Stress Strain Curves
BX-265 / Category 4


002 004 006 008 01 012
Strain (inlin)


Figure 2-4. Summary of stress-strain 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 non-linear

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 11-22 plane in Figure 2-5 is a plane of isotropy.


Source: personal communication, Doug Wells (MSFC)






















Figure 2-5. Coordinate system for the transversely isotropic material

Using the coordinate axes from Figure 2-4 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 2-1 lists the room temperature values of the elastic constants for the orientation

shown in Figure 2-5.

Since one of the objectives of this study is to analyze the effect of temperature, or

thermal loads, on the K-solution 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 (2-19)

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 (2-20)
V13 V31
,3 12 21
E3 El

Table 2-1. 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 2-6 through 2-9. 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.

-0-Ell
Young's moduli --- E22
-- E33
6000

5000

4000

'. 3000

2000

1000

0
-550 -450 -350 -250 -150 -50 50 150 250
temperature (F)



Figure 2-6. 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 2-7. 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 2-8. 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 2-9. Thermal conductivity vs. temperaturet

Fracture Testing for the BX-265 Foam Material

While many foam materials do not exhibit linear-elastic response before yielding,

brittle foams normally exhibit linear-elastic behavior (while in tension) up until fracture

takes place. Thus, we would like to use LEFM concepts to compute the near-tip 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 (2-21)









where Cs is a constant of proportionality and cfs is the failure strength of the struts.

Equation (2-21) can also be re-written as

KI, =C, KIJs (p/p* )3/2 (2-22)

where KIcs is the toughness of the solid material. Equation (2-22) 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 (2-21). 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 well-known 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 single-edge 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, continuum-based relations to determine fracture properties,

such as the plane strain toughness. Pictures of these test specimens are shown in Figure

2-10. 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 2-11

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 2-10. 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 toughness-KIc.

If the material is brittle and the subsequent load vs crack mouth opening

displacement (CMOD) curve (Figure 2-12) is linear-elastic 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 2-13) 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)(A-S) M(T)(A-) M(T)(33-S) C(T)(A.S) C(T)(33-S)


Figure 2-11. Summary K-values for various fracture specimens









P
K B f (a/W) (2-5)




Clip




A CMOD







Figure 2-12. 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 2-13. 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 defect-free, 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

2-14 are used for the 'divot testing.'


















Figure 2-14. 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. Two-dimensional schematics of this process are shown in Figures

2-15 and 2-16.

The applied flux warms the liquid nitrogen and eventually a phase change takes

place. The experiment takes place in a thermo-vacuum 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 2-15. 2-D view of divot test

In Figure 2-17, 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 2-16. 2-D 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 2-16, by calculating near-tip 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 2-17. 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 well-known 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 cylindrical-shaped 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 near-tip 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 near-tip

stresses for materials that are anisotropic, or have direction-dependent 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 (3-1).

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- (3-1)
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

near-tip stress fields for anisotropic materials tried to take advantage of material

symmetry to decouple 'in-plane,' and 'out-of-plane' displacements, which, in turn, makes

the mathematical formulation a little less complicated.

Isotropic Fracture Mechanics

For isotropic materials, the near-tip 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 far-field 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 (3-1)

y



P(x,y)

Crack faces





Figure 3-1. Coordinate system used in equation (3-1)

The coordinate system for equation (3-1) is shown in Figure 3-1, 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 (3-1) are that the near-tip 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 3-2 displays the three modes of deformation.













Figure 3-2. 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
(3-2).

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+l-2cos2 -cos-2 K-l-2sin2j 0 KII (3-2)
u 2, 2 2s 2J 2l 2( 0
Zu KIII
0 0 sin-
2









where K = 3-4v for plane strain and K = (3-v)/(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 direction-dependent. The analytical work

geared at developing the near-tip 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 near-tip stresses for isotropic

materials, Sih et al. use a similar approach for materials with direction-dependent

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 (3-3).

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 in-plane and out-of-plane

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 (3-7)
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 ', (3-8)

q, = S',12 / S'26 + S'22/ (3-9)

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 (3-4)
KIJJ


(3-5)


U
/dzI


(3-6)










S' = S, (3-10)
S Y S33


Q, = Jcos ) +/ sin (3-11)

The C45 and C44 terms in (3-7) are the elastic constants which are in stiffness form

cr = Cec (3-12)

where

C, = S, (3-13)

The aforementioned decouplingg' of the anti-plane displacement is now evident in

equation (3-7). The Uz displacement does not depend on KI and KII. This is similar to

the isotropic case in equation (3-2); 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 anti-plane

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 (3-14)
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 (3-15)


u = Re m NJlK, N (cos0 + j sin)) (3-16)


Here an important distinction is made between the general case and Sih et al's

formulation. Notice in (3-16) how the mode I-III SIFs are to be summed over all three

modes of displacement. Thus, all three components of displacement are dependent on

KI-KIII.

n 1, = S1 ', 2 S'16+S' +A, (S'15 S'14)
=' -S\ 24q (3-17)
m,, = S', -S'26- S2 +A, S S,24(3-17)


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 (3-18). The coordinate system used in equations (3-14) and (3-16)

is the same as the one used in (3-1).

14 (/) ()32 (/)=0 (3-18)

Where

1 ) = SLp-2 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 (3-14) is defined as


N, = -p -2 -J3 (3-20)


Where j is defined as 13 (02)/l(u,).

There are two situations, also called degenerate cases, where this solution is

invalid, however. If anti-plane shear and plane strain displacements becomes decoupled-

for monoclinic materials, say-Hoenig'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 N-matrix, equation (3-20), becomes singular

because l1 and 02 = i. However, we have near-tip solutions for the isotropic cases.

These are shown in equations (3-1) and (3-2).

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









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 3-1.

Table 3-1. 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 K-solutions 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. 2-D problems involve discretizing the line-boundary of the domain and

for 3-D 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 near-tip stress and

displacement fields, Rice's J-integral (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 K-solution. This particular method of solution

would fall under Stage 3 in Table 3-1, and in this case Rice's J-contour 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; K-solutions can be generated for isotropic and fully anisotropic materials

alike.

Since many FE models were built and the subsequent K-solutions computed with

this software, it would be prudent to discuss the pertinent theory (Banks-Sills 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) (3-21)

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 (3-22)


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 3-3

y
n


~~=== ---V-x


Figure 3-3. Path of the J-integral









The J-integral was initially formulated for 2-D problems, however since this

material has properties that are direction dependent, our analysis requires the use of 3-D

FE models. The J-integral, then, needs to be modified so it can be used in integration

that take place over volumes. Following Li et al (1985), the J-integral can be rewritten as


iJ= cO, W8 q ds (3-23)


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 J-integral is now evaluated in a cylindrical domain centered over the crack

front.

Crack front
X2

Virtual crack advance










X3

Figure 3-4. Cylindrical domain and virtual crack advance for a 3-D crack front

Following Figure 3-4 we can define incremental area of virtual crack advance as

A = Aama q (s) ds (3-24)
AL

The definition of J in equation (3-22) 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 2-D problems) because they are

easier to implement numerically. So, the q term in (3-23) is simply a mathematical tool

that enables us to recast the J-integral into a slightly more usable form.

This virtual extension method is the modem-day 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 3-D 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 = (3-25)
Sq,(s)ds Aq

where qt(s) is the crack tip perturbation.

Rice's J-integral 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 J-integral to compute the Ks individually.

Ishikawa's method, though, has not been extended for anisotropic materials. Yau et al.

(1980) uses the M-integral 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) (3-26)
S= (1) u (2) Kill =KIII) + KIII2









The J-integral can now take the form

J= J1 + (2 +M1,2) (3-27)

The M(1,2) term is known as the M-integral 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 (3-28)
F a 9x, a9 ] }9x

This integral describes the interaction between the two equilibrium states, and sometimes

the M-integral 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 3-5.

Now a compressive, or closure stress, is applied in such a fashion to clamp the crack

faces down along the length 6-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 the crack. Substitutions are made for near-tip

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 3-30) 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 lim-a (3 r, 0)-.u (r, r)dr (3-29)
0


















IX





Compressive stress








Figure 3-5. Definition of Irwin's crack closure integral

Now, if the stress and displacement fields are substituted in (3-29), we have

G= = [KI m(m,,N, 1K,)+KIIIm(ml,N, 1K,)+KIIIm(m,,N, 1K)] (3-30)

When we use the relations for SIFs in (3-26), we get

1 (KI( + KI(2) 1M(K2 N ) + (KII() + KII(2))m(m2 N, -1K,)
G 2 +(KIII(1) + KIII2))m(mN,1K) (3-31)


The terms in (3-27) have bars over them because a certain definition of J is used in (3-

23). Let us define an alternate form of (3-27) where J is not normalized with respect to

the extended area as

J = (1) + (2) + M(1,2) (3-32)








Now we substitute (3-31) into (3-32) and for clarity the individual terms of that equation

are listed below


2 (3-33)
(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)) (3-34)
+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 (3-28) and (3-34). 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 (3-14). Now the two

relations for the M-integral are equated

KI(1) Im(m2,Nj-1K (2))+
"(2) d(1)(2) zN1K))+

x 01 0 1 KII() m(mlN 'K (2))
(3-35)
A, 2 KII(2) Im(1,N- 1K (l))+
KIII() im(m3, N 1Kj (2)+
KIIIl(2) m(m3, N-1K (2)









From (3-35) 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 (3-36)
Im m2^N 1 +ml^N1) =Imhn^N^ + N )
Im m2A31 +M3I1) Im MilA31 +M3 21 2Im (M3 331) K-I(1) (1,2) A


Finally, we note that (3-36) 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 (3-36) 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 3-2. 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 (3-2) with co = t7 and Figure 3-6 for a plane strain problem dealing

with an isotropic material, we have

K b- (uyb ,) 2u (3-37)
K+1 V r

(u(uxb-)2l/ 2z (3-38)
K+1 V r









where [s, again, is the shear modulus determined by E/2(1+v) and K = 3-4v for plane

strain We note that Uy and Ux are the opening and sliding modes of displacement. The

above method is sometimes called a 'one-point' matching technique since it involves

extracting displacements from one point near the crack tip. The downside of such a






a
y








Figure 3-6. 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 (3-1), 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 mid-side 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 3-7, 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 mid-side nodes moved to the quarter point locations. Consider the

bottom edge of the element on the right in Figure 3-8. 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) (3-39)

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 mid-side nodes

N, = 1 ( 2)(1+ ) t=-1, (i= 5,7)
2 (3-41)

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+) (3-42)
=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 3-7. Triangular quarter point elements


7 3 4



L 6


1 x 5 2


7 3


8


-/-- 3L/4
L/4


Figure 3-8. Isoparametric element and degenerated element with mid-side nodes moved
to quarter point locations

The global x coordinate is related to the isoparametric coordinate system via


(3-43)


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,


(3-44)









We can now solve for


= 2J -1 (3-45)

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, (3-46)

Again determining the displacements along the edge containing nodes one, two, and five

we write


2 2
u= 1-( l + (l++)u,+(1_- (3-47)


Now we can substitute (3-45) into (3-47) to obtain


u-- 1+2 2-2 u+ -1+2 2 u2+4 5 (3-48)


We note that the displacement shown in (3-48) has a /lx term present. This is consistent

with equation (3-2) where the displacements vary with /r. We need to differentiate

equation (3-48) to obtain the strain in the x direction

u 1 3 4 1 4 U 2 4
SCJ7 q. 4 I+- /_+- u2+ U5 (3-49)
a x 2[xL L 2 f Zj 2L fxL- L

Finally, it is shown in (3-49) how the strain displays the desired singularity (x-1/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 (3-37) and (3-38) are

modified because r = L/4


KI -(ub )4p 2 (3-50)
K+1 L


KI= (U -Ux,)4p 27+l (3-51)
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 3-D elements. This formulation is

particularly useful to the present study since 3-D simulations are used to model the

fracture response of the foam material. Ingraffea and Manu use a special labeling

convention shown in Figure 3-10 for a collapsed 20 node brick element (Figure 3-9), 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 3-9. 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) (3-52)
8

For the nodes at positions i = 17, 18, 19, 20 and ji = 0, ri = 1, (i = 1


N, = (3-53)
4

At i = 10, 12, 14, 16 and i = +1, ri = 0, (i = +1


N, = (3-54)
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, = (3-55)
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 (3-16) are to be employed. The equations to

determine the anisotropic Ks, using Hoenig's displacement equations, are as follows

KI
KII =[B] [A] (3-58)
KIII 2L

where [B] is

mIN1, (cos a + p, sin C)
[B]= Re m2 Nl (cos m+ ,u sin ) (3-59)
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 3-10. Collapsed 20-node 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 (3-57)
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 far-field 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 stress-strain 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 (3-60)

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 11-direction, 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 ) (3-61)
mech
723 Y23
mech
713 Y13
mech
Y12 Y712

or, stated more concisely the total strain is

S mech thermal (3-62)
Total= me + (3-62)









If the thermal expansion is linear with the temperature change, we can rewrite (3-60) 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 near-tip 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 near-tip 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 in-plane and out-of-plane 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 I-III SIFs can be computed. These formulations are different in a few ways.

Computing K with the J-integral 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 (4-1), 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 4-1.

For isotropic materials, the near-tip stresses (in the r, co polar coordinate system)

have been derived by Williams (1957). If the derivative with respect to theta of equation

(4-1) is taken and set to zero, the critical angle, oc, can be determined; equation (4-4).

1 C 3 a) 3
) = COS- KI COS- KIIsin (4-1)
2irr 2 2 2


S= cos- KI 2+sin- +-KI sino-2KIItan- (4-2)
2 r 2 2 2 2
























Figure 4-1. Definition of crack turning angle

o- = 2 r cos [KIsinco+KII(3cos- 1)] (4-3)



) = 2 tan 1l 1 + J (4-4)
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 mixed-mode 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 4-2 in (T,,, r) polar

coordinates. Now T, can be expressed as

T7 = X, sin2 7 + cos2 77 (4-5)


YT





XT XT





YT


Figure 4-2. 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 4-3. 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 4-3. Definition of hoop stress

Mixed mode crack propagation in anisotropic materials is also investigated by

Saouma et al. (1987). The near-tip 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 near-tip 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 Banks-Sills 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 in-plane and anti-plane

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 4-4 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 4-4. In our case, only three toughness values are needed,

however. In Figure 4-4, consider the M(T) specimens that have the K31 and K32 labels.

Here the load is applied in the 3-direction and the crack is propagating in the 1 and 2-

directions, respectively. Since the 3-direction is the rise direction and the 1-2 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 4-5.
































Figure 4-4. Orthotropic toughness values


Crack faces


Figure 4-5. 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 (4-5) we can write a 2D interpolation function of the form


(4-6)









Kk (,k) = K, cos2 (ok) + K, sin2 (ok) (4-7)

We also observe the following trigonometric identities


cos2 tan b -- C-
c b 2 + (4-8)

sin2 tan- b- b22
c b +c

and the components of unit vector a must satisfy the relation

a1 +a2 +a2 =1 (4-9)

Pettit goes on to define equations that can be interpreted as the fracture resistance

components of a in the principal planes. Consider Figure 4-4 once more where the 1-axis

is the loading direction. To estimate an effective K on this plane, we will interpolate

between two toughness values, K12 and K13. Using (4-7) through (4-9) we can write K1

as


K, (a) (K2a22 + K,3a32) (4-10)
1- a2

In a similar fashion we can determine K2 and K3

K2(a)_1 1 2 (K23a32 +K21a12)
K, (a) = -- Ka + Ka
a2 (4-11)
K3 (a)- 1 2 (K31a12 + K32a22)
1 a3
I-a

Finally, Pettit sums the relations in (4-10) through (4-11). 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 ( (4-12)
+ I22 (K23a32 +K2a2)+ -n2(K3,2 + K32a22)
1-a 2 1-a3









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 two-fold. One of the

objectives of this study is to examine the effect of mode mixity on the K-solution, or

stated differently, as the material orientation is changed how does this impact KI-KIII

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

4-6 shows a schematic of the M(T) specimen and Figure 4-7 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,

five-and-a-half inches wide, and an inch-and-half thick. The boxed in area in Figure 4-7

highlights the densely meshed region containing the crack.










i i




2a-







Figure 4-6. Schematic of M(T) specimen


Figure 4-7. 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 4-8, 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, anti-clockwise (of angle zeta) about the x axis, takes place. Table 4-1 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 4-8. Sample transformation

Table 4-1. 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 4-9 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 4-9. 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 3-D 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 K-solution 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 5-1. M(T) models

use though-cracks, 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 5-2.

























a O aCT a

Figure 5-1. Definition of crack inclination angle

To analyze the effect of material orientation of the SIF solution, 17 material

orientation cases are defined; see Figures 5-4 and 5-5. M(T) models are constructed

using ANSYSTM FE software and the flaw is inserted using FRANC3D Next Generation.


Figure 5-2. 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 5-3. M(T) model built in ANSYS

z


x

Figure 5-4. Definition of material orientation




























Figure 5-5. Top view of the cones shown in Figure 5-4

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 5-5.


Portion of External Tank


Figure 5-6. 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 5-6. 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 5-4 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 5-7 is performed. This

can be likened to the knit lines being oriented within the M(T) specimen shown in Figure

5-7. Table 5-1 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 5-7. Definition of case one material orientation for the M(T) FE model









Table 5-1. 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 5-8 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 5-8 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 5-9 through 5-11, the mode I-III SIF is plotted versus the

normalized crack front distance for cases one through four. As expected, when the crack

is inclined KI decreases, and KII-KIII increase. This is consistent with isotropic

assumptions. For the configuration shown in Figure 6-1, the KI and KII are related to p,

or the crack inclination angle (Anderson, 1991)

KI= o cos2 ((P) Q.a
(5-1)
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
- C-ase 2
-e Case 3
-- Case 4


normalized crack front distance


Figure 5-8. Mode I SIF for cases 1-4; 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 5-9. Mode I SIF for cases 1-4; 30 degree crack inclination
















KII vs normalized crack front distance

(p =30 degrees


---Case 0
-A Case 1
- C-Case 2
c-- Case 3
- Case 4


normalized crack front distance


Figure 5-10. Mode II SIF for cases 1-4; 30 degree crack inclination


Kill vs normalized crack front distance


---Case 0
-A-Case 1
C--Case 2
-- Case 3
- Case 4


(P =30 degrees


normalized crack front distance


Figure 5-11. Mode III SIF for cases 1-4; 30 degree crack inclination









The mode I SIF for (p = 0 is listed in tables 5-2 and 5-3. 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 z-axis 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 5-2. KI for cases 0-8, 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 5-3. KI for cases 9-16, 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 near-tip 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 5-12.

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 5-12. 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 5-13. 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 cross-product 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 5-13. 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 y-axis. 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 x-axis. 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 5-13. 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 post-processing of near-tip









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 5-13. If the

inclined crack is oriented in such a fashion so that it rests in the plane of symmetry, the

out-of-plane and in-plane displacements become decoupled and this is a so-called

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 near-tip formulas are utilized.







Knit Line














Figure 5-14. 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-} (5-2)

then the transformation formulas (Ting, 1996) for the elastic constants are

[S,']= [Q][S][QS] (5-2)

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) (5-3
[es]= (5-3)
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 (5-2) and (5-3) 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]{}) (5-4)

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 (5-5)
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 [ ] (5-6)

The two [Q] matrices are related

[Q] = [ l] (5-7)









Equations (5-2) or (5-6) 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 5-15. The near-tip 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 anti-clockwise about the z-axis by co degrees. This is done using standard

transformation equations (Boresi et al 1993).

Y





j ic









Figure 5-15. Cartesian and cylindrical stresses

Calculating co is done via

o)o = aO22 0x 22 + 22 + 2 2/ + 22yz + 2/202"zx + 2a2i2oxy (5-4)

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 5-16. Here we have a close-up 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 5-16. 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 5-2) into a FE model
coordinate system









* Run the model and extract the mode I-III SIFs via FRANC3D

* Rotate material properties into a local crack front coordinate system

* Compute near-tip 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 5-4. 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 5-4 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 K-solution, the measured and predicted angles

are computed at the mid-plane. 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 (4-4). 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 5-4.

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 near-tip 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 5-12)

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 K-solution 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 5-17.

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 mid-plane 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

Banks-Sills et al. (2005) study a similar routine is used to check the M-integral

computations. Recall that the interaction, or M, integral is used to compute the

anisotropic K-solution 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. Banks-Sills et al. use 2-D 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 5-17 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 I-III SIF for case two at both room temperature and with a thermal

gradient is plotted in Figure 5-18. From Table 5-2, the room temperature mid-plane

value for KI is approximately 2.14 psi /in. With the prescribed thermal boundary

conditions, the mid-plane 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 5-18, 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 mixed-mode 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

mid-plane 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 5-17. The next set of results

for this material orientation are for T1 = -200oF and T2 = 1000F.

In Figure 5-19 we see how the larger thermal gradient impacts the mode I SIF.

With these boundary conditions, KI, at the mid-plane, 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

mid-plane reaches -0.122 psi /in and KIII= 0.355 psi /in.

















///////Figure 5-17. Application of thermal loads
( T

Figure 5-17. Application of thermal loads




Full Text

PAGE 1

ANISOTROPIC FRACTURE ANALYSIS OF THE BX-265 FOAM INSULATION MATERIAL UNDER MIXE D-MODE 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 BX-265 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 2-1 Measured material prope rty values from experiment...............................................18 3-1 Various methods for determining the SIF................................................................38 4-1 Direction cosines......................................................................................................66 5-1 Direction cosines for the cas e one material orientation...........................................73 5-2 KI for cases 0-8, 0 degree inclination......................................................................76 5-3 KI for cases 9-16, 0 degree crack inclination...........................................................76 5-4 Measured and predicted tu rning angles (in degrees)................................................83 5-5 Summary of mid-plane KI-KII.................................................................................91 5-6 Summary of mid-plane KIII.....................................................................................91

PAGE 9

ix LIST OF FIGURES Figure page 1-1 Foam loss during Shuttle ascension...........................................................................2 1-2 Spray-on foam insulation (SOFI)...............................................................................3 2-1 Idealized foam microstructure....................................................................................8 2-2 Square specimen tensile test used to determine the Poisson ratio............................15 2-3 Foam test specimens.................................................................................................15 2-4 Summary of stre ss-strain data .................................................................................16 2-5 Coordinate system for the tr ansversely isotropic material.......................................17 2-6 Youngs moduli vs. temperature..............................................................................18 2-7 Shear moduli vs. temperature...................................................................................19 2-8 Coefficient of thermal expansion vs. temperature...................................................19 2-9 Thermal conductivity vs. temperature......................................................................20 2-10 From left to right: M(T), C(T), and SNE(B) fracture test specimens......................22 2-11 Summary K-values for va rious fracture specimens.................................................23 2-12 Clip gauge on C(T) specimen used to measure CMOD...........................................24 2-13 Typical load vs CMOD curve for a brittle material.................................................25 2-14 Divot test panel.........................................................................................................26 2-15 2-D view of divot test...............................................................................................27 2-16 2-D view of test panel after heat flux has been applied...........................................27 2-17 Foam blow out from divot test.................................................................................28 3-1 Coordinate system us ed in equation (3-1)................................................................31

PAGE 10

x 3-2 From left to right: mode I (opening), mode II (shearing), mode III (tearing).........32 3-3 Path of the J-integral................................................................................................39 3-5 Definition of Irwins crack closure integral.............................................................43 3-6 Correlation points typically used for displacement correlation...............................46 3-7 Triangular quarter point elements............................................................................48 3-8 Isoparametric element and degenerated element with mid-side nodes moved to quarter point locations..............................................................................................48 4-1 Definition of crack turning angle.............................................................................57 4-2 Elliptical relationship for T...................................................................................58 4-3 Definition of hoop stress..........................................................................................59 4-4 Orthotropic toughness values...................................................................................62 4-5 Definition of the a and n vectors..............................................................................62 4-7 M(T) model created in ANSYS...............................................................................65 4-8 Sample transformation.............................................................................................66 4-9 Fractured M(T) specimen: the dotted lin e is added to show the location of the knit line.....................................................................................................................67 5-1 Definition of crack inclination angle........................................................................69 5-2 Definition of crack front distance.............................................................................69 5-3 M(T) model built in ANSYS....................................................................................70 5-4 Definition of material orientation.............................................................................70 5-5 Top view of the cones shown in Figure 5-4.............................................................71 5-6 Hypothetical material orie ntation relative to the ET................................................71 5-7 Definition of case one material or ientation for the M(T) FE model........................72 5-8 Mode I SIF for cases 1-4; 0 degree crack inclination..............................................74 5-9 Mode I SIF for cases 1-4; 30 degree crack inclination.............................................74 5-10 Mode II SIF for cases 1-4; 30 degree crack inclination...........................................75

PAGE 11

xi 5-11 Mode III SIF for cases 1-4; 30 degree crack inclination..........................................75 5-12 From left to right: front, left, righ t, and rear sides of the M(T) specimen...............77 5-13 Determination of knit line plane...............................................................................78 5-14 M(T) model with inclined crack...............................................................................79 5-15 Cartesian and cylindrical stresses.............................................................................81 5-16 Crack turning angle..................................................................................................82 5-17 Application of thermal loads....................................................................................87 5-18 Room temperature KI-KIII (F 3D) vs. KI-KIII (DC) with T1 = 0oF and T2 = 100oF........................................................................................................................88 5-19 Room temperature KI (F3D) vs. KI (DC) with T1 = -200oF and T2 = 100oF...........88 5-20 = 30 degrees: room temperature KI (F3D) vs. KI (DC) with T1 = 0oF and T2 = 100oF........................................................................................................................89 5-21 = 30 degrees: room temperature KI (F3D) vs. KI (DC) with T1 = -200oF and T2 = 100oF................................................................................................................90 5-22 Temperature distribut ion along the crack front........................................................92 5-23 KI-KIII 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 BX-265 FOAM INSULATION MATERIAL UNDER MIXE D-MODE 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 closed-cell 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 spray-on 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 re-entry. 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 non-uniform 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 mixed-mode 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 BX-265 foam insulation material.

PAGE 14

1 CHAPTER 1 BACKGROUND Problem Motivation On February 1st, 2003, the Space Shuttle Columbia suffered a catastrophic failure during re-entry. NASA has conduc ted an exhaustive investiga tion of the failure and the consensus now is the breakup was caused by a segment of BX-265 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 re-entry, 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 Figure1-1).

PAGE 15

2 Figure1-1. 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 1-1 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 1-2). Pictures of the Shuttle during liftoff, su ch as those shown in Figure 1-1, 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 1-2. Spray-on 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, crack-lik e, defects for this sort of material is through linear elastic fracture mechanics (LEF M). For isotropic materials, the near-tip stresses are completely characterized by K the stress intensity factor (SIF). The analytical K-solutions for anisotropic materi als take a similar form, but have additional terms that come from the various elastic moduli. The formulations of near-tip stresses for isotropic materials are well established and are now incorporated in many finite element (FE) software packages. Some K-solutions for certain kinds of direction-dependent 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 e-crystal 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 K-soluti 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 crack-like defects in elastic solids and also computes a K-solution 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 K-solutions 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 K-solution is examined 2. The impact of mode mixity on the K-solution 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 K-solution 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 two-dimensional arrays of polygons that fill a planar area and 2) foam s, the three-dimensional 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 2-1 (Gibson and Ashby, 1988) lists several properties of polymer foams and true solids (solid metals and ceramics). Table 2-1. Properties of solids and polymer foams Density (lb/in3) Conductivity (BTU/hr in F)Youngs modulus (psi) Solids 3.6 x 10-2 3.6 x 10-1 4.8 x 10-3 48 3 30 x 106 Polymer foams 3.6 x 10-4 4.8 x 10-3 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 2-1) The equations that arise from this mathematical treatment are in terms of the strut dimensions, which in turn

PAGE 21

8 Figure 2-1: 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 membrane-like 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 (2-1) where t is the cell wall thickness and the l is the edge-length. For open-cell foams 2 2*st C l (2-2) Finally, for the foams of the closed-cell variety we have 3*st C l (2-3) The C terms are numerical constants that de pend on the cell shape, normally taking a value of unity. As mentioned beforehand, certain open-ce 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 (2-4) where E* is the modulus of the foam material. The relative density is related to the dimensions, t and l as shown in (2-2). The second moment of the area, I can also be related to the cell wall dimensions

PAGE 23

10 4 I t (2-5) Now (2-4) can be re-written as 2 1**ssE C E (2-6) Using the same reasoning, we can write down the relationship for the shear moduli 2 2**ssC E (2-7) 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 (2-8) Using (2-8) along with (2-6) and (2-7) we have 1 3 2*1 2 C C C (2-9) The above analysis is for open-celled foam s. When examining closed-cell 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 closed-cell 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 (2-10) Equation (2-10) can be re-wri tten when we note that I t4 e and E* (F/l2)/( /l) which yields 4 4* '' f e st t E Ell (2-11) 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 (2-12) If we substitute (2-12) into (2-11) 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 (2-13) 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 (2-14) 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 closed-cell foam is shown below 2 22*** '1 s ssG CC E (2-15) Finally the last contribution to consider comes from the bending on the cell struts, and that analysis is identical to the ope n-celled formulation presented earlier. We note that the above relations are for foams where the cell walls are equiaxed. Most man-made 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 (2-16)

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 open-celled foams as 2 3 3 1* 2 *(1(1/))E R ER (2-17) For closed cells an additional term 2 1 1(1/)R R appears in (2-17). For the anisotropic shear modulus, we have 31 12* 2 *1G GR (2-18) 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 BX-265 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 all-important 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 un-foamed 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 2-2 and 2-3 display some of the various test specimens used in various tension and torsion tests. A su mmary of typical stress-strain da ta is shown in Figure 2-4. 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 2-2. 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 2-3. Foam test specimens

PAGE 29

16 Figure 2-4. Summary of stress-strain 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 non-linear 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 11-22 plane in Figure 2-5 is a plane of isotropy. Source: personal communication, Doug Wells (MSFC) Summary of Typical Tensile Stress Strain Curves BX-265 / 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 2-5. Coordinate system for th e transversely isotropic material Using the coordinate axes from Figure 2-4 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 2-1 lists the room temperature values of the elastic constants for the orientation shown in Figure 2-5. Since one of the objectives of this study is to analyze the effect of temperature, or thermal loads, on the K-solution 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 (2-19) where the five independent elastic constants are 33 22 11

PAGE 31

18 112233 2313 1331 1221 31; EEE GG EE (2-20) Table 2-1. 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 2-6 through 2-9. 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 -550-450-350-250-150-5050150250 temperature (F)psi E11 E22 E33 Figure 2-6. Youngs moduli vs. temperature Source: personal communication, Doug Wells (MSFC)

PAGE 32

19 Shear moduli0 100 200 300 400 500 600 -500-400-300-200-1000100200300 temperature (F)psi G12 G13 G23 Figure 2-7. 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 -500-400-300-200-1000100200300 temperature (F)in/in/F as11 as22 as33 Figure 2-8. 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 -500-400-300-200-1000100200300 temperature (F)BTU/hr ft F k Figure 2-9. Thermal conductivity vs. temperature Fracture Testing for the BX-265 Foam Material While many foam materials do not exhibit linear-elastic response before yielding, brittle foams normally exhibit linear-elastic beha vior (while in tension) up until fracture takes place. Thus, we would like to use LEFM c oncepts to compute the near-tip 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 (2-21)

PAGE 34

21 where C8 is a constant of proportionality and fs is the failure strength of the struts. Equation (2-21) can also be re-written as 3/2 8'*/cCssKICKI (2-22) where KICs is the toughness of the solid material Equation (2-22) 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 (2-21). 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 well-known 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 single-edge 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, continuum-based relations to determine fracture properties, such as the plane strain toughness. Pictures of these test specimens are shown in Figure 2-10. 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 2-11 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 2-10. 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 toughness-KIc. If the material is brittle and the subsequent load vs crack mouth opening displacement (CMOD) curve (Figure 2-12) is linear-elastic 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 2-13) 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 2-11. Summary K-values 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)(A-S)M(T)(A-S)C(T)(A-S)M(T)(33-S)C(T)(33-S) RT -320F RT -320F RT -320F RT -320F RT -320F 0.1 < a/W < 0.5

PAGE 37

24 (/)Q QP KfaW BW (2-5) Figure 2-12. 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 2-13. 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 defect-free, 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 2-14 are used for the divot testing. Figure 2-14. 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. Two-dimensional sc hematics of this process are shown in Figures 2-15 and 2-16. The applied flux warms the liquid nitrogen and eventually a phase change takes place. The experiment takes place in a thermo-vacuum 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 2-15. 2-D view of divot test In Figure 2-17, 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 2-16. 2-D 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 2-16, by calculating near-tip 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 2-17. 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 well-known 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 al-shaped 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 near-tip 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 near-tip stresses for materials that are anisotropic, or have direction-depe 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 (3-1). 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 (3-1) 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 near-tip stress fields for an isotropic materials tried to take advantage of material symmetry to decouple in-plane, and out-of-plane displaceme nts, which, in turn, makes the mathematical formulati on a little less complicated. Isotropic Fracture Mechanics For isotropic materials, the near-tip 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 far-field loads and are governed by a single parameter, K, the stress intensity factor (Anderson, 1991). (,)()higher order terms 2 ijijK rf r (3-1) Figure 3-1. Coordinate sy stem used in equation (3-1) The coordinate system for equation (3-1 ) is shown in Figure 3-1, 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 near-tip 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 3-2 displays the three modes of deformation. Figure 3-2. 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 r-equation (3-2). 22 x 22 y s zcos -1+2sinsin +1+2cos0 2222 u 1r u=sin +1-2cos-cos -1-2sin0 22 2222 u 00sin 2 KI KII KIII (3-2)

PAGE 46

33 where = 3-4 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 direction-dependent. 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 near-tip stresses for isotropic materials, Sih et al. use a similar approach for materials with direction-dependent 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 (3-3). 11121316 222326 3336 4445 55 6600 00 00 0 0xx yy zz yzyz x zxz x yxySSSS SSS SS SS S symS (3-3)

PAGE 47

34 This configuration is advantageous because the in-plane and out-of-plane 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 (3-4) 31323633 zzxxyyxySSSS (3-5) where the elastic constants in compliance form are defined as iijjS (3-6) 1222112211 1212 1222112211 1212 3 4534411 ReRe0 211 ReRe0 00Rex y zpQpQpQpQ u KI r uqQqQqQqQKII u KIII Q CC (3-7) 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 (3-8) 122622'''jjjqSSS (3-9) Sih et al. enforce a plane strain condition and the reduced compliance matrix becomes

PAGE 48

35 33 33 ij ijijSS SS S (3-10) cossiniiQ (3-11) The C45 and C44 terms in (3-7) are the elastic constants which are in stiffness form iijjC (3-12) where 1 ijijCS (3-13) The aforementioned decoupling of the anti -plane displacement is now evident in equation (3-7). The uz displacement does not depend on KI and KII. This is similar to the isotropic case in equation (3-2); 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 anti-plane 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 (3-14) 313234353633 zzxxyyyzxzxySSSSSS (3-15) 3 1 12 Recossiniijjllj jr umNK (3-16) Here an important distinction is made be tween the general case and Sih et als formulation. Notice in (3-16) how the mode I-III SIFs are to be summed over all three modes of displacement. Thus, all three components of displacement are dependent on KI-KIII. 2 11116121514 2224 2212625 4244 3414645''''' '' ''' '' '''iiii iii ii iii iimSSSSS SS mSSS SS mSSS (3-17)

PAGE 50

37 The i are the roots that occur in conjugate pa irs and they arise from the sixth order polynomial, equation (3-18). The coordinate sy stem used in equations (3-14) and (3-16) is the same as the one used in (3-1). 2 4230 lll (3-18) Where 2 2554544 32 3151456254624 432 41116126626222 222 lSSS lSSSSSS lSSSSSS The N matrix in equation (3-14) is defined as 123 123111ijN (3-20) Where i is defined as l3(i)/ l2(i). There are two situations, also called degenerate cases, where this solution is invalid, however. If anti-plane shear and pl ane strain displacements becomes decoupledfor monoclinic materials, say-Hoenigs 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 N-matrix, equation (3-20), becomes singular because 1 and 2 = i. However, we have near-tip solutions for the isotropic cases. These are shown in equati ons (3-1) and (3-2). 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 3-1. Table 3-1. 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 K-solutions 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. 2-D problems involve discre tizing the line-boundary of the domain and for 3-D 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 near-tip 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 K-solu tion. This particular method of solution would fall under Stage 3 in Table 3-1, and in this case Rices J-contour 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; K-solutions can be generated fo r isotropic and fully anisotropic materials alike. Since many FE models were built and the subsequent K-solutions computed with this software, it would be prudent to discuss the pertinent theo ry (Banks-Sills 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 (3-21) Where = (1-2)/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 (3-22) 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 3-3 Figure 3-3. Path of the J-integral

PAGE 53

40 The J-integral was initially formulated for 2-D problems, however since this material has properties that ar e direction dependent, our analys is requires the use of 3-D FE models. The J-integral, 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 J-integral can be rewritten as ds x q W x u Jj j i ij 1 1 (3-23) 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 J-integral is now evaluated in a cylindrical domain centered over the crack front. Figure 3-4. Cylindrical do main and virtual crack advance for a 3-D crack front Following Figure 3-4 we can define increm ental area of virtual crack advance as max q LAaqsds (3-24) The definition of J in equation (3-22) 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 2-D problems) because they are easier to implement numerically. So, the q term in (3-23) is simply a mathematical tool that enables us to recast the J-integr al into a slightly more usable form. This virtual extension method is the mode rn-day 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 3-D 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 ) ( ) ( ) ( (3-25) where qt(s) is the crack tip perturbation. Rices J-integral 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 J-integral to comput e the Ks individually. Ishikawas method, though, has not been exte nded for anisotropic materials. Yau et al. (1980) uses the M-integral 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 (3-26)

PAGE 55

42 The J-integral can now take the form (1)(2)(1,2)JJJM (3-27) The M(1,2) term is known as the M-inte 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 (3-28) This integral describes the interaction between the two equilibrium states, and sometimes the M-integral 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 3-5. 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 near-tip 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 3-30) 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 (3-29)

PAGE 56

43 Figure 3-5. Definition of Irw ins crack closure integral Now, if the stress and displacement fiel ds are substituted in (3-29), we have 111 2131 Im()Im()Im() 2iijjiijjiijjGJKImNKKIImNKKIIImNK (3-30) When we use the relations for SIFs in (3-26), we get (1)(2)1(1)(2)1 21 (1)(2)1 3()Im()()Im() 1 2 ()Im()iijjiijj iijjKIKImNKKIIKIImNK GJ KIIIKIIImNK (3-31) The terms in (3-27) 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 (3-32)

PAGE 57

44 Now we substitute (3-31) into (3-32) 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 (3-33) (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 (3-34) 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 (3-28) and (3-34). 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 (3-14). Now the two relations for the M-integral 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 (3-35)

PAGE 58

45 From (3-35) 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 (3-36) Finally, we note that (3-36) 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 (3-36) 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 3-2. 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 (3-2) with = and Figure 3-6 for a plan e strain problem dealing with an isotropic material, we have 2 2 1ybyasuu KI r (3-37) 2 2 1xbxasuu KII r (3-38)

PAGE 59

46 where s, again, is the shear modulus determined by E/2(1+ ) and = 3-4 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 3-6. 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 (3-1), 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 mid-side 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 3-7, 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 mid-side nodes moved to th e quarter point locations. Consider the bottom edge of the element on the right in Fi gure 3-8. 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 (3-39) 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 (3-40) We also have shape functions for the mid-side nodes 2 21 11 1,1 (5,7) 2 1 11 1,1 (6,8) 2iii iiiNti Ni (3-41) 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 (3-42)

PAGE 61

48 Figure 3-7. Triangular quarter point elements Figure 3-8. Isoparametric element and dege nerated element with mid-side nodes moved to quarter point locations The global x coordinate is related to th e isoparametric coordinate system via ii x Nx (3-43) where the repeated index is an implied summation. Setting x1 = 0, x2 = L, and x5 = L/4 we have 21 11 24 L xL (3-44)

PAGE 62

49 We can now solve for 21 x L (3-45) 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 (3-46) Again determining the displacements along th e edge containing nodes one, two, and five we write 2 12511 111 22 uuuu (3-47) Now we can substitute (3-45) into (3-47) to obtain 1251 12221224 2 xxxxxx uuuu LLLLLL (3-48) We note that the displacement shown in (3-48) has a x term present. This is consistent with equation (3-2) where the displacements vary with r. We need to differentiate equation (3-48) to obtain the strain in the x direction 12513411424 22xu uuu xLLL xLxLxL (3-49) Finally, it is shown in (3-49) how the strain disp lays the desired singularity (x-1/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 (3-37) and (3-38) are modified because r = L/4 4 2 1ybyasuu KI L (3-50) 4 2 1xbxasuu KII L (3-51) 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 3-D elements. This formulation is particularly useful to the present study si nce 3-D simulations are used to model the fracture response of the foam material. Ingraffea and Manu use a special labeling convention shown in Figure 3-10 for a collapse d 20 node brick element (Figure 3-9), also known as a 15 node wedge. Figure 3-9. 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 (3-52) For the nodes at positions i = 17, 18, 19, 20 and i = 0, i = i = 2111 4ii iN (3-53) At i = 10, 12, 14, 16 and i = i = 0, i = 2111 4ii iN (3-54) Finally, for the nodes at positions i = 9, 11, 13, 15 and i = i = i = 0 2111 4ii iN (3-55) 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 (3-1 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 (3-58) where [B] is 1 1 1 2 1 3cossin []Recossin cossinjjlj jjlj jjljmN BmN mN (3-59)

PAGE 65

52 Figure 3-10. Collapsed 20-node 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 (3-57) 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 far-field 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 stress-stra 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 (3-60) where the total strain is { 11 22 33 23 13 12}T, Th 11(T, Tref) is the free thermal strain in the 11-direction, 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 (3-61) or, stated more concisely the total strain is total = mech + thermal (3-62)

PAGE 67

54 If the thermal expansion is linear with the te mperature change, we can rewrite (3-60) 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 (3-63) 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 near-tip 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 near-tip 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 in-plane and out-o f-plane 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 I-III SIFs can be computed. These fo rmulations are different in a few ways. Computing K with the J-integral 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 (4-1), 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 4-1. For isotropic materials, the ne ar-tip stresses (in the r, polar coordinate system) have been derived by Williams (1957). If the derivative with respect to theta of equation (4-1) is taken and set to zero, the critical angle, c, can be determined; equation (4-4). 213 coscossin 222 2 KIKII r (4-1) 213 cos1sinsin2tan 2222 2rKIKIIKII r (4-2)

PAGE 70

57 Figure 4-1. Definition of crack turning angle 1 cossin3cos1 2 22rKIKII r (4-3) 2 1118() 2tan 4()cKIIKI KIIKI (4-4) 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 mixed-mode 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 4-2 in (T, ) polar coordinates. Now T can be expressed as 22sincosTTTXY (4-5) Figure 4-2. 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 4-3. 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 4-3. Definition of hoop stress Mixed mode crack propagation in anisotr opic materials is also investigated by Saouma et al. (1987). The near-tip 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 near-tip 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 Banks-Sills 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 in-p lane and anti-plane 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 4-4 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 4-4, consider th e M(T) specimens that have the K31 and K32 labels. Here the load is applied in the 3-directi on and the crack is propa gating in the 1 and 2directions, respectively. Since the 3-direction is the rise di rection and the 1-2 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 4-5.

PAGE 75

62 Figure 4-4. Orthot ropic toughness values Figure 4-5. 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 (4-6) The subscript in the theta term denotes the axis normal to the principal plane. Using equation (4-5) we can write a 2D in terpolation function of the form

PAGE 76

63 22cossinkkkikkjkKKK (4-7) We also observe the following trigonometric identities 2 21 22 2 21 22costan sintan bc cbc bb cbc (4-8) and the components of unit vector a must satisfy the relation 222 1231 aaa (4-9) Pettit goes on to define equa tions that can be interpreted as the fracture resistance components of a in the principal planes. Consider Figure 4-4 once more where the 1-axis is the loading direction. To estimate an effective K on this plane, we will interpolate between two toughness values, K12 and K13. Using (4-7) through (4-9) we can write K1 as 22 1122133 2 11 () 1 KaKaKa a (4-10) 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 (4-11) Finally, Pettit sums the relations in (4-10) th rough (4-11). 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 (4-12)

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 two-fold. One of the objectives of this study is to examine the e ffect of mode mixity on the K-solution, or stated differently, as the material orientat ion is changed how does this impact KI-KIII 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 4-6 shows a schematic of the M(T) specime n and Figure 4-7 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, five-and-a-half inches wide, and an inch-and-ha lf thick. The boxed in area in Figure 4-7 highlights the densely meshed region containing the crack.

PAGE 78

65 Figure 4-6. Schematic of M(T) specimen Figure 4-7. 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 4-8, 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, anti-clockwise (of angle zeta) about the x axis, takes place. Table 4-1 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 4-8. Sample transformation Table 4-1. 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 4-9 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 4-9. 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 3-D 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 K-solution 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 5-1. M(T) models use though-cracks, 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 5-2.

PAGE 82

69 Figure 5-1. 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 5-5. M(T) models are constructed using ANSYS FE software and the flaw is in serted using FRANC3D Next Generation. Figure 5-2. 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 5-3. M(T) model built in ANSYS Figure 5-4. 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 5-5. Top view of the cones shown in Figure 5-4 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 5-5. Figure 5-6. 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 5-6. 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 z-axis (local rise direction) will be for each material case. Let us consider an example. Point one in Figure 5-4 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 5-7 is performed. This can be likened to the knit lines being orient ed within the M(T) specimen shown in Figure 5-7. Table 5-1 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 5-7. Definition of case one mate rial orientation for the M(T) FE model

PAGE 86

73 Table 5-1. 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 5-8 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 5-8 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 5-9 through 5-11, the mode I-III 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 KII-KIII increase. This is consistent with isotropic assumptions. For the configur ation shown in Figure 6-1, th e KI and KII are related to or the crack inclination angle (Anderson, 1991) 2cos sincos KIa KIIa (5-1) We can see that KII is zero when = 0 and that KI will decrease for increasing values of

PAGE 87

74 Figure 5-8. Mode I SIF for cases 1-4; 0 degree crack inclination Figure 5-9. Mode I SIF for cases 1-4; 30 degree cr ack inclination

PAGE 88

75 Figure 5-10. Mode II SIF for cases 1-4; 30 degree crack inclination Figure 5-11. Mode III SIF for cases 1-4; 30 degree crack inclination

PAGE 89

76 The mode I SIF for = 0 is listed in tables 5-2 and 5-3. 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 z-axis 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 5-2. 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 5-3. KI for cases 9-16, 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 near-tip 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 5-12. 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 5-12. 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 5-13. 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 oss-product 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 5-13. 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 y-axis. 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 x-axis. 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 5-13. 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 post-processing of near-tip

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 5-13. If the inclined crack is oriented in such a fashion so that it rests in the plane of symmetry, the out-of-plane and in-plane displacements become decoupled and this is a so-called 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 near-tip formulas are utilized. Figure 5-14. 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 (5-2) then the transformation formulas (Ti ng, 1996) for the elastic constants are 'T SSSQSQ (5-2) where S is the matrix of elastic constants in the rotated system and 222 111111111 222 222222222 222 333333333 232323232323232223 313131313131313131 121212121212121212222 222 222SQ (5-3) The order of the terms in and (5-2) 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 (5-4) there is an alternate definition of the [QS] matrix 222 111111111 222 222222222 222 333333333 232323232323232223 313131313131313131 121212121212121212222 222 222CQ (5-5) The transformed elastic constants in sti ffness form are therefore determined by '[][][]T CCCQCQ (5-6) The two [Q] matrices are related 1T SCQQ (5-7)

PAGE 94

81 Equations (5-2) or (5-6) 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 5-15. The near-tip 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 anti-clockwise about the z-axis by degrees. This is done using standard transformation equations (Boresi et al 1993). Figure 5-15. Cartesian and cylindrical stresses Calculating is done via 222 222222222222 x xyyzzyzzxxy (5-4) 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 5-16. Here we have a close-up 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 5-16. 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 5-2) into a FE model coordinate system

PAGE 96

83 Run the model and extract the mode I-III SIFs via FRANC3D Rotate material properties into a local crack front coordinate system Compute near-tip 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 5-4. 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 5-4 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 K-solu tion, the measured and predicted angles are computed at the mid-plane. 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 (4-4). 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 5-4. 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 near-tip 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 5-12) 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 K-solution 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 5-17. 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 mid-pl 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 Banks-Sills et al. (2005) study a similar routine is used to check the M-integral computations. Recall that the interaction, or M, integral is used to compute the anisotropic K-solution 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. Banks-Sills et al. use 2-D 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 5-17 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 I-III SIF for case two at both room temperature and with a thermal gradient is plotted in Figure 5-18. From Table 5-2, the room temperature mid-plane value for KI is approximately 2.14 psi in. With the prescribed thermal boundary conditions, the mid-plane 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 5-18, 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 mixed-mode 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 mid-plane 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 5-17. The next set of results for this material orientation are for T1 = -200oF and T2 = 100oF. In Figure 5-19 we see how the larger ther mal gradient impacts the mode I SIF. With these boundary conditions, KI, at the mid-plane, 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 mid-plane reaches -0.122 psi in and KIII = 0.355 psi in. Figure 5-17. Applicatio n of thermal loads

PAGE 101

88 KI-KIII vs normalized crack front distance = 0 degrees-0.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 5-18. Room temperature KI-KIII (F3D) vs. KI-KIII (DC) with T1 = 0oF and T2 = 100oF KI-KIII vs normalized crack front distance = 0 degrees-1.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 5-19. 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 I-III 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. KI-KIII vs normalized crack front distance = 30 degrees-0.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 5-20. = 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 mid-plane in Figure 5-20, 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 mid-plane, 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 mid-plane is approximately 1.74 psi in, and KII and KIII are -0.455 psi in 0.0394 psi in, respectively. KI-KIII vs normalized crack front distance = 30 degrees-1.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 5-21. = 30 degrees: room temperature KI (F3D) vs. KI (DC) with T1 = -200oF and T2 = 100oF A final summary of the mid-plane 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 mid-plane. The units for KI-KIII are psi in.

PAGE 104

91 Table 5-5. Summary of mid-plane KI-KII 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 5-6. Summary of mid-plane 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 re-use 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 5-22. Here the temperature varies appr eciably. To analyze this problem, the elastic constants are evaluated usi ng Figures 2-5 and 2-6. With those values the sixth order characteristic equation is solv ed at each point, and those roots are used in equation in (3-58) 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 5-22. Temperature dist ribution along the crack front KI-KIII vs normalized crack front distance = 30 degrees-1.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 5-23. KI-KIII vs. normalized crack front distance; TG2 applied in thickness direction

PAGE 106

93 The results of this model are shown in Fi gure 5-23. 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 mid-plane, KI = 1.92 psi in. Comparing this to the mid-plane value in Figure 5-21, KI = 1.74 psi in; the value increased by approximately 10%. Summary The effect of material orientation on th e anisotropic K-solution 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 (5-1) we see that KI should decrease while the remaining modes will increase as the gets larger. This trend is evid ent in Figures 5-9 through 5-11. 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. Mid-plane 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 anti-plane 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 near-tip 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 K-solu tion for a variety of material orientations: isotropy, m onoclinic, and general anisotropy. Alternatively, using the DC method, another approach to determining the anisotropic K-solution 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 K-solution 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 BX-265 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 (A-1) 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 (A-2) Finally we have the compatibility equations 22 2 22 yxyyz xxz y xxyyx (A-3) 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 (A-4) If we substitute (A-4) into the compatibility relations we have a system of equations in terms of U and F 43 320 0LULF LULF (A-5) Where 44444 4222612661611 432234 3333 3242546145615 3223 222 2554544 22'2'2''2'' '''''' '2''LSSSSSS x xyxyxyy LSSSSSS xxyxyy LSSS xxyy (A-6) Through equation (A-5) we can eliminate one of the stress functions, F for example, and it can be shown that 2 4230LLLU (A-7) 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 (A-8) zxy (A-9)

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 (A-10) 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 (A-11) with first and second derivatives 2' '' x x y eye (A-12) If we substitute the derivates into equations (A-12) we get 20ab (A-13) 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 (A-7) which results in 2 243lll (A-14) where 2 2554545 32 3151456254624 432 4111612662622'2'' '''''' '2'2''2'' lSSS lSSSSSS lSSSSSS (A-15)

PAGE 114

101 The coefficients in (A-14) 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 out-of-plane and in-plane 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 (A-16) 2 554544'2''0SSS (A-17) 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 near-tip 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 (A-18) 31323633zzxxyyxySSSA (A-19) 2 111612'''jjj p SSS (A-20) 122622'''jjjqSSS (A-21) cossiniiQ (A-22)

PAGE 115

102 1222112211 1212 1 21222112211 1212 3 3 4534411 ReRe 211 ReRe Re KIpQpQKIIpQpQ u r uKIqQqQKIIqQqQ u Q KIII CC (A-23) where C45 and C44 are the elastic constants in stiffness form C (A-24) where 1 ijijCS (A-25)

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, 120-121, 430-431. ANSYS Elements Reference, ANSYS Rel ease 5.6; ANSYS, Inc. November 1999. Banks-Sills L., Hershkovitz I., Wawrzynek P., Eliasi R., and Ingraffea T., Methods for Calculating Stress Intensity Factors in An isotropic Materials: Part 1-Z = 0 is a Symmetric Plane, Engineering Fracture Mechanics, 2005, Vol. 72, pp. 2328-2358. Barnett D.M., and Asaro R.J., The Fracture Mechanics of Slit-Like Cracks in Anisotropic Elastic Media, Journal of th e Mechanics and Physics of Solids, 1972, Vol. 20, pp. 353-366. Barsom, and Rolfe S. Fracture and Fatigue Control in Structures: Applications of Fracture Mechanics, Woburn, MA: Butterworth-Hei 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. 25-37. Bowie O.L., and Freese C.E., Central Crack in Plane Orthotropic Rectangular Sheet, Engineering Fracture Mechanics, 1972, Vol. 8, pp. 49-58. Bogy D.B., The Plane Solution for Anisot ropic Elastic Wedges Under Normal and Shear Traction, Journal of Applied Mechanics, 1972,Vol. 39, pp. 1103-1109. 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. 185-201. 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, 1061-1065.

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. 544-553. 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. 2033-2058. 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. 1-17. Cherepanov G.P. Mechanics of Brittle Fracture, New York: McGraw-Hill, 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., Mixed-Mode K-Calculations fo r Anisotropic Materials, Engineering Fracture Mechanics, 2002, Vol. 69, pp. 909-922. Dieter G.E. Mechanical Metallurgy, New York: McGraw-Hi ll Inc., 1976, pp. 55-64. 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. 519-527. Fowlkes C.W., Fracture Toughness Tests of a Rigid Polyurethane Foam, International Journal of Fracture, 1974, Vol. 10, pp. 99-108. Freund L.B., Stress Intensity Factor Calcul ations Based on a Conservation Integral, International Journal of Solids and Structures 1978, Vol. 14, pp. 241-250. Gent A.N., and Thomas A.G., The Deforma tion of Foamed Elastic Materials, Journal of Applied Polymer Science, 1959, Vol. 1, pp. 107-113. Gent A.N., and Thomas A.G., Mechanics of Foamed Elastic Materials, Rubber and Chemistry, 1963, Vol. 36, pp. 597-610. Gibson L.J., and Ashby M., The Mechanics of Three-Dimensional Cellular Solids, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 1982, Vol. A 382, pp. 43-59. Gibson L.J., and Ashby, M. Cellular Solids, Elmsford, NY: Pergamon Press Inc., 1988, pp. 76-81, 241-256. 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. 393-403. Huang J.S., and Gibson L.J., Fracture Toughness of Brittle Honeycombs, Acta Metallurgica Materialia, 1991, Vol. 39, pp. 1617-1626. 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. 2-28. Hyer M.W. Stress Analysis of Fiber-Rein forced Composite Materials, New York: McGraw-Hill Inc., 1998, pp. 69-75. Ingraffea A.R., and Manu C., Stress Intensity Factor Computation in Three Dimensions with Quarter-Point Elements, International Journal for Numerical Methods in Engineering, 1980, Vol. 15, 1427-1445. 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. 361-364. 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. 243-246. Ko W.L., Deformations of Foamed Elasto mers, Journal of Cellular Plastics, 1965, Vol.1, pp. 45-50. 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. 197-203. Kreyszig E. Advanced Engineering Mathematics, New York: John Wiley and Sons Inc., 1962, pp. 105-105, 136. Lekhnitskii S.G. Theory of Elasticity of an Anisotropic Body, San Francisco: HoldenDay, Inc., 1963, pp. 3, 104-110, 117-119, 122-123. 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. 413-414, 433-434.

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. 65-86. 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. 111-128. Maiti, S.K., Ashby M.F., and Gibson L.J., Fracture Toughness of Brittle Cellular Solids, Scripta Mettalurgica, 1984, Vol. 18, pp. 213-217. 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. 77-85. McIntyre A., and Anderton G.E., Fracture Properties of a Rigid Polyurethane Foam Over a Range of Densities, Polymer, 1979, Vol. 20, pp. 247-253. Moran B., and Shih C.F., A General Treat ment of Crack-Tip Contour Integrals, International Journal of Fracture, 1987, Vol. 35, pp. 295-310. 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. 37-42. Nobile L., and Carloni C., Fr acture Analysis for Orthotropic Cracked Plates, Journal of Composite Structures, 2005, Vol. 68, pp. 285-293. Patel M.R., and Finnie I., Str uctural Features and Mechan ical Properties of Rigid Cellular Plastics, Journal of Materials, 1970, Vol.5, pp. 909-932. 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. 379-386. 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. 87-93. 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. 115-121. Sih G.C., Paris P.C., and Irwin G.R., On Cr acks in Rectilinearly Anisotropic Bodies, Engineering Fracture Mechanics, 1965, Vo1. 1, pp. 189-203. Sih G.C., Strain Energy Density Factor Applied to Mixed Mode Crack Problems, International Journal of Fracture, 1974, Vol. 10, pp. 305-321.

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. A-187, pp. 229-260. 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. 54-55, 118-133. Timoshenko S.P., and Goodier J.N. Theory of Elasticity, Tokyo: McGraw-Hill, 1982, pp. 354-380. Tupholme G.E., A Study of Cracks in Orthotro pic Crystals Using Di slocation Layers, Journal of Engineering Mathematics, 1974, Vol. 8, pp. 57-69. Wells A., Unstable Crack Propagation in Metals-Cleavage 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. 49-53. Williams M.L., On the Stress Distribution at th e Base of a Stationary Crack, Journal of Applied Mechanics, 1957, Vol. 24, pp. 109-114. Yau J.F., Wang S.S., and Corten H.T., A Mixed-Mode Crack Anal ysis of Isotropic Solids Using Conservation Laws of Elasti city, Transactions of the ASME, 1980, Vol. 47, pp. 335-341.

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 full-time teacher at the college level or perhaps retain a positi on in the aerospace industry.