UFDC Home  myUFDC Home  Help 



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