Citation |

- Permanent Link:
- https://ufdc.ufl.edu/UFE0021260/00001
## Material Information- Title:
- Prediction of Gas Permeability in Composite Laminates Using Three-Dimensional Finite Elements
- Creator:
- Xu, Jianlong
- Place of Publication:
- [Gainesville, Fla.]
- Publisher:
- University of Florida
- Publication Date:
- 2007
- Language:
- english
- Physical Description:
- 1 online resource (105 p.)
## Thesis/Dissertation Information- Degree:
- Doctorate ( Ph.D.)
- Degree Grantor:
- University of Florida
- Degree Disciplines:
- Mechanical Engineering
Mechanical and Aerospace Engineering - Committee Chair:
- Sankar, Bhavani V.
- Committee Members:
- Cristescu, Nicolae
Ifju, Peter Mecholsky, John J. - Graduation Date:
- 8/11/2007
## Subjects- Subjects / Keywords:
- Composite materials ( jstor )
Delamination ( jstor ) Laminates ( jstor ) Material properties ( jstor ) Microcracks ( jstor ) Modeling ( jstor ) Resultants ( jstor ) Stitches ( jstor ) Stress ratio ( jstor ) Two dimensional modeling ( jstor ) Mechanical and Aerospace Engineering -- Dissertations, Academic -- UF composite, fem, permeability City of Gainesville ( local ) - Genre:
- Electronic Thesis or Dissertation
born-digital ( sobekcm ) Mechanical Engineering thesis, Ph.D.
## Notes- Abstract:
- Gas permeability in composite laminates is investigated based on Darcy's law for porous materials. The composite with transverse cracks and delaminations is treated as a porous material. The permeability is derived in terms of crack densities in each ply, interfacial intersection area and a material constant to be determined by experiments. The permeability in cross-ply laminates is investigated first. A three-dimensional finite element model in conjunction with the concept of strain energy release rate is used to determine the crack densities in each ply of the laminate. Then a three-dimensional model with delaminations is used to obtain the intersection area that forms the leakage path. Parametric studies were performed, and it has been found that the permeability is related not only to the resultant forces and thermal loads but also to the delamination length, the delamination shape, the stacking sequence and the temperature dependent material properties. Finally, the progressive permeability in cross-ply laminate is predicted as function of applied loads. Multi-directional laminates are also investigated. Unlike the transverse matrix cracks found in cross-ply laminates, so-called stitch cracks were observed in angle-plies in multi-directional laminates. A unit cell is taken based on experimental observation and the evolution of the stitch crack is studied based on the strain energy release rate. The stitch crack length is predicted as a function of applied strain under uniaxial loading. The effects of ply angle and ply thickness on the formation and propagation of stitch crack have been investigated. Some experimentally observed phenomena are confirmed by the modeling results. A special cracking scenario in laminate [0/60/90]s is considered for the permeability prediction including progressive changes in permeability. ( en )
- General Note:
- In the series University of Florida Digital Collections.
- General Note:
- Includes vita.
- Bibliography:
- Includes bibliographical references.
- Source of Description:
- Description based on online resource; title from PDF title page.
- Source of Description:
- This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
- Thesis:
- Thesis (Ph.D.)--University of Florida, 2007.
- Local:
- Adviser: Sankar, Bhavani V.
- Electronic Access:
- RESTRICTED TO UF STUDENTS, STAFF, FACULTY, AND ON-CAMPUS USE UNTIL 2009-08-31
- Statement of Responsibility:
- by Jianlong Xu.
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- Copyright Xu, Jianlong. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
- Embargo Date:
- 8/31/2009
- Resource Identifier:
- 439096401 ( OCLC )
- Classification:
- LD1780 2007 ( lcc )
## UFDC Membership |

Downloads |

## This item has the following downloads: |

Full Text |

PAGE 1 1 PREDICTION OF GAS PERMEABILITY IN COMPOSITE LAMINATES USING THREE-DIMENSIONAL FINITE ELEMENTS By JIANLONG XU A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2007 PAGE 2 2 2007 Jianlong Xu PAGE 3 3 To my mom and to my wife. PAGE 4 4 ACKNOWLEDGMENTS First, I would like to thank my family, especi ally my wife Yanzhi Wang, for their support and love through my life and my graduate career I would like to thank my colleagues at the Center for Advanced Composites for their effort to create a friendly working environment. I also would like to thank my committee members, Dr. John Jr. Mecholsky, Dr. Nicolae D. Cristescu, and Dr. Peter G. Ifju for evaluating my work. Especially, I would like to thank my committee chair, Dr. Bhavani V. Sankar, for both his acad emic and financial support. Without his guidance and mentoring, I could not have finished th is work. His profound knowledge, creative thinking, and positive attitude always motivated me thr ough out my research ca reer. Finally, I thank NASA for the financial support under the H ydrogen Research and Education Program. PAGE 5 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS...............................................................................................................4 LIST OF TABLES................................................................................................................. ..........7 LIST OF FIGURES................................................................................................................ .........8 ABSTRACT....................................................................................................................... ............11 CHAPTER 1 INTRODUCTION..................................................................................................................13 Background..................................................................................................................... ........13 Literature Review.............................................................................................................. .....15 Objectives..................................................................................................................... ..........19 2 MODELS AND METHODS..................................................................................................22 Permeation Model............................................................................................................... ....22 Preliminary Two-Dimensional FE Modeling.........................................................................25 3 CRACK DENSITY PREDICTION........................................................................................31 Methodology.................................................................................................................... .......31 Predicted Results.............................................................................................................. ......34 Uniaxial Loading.............................................................................................................34 Biaxial Loading...............................................................................................................36 4 INVESTIGATION ON INTERSECTION AREA.................................................................45 Model Description.............................................................................................................. ....45 Intersection Area Computation...............................................................................................46 Effects of Delamination Shap e and Delamination Length.....................................................47 Single Delamination Model....................................................................................................48 Effect of Lay-Ups.............................................................................................................. .....49 Temperature Dependent Properties........................................................................................50 Progressive Permeability Prediction.......................................................................................51 PAGE 6 6 5 DAMAGE AND PERMEABILITY PRED ICTION IN MULTI-DIRECTIONAL LAMINATES...................................................................................................................... ...62 Damage Prediction.............................................................................................................. ....62 Effects of Stitch Crack Length a .....................................................................................63 Effects of Spacing lx........................................................................................................66 Effects of Spacing ly........................................................................................................66 Permeability Prediction Model...............................................................................................67 6 CONCLUSION AND DISCUSSION....................................................................................79 Conclusions.................................................................................................................... .........79 General Recommendations.....................................................................................................80 Future Work.................................................................................................................... ........80 APPENDIX A SOME FLUID FLOW LAWS RE LATED TO PERMEABILITY........................................82 Darcys Law.................................................................................................................... ........82 Poiseuille Law................................................................................................................. .......82 Darcy-Weisbach Equation......................................................................................................82 B MATLAB CODES FOR CRACK DENSITY PREDICTION...............................................84 Matlab Codes................................................................................................................... .......84 Least Square Approximation for [A] Matrix..........................................................................90 C MATLAB CODES TO PREDICT PR OGRESSIVE PERMEABILITY...............................92 D SUPERPOSITION OF STRAIN ENERGY RELEASE RATE.............................................99 LIST OF REFERENCES.............................................................................................................102 BIOGRAPHICAL SKETCH.......................................................................................................105 PAGE 7 7 LIST OF TABLES Table page 2-1 Orthotropic material properties for the epoxy/graphite composite material in the present study. (Bapanap alli, et al. 2006)............................................................................27 4-1 Comparison between computed and tested permeability...................................................53 4-2 Typical intersection areas for [0/90n/0] laminates. (unit: m2)............................................53 5-1 Intersection areas for unit cells with different numbers of stitch cracks...........................70 PAGE 8 8 LIST OF FIGURES Figure page 1-1 Lockheed Martin X-33 (http://www.makelengin eering.com/dir/Company/ History/History.htm).........................20 1-2 Microcracks in the 90o-ply and interfacial delaminations between 0oand 90oplies due to thermal/mechanical loading....................................................................................20 1-3 Gas permeation pathway provided by microc racks and delaminations in cross-ply laminates...................................................................................................................... ......21 2-1 Pressure distribution along thic kness direction of laminate..............................................28 2-2 Laminar flow in a rectangle tube.......................................................................................28 2-3 Two-dimensional delamination model. A) Delamination from surface plies [90/0/90]. B) Delamination from middle ply [0/90/0].......................................................29 2-4 J-integral variation with respect to normalized delamination length a/h A) [90/0/90] laminate. B) [0/90/0] laminate...........................................................................................29 2-5 COD variation with respect to crack dens ity. A) [90/0/90] laminate. B) [0/90/0] laminate....................................................................................................................... .......30 3-1 Two microcracking scenarios in composite laminate [0/ 90/0] under uniaxial loading. A) Unit cell for microcracks in middle ply and new crack formation. B) Unit cell for microcracks in surface ply and new cracks fo rmation. (Courtesy: Bapanapalli, et al., 2006) ......................................................................................................................... .....37 3-2 Three-dimensional unit cell with microc racks in both middle and surface plies..............37 3-3 Response surface for [A] matrix of laminate [0/90/0]. A) A11. B) A12. C) A22.................38 3-4 Crack density evolution in the middle ply under loading in x -direction............................39 3-5 Crack density evolution in th e surface ply under loading in y -direction...........................39 3-6 Crack density prediction for laminate [0/90/0]. Nx/ Ny = 0.5..............................................40 3-7 Crack density prediction for laminate [0/902/0]. Nx/ Ny = 0.5............................................40 3-8 Crack density prediction for laminate [0/904/0]. Nx/ Ny = 0.5............................................41 3-9 Crack density prediction for laminate [0/90/0]. Nx/ Ny = 1.0..............................................41 3-10 Crack density prediction for laminate [0/902/0]. Nx/ Ny = 1.0............................................42 PAGE 9 9 3-11 Crack density prediction for laminate [0/904/0]. Nx/ Ny = 1.0............................................42 3-12 Crack density prediction for laminate [0/90/0]. Nx/ Ny = 2.0..............................................43 3-13 Crack density prediction for laminate [0/902/0]. Nx/ Ny = 2.0............................................43 3-14 Crack density prediction for laminate [0/904/0]. Nx/ Ny = 2.0............................................44 4-1 Three-dimensional delamination mo del and typical deformation.....................................53 4-2 COD along specific path s hows the local effect................................................................54 4-3 Variation of intersecti on area with respect to Nx, Ny, and T. ...........................................54 4-4 J-integral distribution along crack fr ont with different connect radius..............................55 4-5 Effect of connecting radius of delamination on intersection area......................................55 4-6 Effect of delamination length on intersection area............................................................56 4-7 Typical deformation of unit cell from different views.......................................................56 4-8 Gas leakage path formed by the de lamination from surface ply only...............................57 4-9 Comparison between two different scenarios....................................................................57 4-10 Geometric model and typical deforma tion for multiple-layer laminates...........................58 4-11 Effect of temperature dependent materi al properties on the intersection area...................58 4-12 Unit cell used to compute typical intersection area. (x = y = 2 cm-1).............................59 4-13 Typical deformation of the unit cell of [0/90/0]................................................................60 4-14 Progressive permeability in laminates with loading ratio = 0.5.....................................60 4-15 Progressive permeability in laminates with loading ratio = 1.0.....................................61 4-16 Progressive permeability in laminates with loading ratio = 2.0.....................................61 5-1 Unit cell taken from laminate [0/ /90]s. A) Top view of laminate. B) Top view of the unit cell. C) Three-dimensional view of the half of the unit cell.......................................71 5-2 Comparison of strain energy release rate computed from two different methods.............71 5-3 Thermal stress effect due to difference in curing temperature and service temperature. A) [0/60/90]s. B) [0/30/90]s...............................................................................................72 5-4 Effect of angle on the formation of stitch crack.............................................................73 PAGE 10 10 5-5 Progressive stitch crack growth in terms of strain in [0/ /90]s laminates..........................73 5-6 Effect of number of angle plies, n, of lay-up [0/60n/90]s on the strain energy release rate. ......................................................................................................................... .....74 5-7 Effect of number of angle plies, n, of lay up [0/30n/90]s on the strain energy release rate. ......................................................................................................................... .....74 5-8 Effect of spacing lx on the strain ener gy release rate of stitch crack.................................75 5-9 Strain energy release rate related to the formation of new stitch crack.............................75 5-10 Interlaminar cross-sectional area formed by crack openings of stitch crack and transverse crack............................................................................................................... ...76 5-11 Special cracking scenarios in multi-directional laminates.................................................76 5-12 Unit cell for laminate with stitch cracks............................................................................77 5-13 Effect of the stitch crack le ngth on the inters ection areas.................................................77 5-14 Permeability in [0/60/90]s laminate with stitch cracks versus average stress compared to the permeability in [0/904/0] laminate...........................................................................78 5-15 Permeability in [0/60/90]s laminate with stitch cracks versus applied strain compared to the permeability in [0/904/0] laminate...........................................................................78 D-1 Angled crack subjected to biaxial loading.......................................................................101 PAGE 11 11 Abstract of Dissertation Pres ented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy PREDICTION OF GAS PERMEABILITY IN COMPOSITE LAMINATES USING THREE-DIMENSIONAL FINITE ELEMENTS By Jianlong Xu August 2007 Chair: Bhavani V. Sankar Major: Mechanical Engineering Gas permeability in composite laminates is investigated based on Darcys law for porous materials. The composite with transverse cr acks and delaminations is treated as a porous material. The permeability is derived in terms of crack densities in each ply, interfacial intersection area and a material c onstant to be determined by experiments. The permeability in cross-ply laminates is investigated first. A thr ee-dimensional finite element model in conjunction with the concept of strain energy release rate is used to determine the crack densities in each ply of the laminate. Then a three-dimensional mode l with delaminations is used to obtain the intersection area that forms the leakage path. Para metric studies were performed, and it has been found that the permeability is relate d not only to the resultant forces and thermal loads but also to the delamination length, the delamination shap e, the stacking sequence and the temperature dependent material properties. Finally, the progressive permeabil ity in cross-ply laminate is predicted as function of applied loads. Multi-directional laminates are also investigated. Unlike the transverse matrix cracks found in cross-ply laminates, so-called stitch cracks were observed in angle-plies in multi-directional laminates. A unit cell is taken based on experime ntal observation and the evolution of the stitch crack is studied based on the stra in energy release rate. The stitch crack length is predicted as a PAGE 12 12 function of applied strain under uniaxial loading. The effects of ply angle and ply thickness on the formation and propagation of stitch crack ha ve been investigated. Some experimentally observed phenomena are confirmed by the modeli ng results. A special cracking scenario in laminate [0/60/90]s is considered for the permeability pr ediction including progressive changes in permeability. PAGE 13 13 CHAPTER 1 INTRODUCTION Background High strength/weight and stiffn ess/weight ratio makes grap hite/epoxy composite materials much more favorable to other conventional ma terials in aerospace structures. Reduction of structural weight of space vehicl es is a critical issue to reduce the cost of sending payload to space. It is estimated that various gas storage tanks account for about 50% of the dry weight of space vehicles. Currently, the cost to transp ort payload into orbit is about $10,000 per pound. NASA initiated the next generation single-stage-to-o rbit reusable launch vehicle (RLV) program trying to reduce the cost to $1,000 per pound. One su ch vehicle is the Lockheed Martin X-33 (Figure 1-1). The linear aerospike engines in X-33 opera te using hydrogen mixed with an oxidizer, which needs to be stored at cryogenic temperatur es. It had been estimated that composite fuel tank would reduce its weight by 40% and the overall weight of the vehicle by 14% compared to the conventional metallic fuel tank (Kessler, Matuszeski, & Mc Manus, 2001). The honeycomb sandwich structure was used in liquid hydroge n (LH2) tank in X-33. In November 1999, the outer wall of the sandwich structure explosiv ely disbonded when the tank was ground-tested at the NASA Marshall Space Flight Center. The cause of the failure was believed to be excessive permeation of liquid hydrogen from inner face sh eet to honeycomb. When the tank was warmed up, the gas in honeycomb could not escape quic kly, and the corresponding pressure increase caused the failure. The failure has caused a setback to further app lication of composite laminate in cryogenic tanks. Hence, understanding the mechanism of gas leakage through composite laminates is very critical to the future application of compos ite materials in cryogenic storage tanks. PAGE 14 14 Composite laminate usually include several layers with different fiber directions. Generally, the cryogenic tank is under high mechanical/thermal loading and also subjected to cyclic loads, which could cause damage in the laminate such as microcracks and delaminations. From experimental observation (Berthelot, 2003), transverse mi crocracks developed in the 90o-ply of cross-ply laminate when the laminate was loaded in the 0o-direction. The microcracks in the present work refers to those matrix cracks fully developed in the thickness direction of the ply and fully extended through the laminate. Th ese microcracks induced local stress concentrations at the cr ack tip and can cause the interl aminar delaminations between 0oand 90o-plies (Figure 1-2). In reality, the composite laminate is under bi-a xial cyclic loading as the wall of the cryogen tank. In this case, the matrix mi crocracks will form in both 90o-plies and 0o-plies of the cross-ply laminates. Because of the interlaminar delami nation, microcracks in different plies can be connected by intersection areas (s haded areas) as shown in Figure 1-3. The network formed by microcracks and intersection areas provide a pa thway for hydrogen to leak through the laminate. Of course, there would be some gas leakage du e to diffusion mechanism even no microcracks are generated. But the amount would be negligibly sm all compared to the am ount leaked through the cracks and delamination. Hence, in the present work, we focus on the gas leakage due to cracks and delaminations in the composite laminate. The goal of the present work is to develop an analytical numerical tool enabling us to investigate the mechanism of gas permeability in composite laminates, to quantify the gas permeability under given mechanical and thermal lo ads and to provide information for future design of laminates. PAGE 15 15 Literature Review There are a lot of experimental and analytic al investigations on microcracking formation available in the l iterature, e.g., see (Nairn, 2000) and (Berthelot, 2003) and re ferences therein. In the current research, we focus on the analytic al modeling of microcr acking. In conventional stress analysis it is assumed that microcracks form when stress in the 90o-plies reaches the strength of matrix material. This method, howev er, did a poor job to pr edict the microcracking process compared to the experimental observation (Nairn, 2000). Nairn (2000) presented methods based on fracture mechanics and ener gy balance to analyze microcracking, which assumed that new microcracking formed if the total energy released due to the formation of the microcrack reached the microcracking fracture toughness Gmc. These energy-release-rate-based methods did a good job of predicting the microc racking in cross-ply laminates under uniaxial loading. The energy release rate method was then extend by Bapanapalli, Sankar and Primas (2006) to predict crack density evolution in cross-ply laminates under bi-axial and thermal loading using three-dimensional finite elements. This method will be used in the current research to predict the crack density evolution in each ply of the composite laminate. After the failure of the liqui d hydrogen tank of X-33, many e xperimental and analytical works have been done to investigate the gas permeability in composite laminates. Stokes (2002) performed permeability testing on IM7/BMI lamina ted composites under bi-axial strains based on ASTM standard (ASTM D1434-82, 1992). The permeability was measured under constant strain at room temperature and it has been found that the permeability was a time dependent property. Nettles (2001) tested permeability of co mposite laminates afte r impact testing and found that the flow rate usually had a non linear dependence on th e applied pressure. McManus, Faust, and Uebelhart (2001) inve stigated the influence of load ing conditions, crack density and ply orientation on the permeability of graphite/epoxy laminates. Choi (2005) performed PAGE 16 16 permeability tests on various composite material systems subjected to certain cryogenic cycles, and found that for the laminates with same thickness the one with grouped lay-up had higher permeability than the one with dispersed layup. Among all the specimens, textile composite yielded the lowest permeability. Kumazawa, Aoki, and Susuki ( 2003) investigated Helium gas leakage through composite laminates, and compared the experimental results with the analytical models. In their later work (Kumazawa, Susuki, & Aoki, 2006), they measured the gas leakage rate of cross-ply laminates under biaxial load ing and found that besides the amount of the damage the gas permeability was also dependent on both load level a nd biaxial load ratio. Grenoble and Gates (2006) performed permeability tests on composite mate rial IM-7/977-2 after specimen being mechanically cycled at room te mperature and found that the leak rate depended on applied mechanical strain, crack density and test temperature. Bechel, Camping, and Kim (2005) measured the evolution of microcrack density in composites du e to cryogenic thermal cycling. Bechel and Arnold (2006) investigated th e permeability in multi-directional laminates after certain thermal cycling fr om room temperature to cryogeni c temperature and found that the permeability increased significantly as the number of cycling increased. Besides the experimental measurement of the permeability in composite laminates, efforts were also devoted to develop analytical and numerical models to understand the mechanism of permeation. As shown in Figure 1-3, the matrix microcracking and the interlaminar delaminations form the leakage pathways for the cryogenic gas. The intersection areas (the shaded areas shown in Figure 1-3) formed by the crack opening in adjacent plies are the neck areas in the leakage pathway. If we treat the damaged composite laminate as a porous medium, then Darcys law for porous material can be used in modeling gas permeability (Aoki, Kumazawa, & Susuki, 2002, Kumazawa, et al. 2003, Roy & Benjamin, 2004, and Roy & Nair, PAGE 17 17 2006). In Darcys law, the volume fl ow rate of fluid is proportiona l to the pressure gradient, and the proportionality constant is gi ven by the ratio of permeability c onstant to the viscosity of the fluid. The gas permeability constant is determined by the micro-structure of the porous material. In our case, it is determined by crack densities in each ply, ply thicknesse s, intersection area and a material constant that can be determined from experiments. For a given laminate, the ply thicknesses and the material constant are fixed. The crack densities can be predicted using the energy release rate method mentioned before, or di rectly measured from experiments. However, the crack opening of the microcracks is hard to measure in experiments. Hence, analytical models are required to pred ict the intersection area. Kumazawa et al. (2003) only considered microcracking in their model and took mean crack opening displacement computed from two-dimensional shear-lag analysis to get the intersection area. Theoretically, however, microcracking along cannot form the intersection ar ea; delaminations have to be introduced in the model. Roy and Benjamin (2004) developed a mode l with delamination and computed crack opening displacement based on firstorder shear deformable laminate theory. They compared the analytical results for crack opening displacem ent with finite element (FE) results and investigated the permeability under mechani cal loading and thermal loading. Later on, Roy and Nair (2006) extended this model to include the st itch cracks developed in the oblique layers of the laminate by introducing effec tive spring stiffness. All these models are 2-D models. In the current work, we have found that the intersec tion areas were largely dependent on the local geometries, and a full three-dimensional analysis is required to obtai n accurate result (Xu & Sankar, 2007). Some other researchers predicted the gas pe rmeability in composite laminate based on fluid mechanics principles (Noh, Whitcomb, Peddiraju, & Lagoudas, 2004; Peddiraju, Lagoudas, PAGE 18 18 Noh, & Whitcomb, 2004; Peddiraju, Grenoble, Fried, Gates, & Lagoudas, 2005; and Peddiraju, Noh, Whitcomb, & Lagoudas, 2007). In their models, the opening due to transverse microcracks and delaminations and the inte rsection area are first investig ated. Then a damage network through laminates with assumed crack opening a nd delamination was set up and the effective conductance through this network was predicte d using the computational fluid dynamics program FLUENT. In these studies, the effective conduc tance was related to the ply thickness and the crack density in each ply. Also, the cond uctance of gas through a representative volume element (RVE) of intersection and microcracks were computed and a simplified model was set up, which agreed well with the previous simulati on results. In these models, however, since the damage network was assumed with fixed opening and intersection area, each model can only represent one single loading case because the op ening and the intersection area are dependent on the applied loads. Hence these models have diffi culties to relate conducta nce to applied loads. Most of the analytical models focused on th e cross-ply laminates, in which transverse matrix cracks developed in each ply and extende d completely in both thickness direction and width direction. In general, it is good to include angle-plies in composite laminate design to gain better performance. Different from extensive tran sverse matrix cracks in cross-ply laminates, short cracks (or called stitch cracks) emanated in angle plies in multi-directional laminates. This phenomenon was first identified by Jamison, Schulte, Reifsnider, & Stinchcomb (1984) in the multi-directional laminates under fatigue loading. Lately, the stitch cracks formed in laminates with stacking sequence [+ n/n/902 n]s were systematically investigated by Lavoie and Adolfsson (2001). It has been found that th e stitch cracks not only emanated under fatigue loading but also emanated due to monotonic loading and thermal residual stresses. Stit ch cracks were also observed in [0/+45/-45/90]s laminates by Bechel et al. (2005) afte r certain thermal cycling PAGE 19 19 between room temperature and cryogenic te mperature. Yokozeki and his colleagues (Yokozeki, Aoki, Ogasawara, & Ishikawa (2005); Yokozeki & Ishikawa (2005); and Yokozeki, Aoki, & Ishikawa (2005)) performed experiments on laminates of type [0/ n/90]s under uniaxial loading and they found that whether stitch cracks or developed cracks will emanate in the o-ply depended on the value of and the thickness of the o-ply. Some of the aforementioned analytical models included stitch crack s in the modeling. As mentioned before, Roy and Nair (2006) introduced an effec tive spring stiffness to incorporate the stitch crack into their model and still used the first-order shear deformable laminate theory to compute the crack opening displacement of stitch crack. In the damage network, Peddiraju et al. (2005) included the stitch cracks with assumed crack length and crack opening. Note that all these models either eliminated the difference betw een stitch crack and tran sverse crack using an effective spring stiffness or assumed damage st ate of stitch cracks. They were not able to represent the relationship between stitch crack and lo ads. In reality, the crack densities, the crack length and the opening displacement of stitch cracks are functions of the applied loads. Objectives In the current work, the objective is to develop a method based on Darcys law and three-dimensional finite element analysis (FEA ) to predict the gas permeability in composite laminates. By using this method, the permeability in both cross-ply laminates and multi-directional laminates will be predicted as a function of applied load. Also, the effects of different factors, such as ply thickness, delamination length, delami nation shape, lay-ups, etc., on the permeability will be investig ated and some useful suggestions will be provided for the future design of composite cryogenic tanks. PAGE 20 20 Figure 1-1. Lockheed Martin X-33 (http://www.makelengin eering.com/dir/Compa ny/History/History.htm ) Figure 1-2. Microcracks in the 90o-ply and interfacial delaminations between 0oand 90oplies due to thermal/mechanical loading. PAGE 21 21 Figure 1-3. Gas permeation pathway provided by microcracks and delaminations in cross-ply laminates. Intersection Area Crack Opening Displacement PAGE 22 22 CHAPTER 2 MODELS AND METHODS There are several theories that might be us eful in developing the permeation model for composite laminates: Darcys Law for porous ma terials, Poiseuilles Law for capillary and Darcy-Weisbach Equation for pipe pressure loss These equations are listed in Appendix A and they can be derived from each other. In the present work, the permeation model in composite laminate is mainly derived from Darcys law fo r porous material and veri fied by other laws. The permeability is represented in terms of crack density, intersection area and a material constant to be determined by experiments. The derivation of the permeation model is based on the cross-ply laminate, where only 0o-ply and 90o-ply present, but it can be extended to multi-directional-ply laminate as well. Permeation Model From experimental observations, matrix transverse micro-cracking in the 90o-ply is the first form of damage in cross-ply composite lamina tes. Due to stress singularity at the crack tip when transverse crack reaches the interface between 90 o-ply and 0o-ply, interlaminar delamination may occur at the interface. Theo retically, microcracks alone cannot form a contiguous path for gas leakage since no intersection area is form ed. Hence, the delamination has to be introduced in the numerical model as shown in Figure 1-3. The shaded intersection areas formed by the crack opening displacements (CODs) of adjacent plies connect the micro-cracks, and allow gas to leak out. If the composite lamina te with micro-cracks and delaminates is treated as a porous material, Darcys Law fo r porous materials can be applied. Darcys law for viscous flow of gases through porous media is given as P q (2-1) PAGE 23 23 where q is volume flow rate per unite area, is the permeability tensor, is the viscosity and P is pressure gradient vector. If we only consider the gas leakage in the thickness direction of the laminate, following one-dimensional vers ion of Darcys law can be used: z P k q (2-2) where k is the overall permeability constant in the thickness direction and z is the thickness coordinate. The perm eability constant k is dependent on the microstructure of the laminate. Integrating both sides of Equation 2-2 over laminate thickness gives 0P P k qhN (2-3) or k qh P PN 0 (2-4) where h is the laminate thickness, PN and P0 are pressure at the top and bottom respectively (Figure 2-1). Also, if the integration is done layer wise, we have N i k qh P Pi i i i,..., 1 ,1 (2-5) where hi and ki are thickness and permeability for each lamina. Summing Equation 2-5 from i = 1 to i = N yields N i i i Nk h q P P1 0) ( (2-6) Comparing Equation 2-6 with Equation 2-4, it can be seen that N i N i i i i ik h h k k h k h1 1or (2-7) For each lamina, most of the resistance come s from the intersection area. According to Poiseuilles Law, the volume flow ra te in a rectangular tube (Figure 2-2) for the laminar flow is given by: PAGE 24 24 P L d Q 124 (2-8) Dividing both sides by the cros s-section area of the tube, d 2, gives: L P d q 122 (2-9) Comparing Equation 2-9 with Equation 2-2, it can be seen that the permeability constant k for the rectangular tube is propor tional to the crosssection area. Based on this fact, we can assume that the permeability for each lamina is proportional to the total intersection area. This assumption is also confirmed by Peddiraju, et al. (2007) using FLUENT simulation. Noting that first and last lamina only has one interface with other laminas while internal laminas have two interfaces, so the individual permeability is assumed to take following form: ) 1 ( 1 ) 1 ( 1 ) 1 ( 1 ) 2 1 ( 2 1 11 ,..., 2 2 /N N N N N i i i i i i i i iA C k N i A A C k A C k (2-10) where i is the crack density in ith ply, ( i i +1) is the individual intersection area between ith and i +1th ply, and C is a material property to be determined by experiments. If we substitute ki back into Equation 2-7, the overall permeability of the laminate can be written as: ) 1 ( 1 1 2 ) 1 ( 1 ) 1 ( 1 ) 2 1 ( 2 1 12N N N N N N i i i i i i i i i iA h A A h A h Ch k (2-11) It is obvious from Equation 2-11 that the permeability is determined by the crack densities i and the intersection areas A(i,j) once the laminate lay-up is given and the material constant C is determined by experiments. In another word, the overall permeability is determined by the total intersection area in each ply because that the crack densities give the number of intersection area in each ply. Hence, it is easy to extend Equation 2-11 to multidirectional laminates if the total intersection area in each ply is obtained. PAGE 25 25 Preliminary Two-Dimensional FE Modeling Before using a three-dimensional finite elem ent model, two-dimensional FE modeling is presented first to obtain some preliminary results. A three layer cross-ply laminate is considered in this modeling. Figure 2-3 shows the unit cells of two differe nt cases of damage scenarios that could occur in cross-ply composite laminates. In the first case, [90/0/90], the delamination emanates from the surface ply of the laminate wh ereas in the second case, [0/90/0], it emanates from the middle ply. Due to symmetry, only one-f ourth of the unit cell is considered in the simulation. Eight-node isoparametr ic elements are used in the simulation and the orthotropic material properties used in the simulation are shown in Table 2-1. The simulation has been done using the commercial finite element package ABAQUS 6.5. From Figure 2-3, it can be seen that the unit cell can be described by the crack spacing (reciprocal of the crack density), the ply thickne ss and the delamination length. The ply thickness h is taken as 0.33mm in the current model. Firs t, the variation of stra in energy release rate (SERR) at the delamination tip with respect to the delamination length a is investigated. A fixed crack density of 1 cm-1 is taken. Even though the individual Mode I and Mode II strain-energy release rates do not exist due to the oscillatory ch aracteristic of stresses and displacements near the crack tip for bi-mat erial interface cracks (Sun & Jih, 1987), the total strain-energy release rate is well defined (Hutchinson & Suo, 1992). Under a line ar elastic assump tion, the path independent J-integral is identical to the strain energy release rate. Since J-integral is insensitive to the FE mesh, uniform mesh is accurate enough for this simulation. Figure 2-4 shows the J-integrals with re spect to the delamination length a The delamination length a is normalized by the ply thickness h in Figure 2-4. The unit cell is under constant displacement loading. It can be seen that the J-integrals reach maximum values when the delamination length is about the ply thickness, and then decrease slowly. This implies that PAGE 26 26 the delamination propagates steadily when the delamination length exceeds the ply thickness under displacement controled loadin g. When the delamination length is too small, less than the ply thickness, the boundary effects are obvious. From Figure 2-4, it seems to be a proper choice to take ply thickness as the delamination length since the delamination propagates steadily after its length exceeds ply thickness. Hence, in three-dimensional FE analysis, a length compar able to the ply thickness is taken to be the delamination length and the sensitivity of the intersection area to the delamination length is investigated. Next, the COD of the microcrack is com puted with a fixed delamination length and varying crack density. Even though both the crack density and intersection area are functions of applied loads, fortunately, the following results show that thes e two factors can be decoupled, since crack density does not ha ve much influence on the COD. Figure 2-5 shows the computed COD with respect to crack density. The CODs are expressed for unit force by normalizing w ith respect to the force resultant Nx, which can be computed directly from FEA by volume average method (average stress e quals to the summation of production of stress and element volume divi ded by the total element volume, and force resultant is given by the production of average stress and laminate th ickness). It can be seen from these two figures that the CODs are almost cons tants for a given force resultant and delamination length. From the above discussion, we can argue th at the crack density does not influence the crack opening displacement much. Hence the perm eability prediction can be decoupled into two parts, first find crack densities for given for ce resultants and then compute intersection area under these loads. Then, the permeability for a gi ven set of force resultants can be computed PAGE 27 27 using Equation 2-11. The crack densities in each ply can be either predicted using analytical methods (Nairn, 2000 and Bapanapalli, et al., 2006) or meas ured directly from experiments (Nairn, et al., 1993 and Nairn, 2000). However, most of the av ailable experimental data are for laminates under uniaxial loading. In reality, the composite lami nate in a cryogenic tank is under biaxial loading. Hence, the analytical method developed by Bapanapalli, et al. (2006) will be adopted in the current work to predict crack densities in cross-ply laminates under biaxial loading. The intersection area will be computed using three-dimensional finite element method (FEM) and a parametric investigation will be performed. Table 2-1. Orthotropic material properties for the epoxy/graph ite composite material in the present study. (Bapanapalli, et al. 2006) Properties Value E1 169 GPa E2, E3 8.62 GPa G12, G13 5.0 GPa G23 1.22 12, 0.355 23 0.41 -3.384e-9/oC 28.998e-6/oC PAGE 28 28 Figure 2-1. Pressure distribution alo ng thickness direction of laminate. Figure 2-2. Laminar flow in a rectangle tube. P 0 P 1 P 2 P N PAGE 29 29 Figure 2-3. Two-dimensional delamination model. A) Delamination from surface plies [90/0/90]. B) Delamination from middle ply [0/90/0]. J-integral versus a (90/0/90) Delamination length ( a / h )J-integral, J/m2 J-integral versus a (0/90/0) Delamination length ( a / h )J-integral, J/m2 A B Figure 2-4. J-integral vari ation with respect to norm alized delamination length a/h A) [90/0/90] laminate. B) [0/90/0] laminate. (a) (b) a a A B PAGE 30 30 90/0/90 Crack Density, cm-1COD/Nx, m/(N m) 0/90/0 Crack Density, cm-1COD/Nx, m/(N*m) A B Figure 2-5. COD variation with respect to crack density. A) [ 90/0/90] laminate. B) [0/90/0] laminate. PAGE 31 31 CHAPTER 3 CRACK DENSITY PREDICTION In this chapter, the method developed by Bapanapalli, et al. (2006) is briefly summarized; predicted crack densities under uniax ial loading is compared with e xperimental results; and then crack densities under biaxial loadi ng are predicted for various cases. Methodology Cross-ply composite laminate [0/90/0] is considered as an example in this section to understand the methodology. From experimental observations (Nairn, 1993), two different microcracking scenarios could develop in the la minate when the laminate was subjected to uniaxial loading as shown in Figure 3-1. If the laminate is loaded in the 0o-direction, microcracks emanated in the 90o-ply (middle ply) and next new microc rack usually formed in between two existing cracks (Figure 3-1A). On the other hand, if th e laminate is loaded in the 90o-direction, microcracks emanated in the 0o-plies (surface plies). The asy mmetric crack pattern shown in Figure 3-1B has been observed by Nairn (2000). Hence, four new cracks formed to keep the same pattern. A finite fracture mechanics approach was used by Nairn (2000) to pred ict new microcrack formation. In this approach, the strain energy release rate due to the formation of new microcracks is expressed as a function of applie d load, mechanical properties and crack density. Then this strain energy release rate is equated to microcracking fracture toughness, Gmc, to solve for the required load to form new microcracks. It has been shown that the microcracking fracture toughness Gmc is a material property and it can be de termined through experimental techniques (Nairn, 2000). The predicted results using this ap proach matched the experimental results quite well. However, this approach can only predict th e crack density of compos ite laminate subjected PAGE 32 32 to uniaxial loading. Bapanapalli, et al. (2006) extended this approach to predict crack densities in both 0o-ply and 90o-ply at same time when the laminate is subjected to biaxial loading. In the current work, the extended finite fracture mechanics approach (Bapanapalli, et al., 2006) is used to predict crack densities in ea ch ply of composite laminate. Based on the two scenarios shown in Figure 3-1, a three-dimensional unit cel l is taken as shown in Figure 3-2. In this model, crack surfaces are shown in shaded area. The dimension of the unit cell is determined by the crack densities of surface and middle plies. The strain energy of the unit cell fo r load control condition is given by N A N GT y x 1 1 2 1 (3-1) where 1 22 21 12 11 A A A A A is the inverse of laminate extensional stiffness matrix (or the compliance matrix), y xN N N is force resultant matrix, and x, y are crack densities in xand y-direction. The formation of new microcracks will cause the degradation of the stiffness matrix of the laminate. Hence, the released strain en ergy due to the formation of new microcracks is given by N A N GT y x m 1 1 2 1 (3-2) where 1 2A A A is the difference of the compliance matr ix before and after the formation of new microcracks. As the wall of cryogenic tank, the composite laminate is unde r biaxial loading, and the hoop stress and the longitudinal stre ss are proportional to each othe r. Hence, it is reasonable to assume that the force resultants in the laminate in the form: y xN N (3-3) PAGE 33 33 where is a proportional factor, whose value equals to 2 or 0.5 for cylindrical body and 1 for the spherical cap of a pressure vessel. Substituting Equation 3-3 into Equation 3-2, the released strain energy can be rewritten as 2 22 12 11 22 1 1 2 1y y x mN A A A G (3-4) where 11A 12A and 22A are the entries of matrix A According to the finite fracture mechanics, the strain energy releas ed due to new crack formation should equal to the energy required to create new crack surfaces There are three cases could happen for the microcrack formation: Case I: new microcracks formed onl y in surface plies. In this case, x remains the same but y is tripled. Case II: new microcracks formed onl y in middle ply. In this case, x is doubled and y remains the same. Case III: new microcracks formed in both middle and surface plies. In this case, x is doubled and y is tripled. The required surface energies for these three cases are given by mc x mc y III mc y II mc x IG h G h U G h U G h U0 90 90 01 4 1 1 1 4 (3-5) where h0 and h90 are the thickness of 0o-ply and 90o-ply respectively. Equating these surface energies to the released strain energy in Equation 3-4, one can solve for the required force resultants for the three cases as follows 22 12 11 2 0 ) (2 8A A A G h Nmc y I y (3-6) PAGE 34 34 22 12 11 2 90 ) (2 2A A A G h Nmc x II y (3-7) 22 12 11 2 0 90 ) (2 8 2A A A h h G Ny x mc III y (3-8) In reality only one case could happen, so the case that gives the smallest value of Ny will happen. In another word, for a given unit cell (x and y fixed), if Ny(II) is the least values among three computed force resultants, then next microc rack will form in surface ply. Note that the values of A are different for these three cases. This method can also be extended to displa cement control condition. For details, please refer to Bapanapalli, et al., (2006). Predicted Results Uniaxial Loading First, the crack density in a composite la minate is predicted under uniaxial loading condition. To compare with the availabl e experimental data, laminate [0/90n/0] (n=1, 2, 4) are investigated. The material properties ar e taken to be the same as in Table 2-1. In Equation 3-6 to Equation 3-8, for each combination of (x, y), the stiffness matrix [A] is different. Matrix [A] is computed using FE an alysis, but if we run the FE analysis for every unit cell (given x and y), there will be too much work si nce 3-D FE analysis is very time consuming. Hence, in the current work, the [A] matrices for certain (x, y) combinations are computed and then stiffness matrix is fitted as a complete cubic function of x and y based on the computed values using th e least square approximation. For a given unit cell, the [A] matrix is computed based on constitutive relation: y x y xA A A A N N 22 21 12 11 (3-9) PAGE 35 35 Displacement boundary conditions are applied to th e unit cell to make the strain status be x = 1 and y = 0. The force resultants Nx and Ny are computed using volume average method. Then the first column in [A] matrix equals to the computed force result ants. Similarly, applying the boundary condition x = 0 and y = 1 will give the second column of [A] matrix. For each pair of (x, y), two FE analysis are needed to determin e [A] matrix. In our case, seven values of x (ranging from 0 to 9.5 cracks/cm) and eight values of y (ranging from 0 to 12.6 cracks/cm) are taken. So totally 112 FE analysis are needed fo r each laminate to dete rmine response surface of [A] matrix. It can be seen that this is a quite time consuming procedure. Figure 3-3 shows the response surfaces for the entries in [A] matrix of laminate [0/90/0]. The dot marks are data points. Other laminates also have similar response surfaces. A MATLAB program is developed to pred ict the progressive crack densities x and y of laminate subjected to biaxial loading as a function of Nx or Ny. In this program, initial values of crack densities and the ratio between Nx and Ny are the inputs and crack densities in both 0o-ply and 90o-ply are outputs. The crack densities are plot ted with respect to the applied loads. For uniaxial loading, the ratio is taken as 0. The MATLAB code s and the derivation of the least square approximation of [A] matrix are enclosed in Appendix B. From experimental observations (Nairn, 2000) in laminates [0/90n/0], it has been found that the initiation stress for microcracks in the middle ply is less for the laminate with a thicker middle ply, and the initiation stress for microcracks in the surface ply is greater for the laminate with a thicker middle ply. However, at high stress levels, a larger crack density was observed in the middle ply of the laminate with the thinne r middle ply. The predic ted crack densities of laminates [0/90n/0] (n = 1, 2, 4) under uniaxial loading are shown in Figure 3-4 and Figure 3-5. The microcracking fracture toughness Gmc is taken as 240 J/m (Nairn, 2000). In Figure 3-4, the PAGE 36 36 stress to initiate first microcrack increases as the middle ply thickness decreases. However, once the initiation stress is reached, the crack density increases faster in the laminate with a thinner middle ply. For a [0/90/0] laminate, a negati ve slope is obtained, which means that the microcracks accumulate simultaneously in the laminate. As shown in Figure 3-5, as the middle ply thickness of the laminate increases, the initiation stress increases as well. Hence, our prediction and experimental observation matches quite well qualitativ ely, which gives us confidence to use this approach to predict the crack densities of laminates under biaxial loading. Biaxial Loading In reality, the laminate used in propellant tank is under bi-axial lo ading. The proportional constant between force resultants Nx and Ny is taken as 0.5, 1.0 and 2.0 to represent different possible stress states in hydrogen tank. Figure 3-6 to Figure 3-14 show the variations of crack densities in different laminates under different biaxial loading. Fo r each laminate, crack densities x and y are plotted in the same figure with respect to the average stress in the y-direction (Ny divided by the laminate thickness). Figures 3-6 through 3-8 represent the cases when = 0.5, Figures 3-9 through 3-11 represent = 1.0, and Figure 3-12 through 3-14 represent = 2.0. From these figures, it can be seen that the ini tiation stresses for first microcrack in both xand y-direction are lower than their corresponding valu es under uniaxial loading. With the same proportionality constant the initiation stress for x decrease dramatically as the middle ply thickness increase. On the contrary, the initiation stress for y increases as the middle ply becomes thicker, but the increase step is not as great as that under uniaxial loading. Loads in the x-direction affect the crack densities in the y-direction. For the same laminate, say [0/90/0], when increases (x increases), the initiation stress for y decreases. No microc racks are generated in PAGE 37 37 [0/90/0] for = 0.5 and 1.0 when y already meet the preset limit (saturation status). These predicted crack densities will be used to pred ict the progressive gas permeability later on. A B Figure 3-1. Two microcracking s cenarios in composite laminate [0/90/0] under uniaxial loading. A) Unit cell for microcracks in middle ply and new crack formation. B) Unit cell for microcracks in surface ply and new cracks formation. (Courtesy: Bapanapalli, et al., 2006) Figure 3-2. Three-dimensional unit cell with microcracks in both middle and surface plies. 0o 90o0o PAGE 38 38 A B C Figure 3-3. Response surface for [A] ma trix of laminate [0/90/0]. A) A11. B) A12. C) A22. PAGE 39 39 (MPa)x (cm-1) Figure 3-4. Crack density evolution in the middle ply under loading in x-direction. y (MPa)y (cm-1) Figure 3-5. Crack density evolution in the surface ply under loading in y-direction. PAGE 40 40 = 0.5, [0/90/0]0 2 4 6 8 10 12 14 02004006008001000y (MPa) (cm-1) x y Figure 3-6. Crack density pred iction for laminate [0/90/0]. Nx/Ny = 0.5 = 0.5, [0/902/0]0 2 4 6 8 10 12 14 02004006008001000y (MPa) (cm-1) x y Figure 3-7. Crack density prediction for laminate [0/902/0]. Nx/Ny = 0.5 PAGE 41 41 = 0.5, [0/904/0]0 2 4 6 8 10 12 14 02004006008001000 (MPa) (cm-1) x y Figure 3-8. Crack density prediction for laminate [0/904/0]. Nx/Ny = 0.5 = 1.0, [0/90/0]0 2 4 6 8 10 12 14 0100200300400500600700 (MPa) (cm-1) x y Figure 3-9. Crack density pred iction for laminate [0/90/0]. Nx/Ny = 1.0 PAGE 42 42 = 1.0, [0/902/0] (MPa) (cm-1) x y Figure 3-10. Crack density prediction for laminate [0/902/0]. Nx/Ny = 1.0 = 1.0, [0/904/0]0 1 2 3 4 5 6 7 8 0100200300400500600 (MPa) (cm-1) x y Figure 3-11. Crack density prediction for laminate [0/904/0]. Nx/Ny = 1.0 PAGE 43 43 = 2.0, [0/90/0]0 1 2 3 4 5 6 7 8 9 10 0100200300400500 (MPa) (cm-1) x y Figure 3-12. Crack density pred iction for laminate [0/90/0]. Nx/Ny = 2.0 = 2.0, [0/902/0]0 2 4 6 8 10 050100150200250300350400 (MPa) (cm-1) x y Figure 3-13. Crack density prediction for laminate [0/902/0]. Nx/Ny = 2.0 PAGE 44 44 = 2.0, [0/904/0]0 2 4 6 8 10 050100150200250300350400 (MPa) (cm-1) x y Figure 3-14. Crack density prediction for laminate [0/904/0]. Nx/Ny = 2.0 PAGE 45 45 CHAPTER 4 INVESTIGATION ON INTERSECTION AREA After predicting crack density in the previous chapter, we will focus on the intersection area and its variation due to diffe rent factors in this chapter. Model Description Three layer cross-ply composite laminate [0/90/0] is considered in the current model and the following assumptions have been made for simplicity: the crack densities and the position of microcracks in the two surface plies are identi cal; the microcracks in both surface plies and middle ply are evenly distributed. Once crack dens ities are given, the dimensions of the unit cell are fixed as shown in Figure 4-1. The upper figure in Figure 4-1A represents th e top view of the laminate with transverse matrix microcracks. Th e horizontal and vertical solid lines represent the microcracks in 0oand 90oplies, respectively. Due to both in-plane symmetry and symmetry in thickness direction, only one-eight h of the unit cell is analyzed as shown in the bottom figure in Figure 4-1A. Two free surfaces are the microcra ck faces. Symmetry boundary conditions are applied to the surfaces opposite to the free surfaces and the bottom surface. Constant displacement boundary conditions are applied as shown by the arrows. Delamination is introduced in this model since no contiguous gas-leak-pathway will be formed without delamination. The delaminations emanated from both surface ply and middle ply form a rectangular delamination front. In the curr ent model, the crack densities in the 0oand 90o-plies are taken as equal, x = y = 1cm-1. The material properties are ke pt identical to those of the two-dimensional model as shown in Table 2-1. Ply thickness is still taken as 0.33 mm. Twenty-node 3-D isoparametric solid elements a nd uniform mesh are used in the simulation. About 28,000 elements were used in the model. Figure 4-1B shows a typical deformation of the model under bi-axial loading. Note that the deformation is enlarged by a factor of 10. PAGE 46 46 Kumazawa et al. (2003) and Roy and Benjamin (2004) proposed to superpose the two-dimensional CODs obtained from analytical methods to get the intersection area. In fact, if we monitor the COD along the path shown by the arrow in Figure 4-2A, the COD first remains as a constant at the straight edge, and then in creases remarkably close to the rectangular crack front. The constant value along the straight e dge is identical to th e result predicted by two-dimensional analytical model. The COD near the rectangular delamination front, however, is almost 50% larger than the two-dime nsional result as shown in Figure 4-2B. Therefore, the local effects cannot be neglected and COD near the rectangular front computed from three-dimensional analysis should be us ed to compute the intersection area. Intersection Area Computation Due to the linear relation between COD and load, the intersection area under any loading case (Nx, Ny, T) can be obtained once the CODs for the following three basic cases are computed: ; 0 0 %, 1 Ty x ; 0 %, 1 0 Ty x re). temperatu cryogenic to re temperatu room (From 223 0 0 C T o y x For any combination of (Nx, Ny, T), the corresponding (x, y, T) can be obtained from the basic formula for composite laminate shown in Equation 4-1, in which the matrices [A] and [] can be determined from FE analysis of thos e three basic cases using the approach mentioned in Chapter 3. y x y x y xT A A A A N N 22 21 12 11 (4-1) The corresponding CODs for (Nx, Ny, T) are then obtained by superposing the results of three basic cases with different weighting factor s. The intersection area is computed by assuming PAGE 47 47 that the cross-section formed by the C ODs at corner is a rectangle. Figure 4-3 shows the variation of the intersection area with respect to the mechanical and th ermal loads. In Figure 4-3A, different symbols repres ent different ratios between Ny and Nx, the solid lines represent the cases with thermal load and the dotted lines repr esent the cases without thermal load; in Figure 4-3B, each surface represents the intersec tion area variation with respect to Ny and Nx under certain thermal load. Actually, from those three basic numerical test s, the CODs can be expressed explicitly in terms of Nx, Ny, and T. T e N e N e uy x 8 11 11082 1 446 0 592 0 (4-2) T e N e N e vy x 8 11 11036 1 156 1 204 0 (4-3) The corresponding intersection area A is given by T N e T N e N N e T e N e N e v u Ay x y x y x 18 19 22 2 15 2 22 2 22171 0 834 0 775 0 112 0 516 0 121 0 4 (4-4) From Equation 4-4, it is obvious that the intersection area shows pa rabolic variations with respect to Nx, Ny, and T, which can be seen clearly in Figure 4-3. Equations 4-2 through 4-4, however, are only valid for the current model. If any condition of the model, such as delamination length, material properties and ply thickness, is changed, we have to redo the computations of those three basic cases. Effects of Delamination Shape and Delamination Length The effects of delamination shape and the de lamination length on the intersection area are investigated in this section. In the previous section, the dela mination fronts developed from the surface ply and from the middle ply were connect ed by a pair of straight lines forming a rectangular delamination front. Ho wever, if three-dimensional Jintegrals are computed along the crack front as shown in Figure 4-4, it can be seen that the greates t value occurs at the corner, i.e., PAGE 48 48 the crack will first propagate at the corner. Hence, the delamination with a circular front might be a more appropriate scenario. Figure 4-4A shows the top view of quarter of th e unit cell with delaminations, and Figure 4-4B shows the distribution of Jintegrals with different connecti ons. The model is subjected to thermal load only. As we can see the J-integral va lues remain as constant at the straight crack front and increase near the connecting corner. For the case of the rectangular front, the J-integral value computed at the corner, which is extremely high due to the discontinuity of the first order derivative of the geometrical path, may not be meaningful and is eliminated from Figure 4-4B. As we expected, the circular dela mination front reduces the highest value of J-integral at corner, and bigger is the radius, lower is the value. The intersection areas in the models with different connec ting radii are computed. Figure 4-5 shows the effects of the ra dius of the connecting circle. In all three basic cases, the cross-sectional area increases quadratic ally with respect to the radius. Moreover, the effect of delamination length is also investigated, as shown in Figure 4-6. In this case, the rectangular front is still used and the delamination length changed from half of the ply thickness to twice the ply thickness. It can be seen that the intersecti on area is more sensitive to the delamination length than to the circle radi us. It does not show the tendency to be a plateau as mentioned by Roy and Benjamin (2004). Single Delamination Model Another phenomenon shown in Figure 4-4B is that the steady state values on different sides are quite different even though both directi ons are under the same loading condition. Hence, delaminating from the surface ply (high J-integral ) is much easier than delaminating from the middle ply (low J-integral). From a different view of the typical deformation of the unit cell as shown in Figure 4-7, it can be seen that delami nation from the middle ply (Figure 4-7A) tries to PAGE 49 49 close while the delamination from the surface ply (Figure 4-7B) tries to open up. Therefore, it is much more possible that the delamination only emanates from surface plies. Even if only one-side of the delamination appears, the path fo r gas leakage can still be formed as shown in Figure 4-8. We call this model as the Si ngle Delamination Model (SDM), and, correspondingly, the previous model is ca lled the Double Delamination Model (DDM). In this case, the intersection area is smalle r than the one in the previous model. Figure 4-9 shows the comparison between the intersection areas computed using SDM and DDM under the same loading condition ( Ny/Nx = 0.5 and T = -223oC ). The intersection area A computed using SDM also shows a parabolic variation, bu t the magnitude is only about half of that computed using DDM. Effect of Lay-Ups Until now, only three-ply laminates are investig ated. In reality, there are various lay-up possibilities. Choi (2005) performed permeability tests for different laminates and found that for the same thickness the specimen with layers gro uped together has greater permeability than the one with layers dispersed, i.e., laminate [0/90/02/90/0] performs better than laminate [02/902/02]. Similar conclusion was also obtained by Kumazawa, Hayashi, Susuki & Utsunomiya (2006) in the permeability test for laminates [0/90/0/90]s and [0/0/90/90]s under bi-axial loading. In this section, the current modeling technique will be used to demonstrate this phenomenon. In this investigation, the exact same lay-up and material properties ar e used as in Chois tests (Choi, 2005). For laminate [02/902/02], the unit cell is kept same as the one described in previous sections. For laminate [0/90/02/90/0], in order to describe the unit cell, the micro-cracks in the 0o-plies and in the 90o-plies are assumed to appear in same position and the crack densities in all 0o-plies or in all 90o-plies are equal, as shown in Figure 4-10A. Due to the symmetry of the PAGE 50 50 model, only one-eighth of the un it cell is considered. Figure 4-10B shows the typical deformation of the unit cell. The crack densities in laminate [02/902/02] and in laminate [0/90/02/90/0] are kept the same. The permeabilities in both lamina tes are computed using Equation 2-11. Since the material constant C is to be determined, the permeability is normalized by the constant C. In Table 4-1 the normalized permeability k/C for both specimens are tabulated and compared with the experimental results by Choi (2005). From this table, it can be seen clearly th at the specimen with grouped layers has a much higher permeabil ity than the specimen w ith dispersed layers. Note that the ratios between the permeability of the two specimens are quite different in the models and the experiments. This is due to the fact that the modeling procedure does not precisely reflect the testing pro cedure. In Chois tests, the sp ecimens were tested at room temperature after certain cryogeni c cycles (cooling specimen down to cryogenic temperature and then warming it up to room temperature is define d as one cryogenic cycle) without mechanical loading. While in our modeling, crack density and delamination length is assumed and bi-axial loading is applied. Hence, the two sets of re sults can only be compared qualitatively for now. Temperature Dependent Properties In this section, the effect of temperature dependent properties is investigated. In previous sections all material properties are assumed to be temperature inde pendent. For graphite/epoxy composite laminates, the matrix properties la rgely depend on temperature, while the fiber properties do not. Hence, the composite properties that depend mainly on matrix properties, such as transverse modulus E2 and shear modulus G12, vary remarkably with respect to temperature. Schulz (2005) and Speriatu (2005) have performed some ex periments to measure the temperature dependent properties for graphite/epoxy co mposite IM7/977-2. Transverse modulus E2, shear modulus G12, and coefficient of thermal expansion were measured at different temperature and PAGE 51 51 fitted with polynomial functions. In the following simulation, these experimental results are used. Longitudinal modulus E1 and Poissons ratios are assumed to be temperature independent, and shear modulus G23 is assumed to follow the relation 23 2 231 2 E G (4-5) In the simulation, constant strain boundary conditions are applie d. Strains in both x-direction and y-direction are 1%. Thermal load is appl ied from room temperature to cryogenic temperature. Figure 4-11 shows the computed intersection ar ea. The line with square symbols is the result with temperature independent material properties and the line with diamond symbols is the result with temperature dependent properties. From this figure, it can be seen that at cryogenic temperature the result with temperature dependent prope rties is almost 15% less than that with temperature independe nt properties. Hence, if temperature independent material properties were used, the simulation will give the upper bound of the permeability. Progressive Permeability Prediction Based on the argument drawn from the two-dimens ional FE analysis in Chapter 2, the two determining factors of gas perm eability can be decoupled. In th is section, progressive gas permeability in laminates [0/90n/0] (n = 1, 2, 4) is predicted in te rms of applied load. In Chapter 3, the crack densities in laminates [0/90n/0] are predicted under different biaxial loadings. In this section, a typical unit cell is taken to compute th e intersection area. To keep the same as the unit cell in crack density prediction, the surface cracks in the unit cell are not symmetric. The crack densities in this unit cell is taken as x = y = 2 cm-1. Delaminations from both surface ply and middle ply are introduced in this model. Due to symmetry, only one-fourth of unit cell is taken for analysis, as shown in Figure 4-12. In this section, only mech anical loads are considered, so only two basic cases are need ed for each unit cell: (x = 1%, y = 0) and (x = 0, y = 1%). Figure PAGE 52 52 4-13 shows the deformation of the unit cell under loading case (x = 1%, y = 0). The deformation is enlarged by a factor of 10. The intersection areas of the two basic cases for laminates [0/90n/0] are tabulated in Table 4-2. From this table, it can be seen that in both cases the intersection area increases as the thickness of the middle ply increases. A MATLA B program is then developed to predict progressive permeability and to plot the permeability versus the applied load. The program is attached in Appendix C. In this program, th e predicted crack densit ies shown in Figures 3-6 through 3-14 are fitted with cubic polynomials and th en the intersection areas are computed from the two basic cases shown in Table 4-2. Finally, the progressive permeability is computed using Equation 2-11. Figures 4-14 through 4-16 show the predicted permeability as a function of average stress y. Each figure represents one loading ratio. In Figure 4-14, Nx/Ny = 0.5, in Figure 4-15, Nx/Ny = 1.0, and in Figure 4-16, Nx/Ny = 2.0. From these figures, it can be seen that the initiation stresses of gas permeability increase as the thickness of middle ply decrea se. After initiation, however, the gas permeability in the laminate with thinner middle ply increases faster than in the one with thicker middle ply. This follows the same tre nd as the evolution of crack density. Even though the intersection area in [0/904/0] is the largest, th e permeability in [0/904/0] increases at the slowest rate comparing to the other laminates. He nce, it can be concluded that the permeability is dependent more on the crack densit y than on the intersection area. For laminate [0/90/0], no permeability is predicted for Nx/Ny = 0.5 and Nx/Ny = 1.0 under available stress level because the middle ply remains undamaged (Figure 3-6 and Figure 3-9). PAGE 53 53 Table 4-1. Comparison between co mputed and tested permeability Laminates Normalized Permeability k / C Test results by Choi (2005) (mol/sec/m/Pa) [02/902/02] 1.2610-8 2.2310-18 [0/90/02/90/0] 3.3010-9 1.1610-20 Table 4-2. Typical inte rsection areas for [0/90n/0] laminates. (unit: m2) Lay-ups [0/90/0] [0/902/0] [0/904/0] Case I: x = 1%, y = 0 8.3210-11 10.510-11 12.510-11 Case II: x = 0, y = 1% 7.7210-11 10.110-11 11.810-11 A B Figure 4-1. Three-dimensional delamina tion model and typical deformation. PAGE 54 54 A B Figure 4-2. COD along specific pa th shows the local effect. A B Figure 4-3. Variation of inte rsection area with respect to Nx, Ny, and T. COD along the crack front Distance along crack front (s/ h )COD/ h PAGE 55 55 A B Figure 4-4. J-integral dist ribution along crack front with different connect radius 0.E+00 2.E-04 4.E-04 6.E-04 8.E-04 1.E-03 1.E-03 1.E-03 2.E-03 2.E-03 2.E-03 00.511.522.533.54Radius (R/ h ) A / h2 load only in x-dir load only in y-dir thermal load Figure 4-5. Effect of co nnecting radius of delami nation on intersection area. s 0 20 40 60 80 100 120 140 160 051015202530Normalised distance (s/ h )J-integral, J/m2 rectanlge connection small radius circle large radius circle PAGE 56 56 0.E+00 5.E-04 1.E-03 2.E-03 2.E-03 3.E-03 00.511.52Delamination length ( a / h ) A / h2 load in x-dir load in y-dir thermal load Figure 4-6. Effect of delamina tion length on intersection area. A B Figure 4-7. Typical deformation of unit cell from different views. PAGE 57 57 Figure 4-8. Gas leakage path formed by the delamination from surface ply only. Figure 4-9. Comparison between two different scenarios. PAGE 58 58 A B Figure 4-10. Geometric model and typical deformation for multiple-layer laminates. 0.E+00 1.E-03 2.E-03 3.E-03 4.E-03 5.E-03 6.E-03 -250-200-150-100-500Thermal load T, oC A / h2 Temperature dependent Temperature independent Figure 4-11. Effect of temper ature dependent material prope rties on the intersection area. PAGE 59 59 Figure 4-12. Unit cell used to com pute typical intersection area. (x = y = 2 cm-1). PAGE 60 60 Figure 4-13. Typical deformation of the unit cell of [0/90/0]. Figure 4-14. Progressive permeability in laminates with loading ratio = 0.5. PAGE 61 61 Figure 4-15. Progressive permeability in laminates with loading ratio = 1.0. Figure 4-16. Progressive permeability in laminates with loading ratio = 2.0. PAGE 62 62 CHAPTER 5 DAMAGE AND PERMEABILITY PREDICTION IN MULTI-DIRECTIONAL LAMINATES In previous chapters, the permeability in cros s-ply laminates was investigated. Usually, it is good to include angle plies in realistic composite laminate design to gain better performance. The damage pattern in multi-directional laminates is quite different from the one in cross-ply laminates (Lavoie & Adolfsson, 2001; Yokozeki, Aoki, Ogasawara, & Ishikawa, 2005; Yokozeki & Ishikawa, 2005 and Yokozeki, Aoki, & Ishikawa, 2005). In th is chapter, the damage evolution and the permeability in multi-directional laminates is investigated. Laminates with stacking sequence [0/n/90]s (same as the laminates investigated by Yokozeki and his colleagues (Yokozeki, Aoki, Ogasawar a, & Ishikawa, 2005; Yokozeki & Ishikawa, 2005 and Yokozeki, Aoki, & Ishikawa, 2005)) are investigated in order to compare with their experimental observations. In their experiments, was taken as three different values (30o, 45o, and 60o). The specimen was loaded in 0o-direction. Depending on the value of and the thickness of o-ply, different patterns of cracking appe ar in the laminates. For large angle and thick angle ply, develope d cracks were generated in a ngle ply, and for small angle and thin angle ply, stitch cracks were generated. For example, developed cracks were observed in 30o-ply and stitch cracks were observed in 60o-ply in [0/2/90]s laminates. In the current work, = 30o and 60o are considered. Damage evol ution is first investigated and then gas permeability in laminate [0/60/90]s is predicted. Damage Prediction In the current model ply thickness is taken as 0.14 mm. The orthotropi c material properties are the same as those listed in Table 2-1. Figure 5-1 depicts the unit cell taken from the laminate. As shown in Figure 5-1A, the vertical lines represent the tr ansverse matrix cracks developed in the 90o-ply and the short-oblique lines represent the stitch cracks developed in the o-ply. In the PAGE 63 63 current model, it is assumed that both the transverse cracks in the 90o-ply and the stitch cracks in the o-ply are evenly distributed. Hence, a unit cell containing one stitch cr ack can be taken as shown in Figure 5-1B. It can be seen that once the ply th ickness is fixed this unit cell can be described using three parameters: the spacing of transverse crack lx, the spacing of stitch crack ly, and the stitch crack length 2a. Due to symmetric lay-up, only half of the unit cell is analyzed, as shown in Figure 5-1C. The dependence of strain energy release rate (SERR) on the three parameters (a, lx, and ly) is investigated in the current work. Two methods are used to compute the SERR: definition of SERR, and the J-integral at the cr ack front. In the first method, th e total strain energy in the unit cell is computed before and afte r stitch crack propagation, and then the SERR is obtained as the difference of the strain energy divided by the ge nerated crack surface area during propagation. In the second method, three-dimensional J-integrals along the stitch crack fr ont are computed. The computed J-integrals do not remain constant along the th ickness direction. Hence, the average of the J-integrals is taken as the SERR. Figure 5-2 shows the comparison of the two methods with the same unit cell. It can be seen that the SE RR computed from the two methods agree with each other quite well except for the first point. When the stitch crack is very short, there will be very few elements after the crack front, which may cause the error in the J-integral computation. In the following investigation, the J-integral me thod will be used in most cases, and the method using the definition of SERR will be used when the J-integral method is not applicable. Effects of Stitch Crack Length a First, spacing lx and ly are fixed and the variation of SE RR with respect to stitch crack length a is investigated as shown in Figure 5-3. Lay-ups [0/60/90]s and [0/30/90]s are considered. In the present model, periodic boundary conditions are applied to obtain unidirectional loading state: x = 0 and y = 0%. First, only mechanical loads are applied. The SERR computed from PAGE 64 64 J-integrals are plotted with s quare symbols. It can be seen that for a strain of about x = 1% the value of SERR is around 100 J/m2, which is much less than th e microcracking fracture toughness Gmc. It means that there is no crack formed at this strain, which is not true from experimental observations (Yokozeki, et al., 2005; Yokozeki & Ishikawa, 2005 and Yokozeki, Aoki, & Ishikawa, 2005). Hence the effect of thermal residual stress has to be included in the calculation. For these calculations, the un it cell is under zero strain (x = y = 0%) and the temperature decreases from the curing temperature (177oC) to room temperature (27oC). The computed SERR is denoted as GT. From mixed-mode fracture mechan ics (Appendix D) the SERR due to thermal residual stresses is include d in the total SERR using Equation D-8. It has been found that GT itself is very low (around 10-20 J/m2) but it has a great effect on the total SERR. From Figure 5-3, it can be seen that after in cluding thermal effects, the total SERR increases to around 250 J/m2, which is comparable to the microcracking fracture toughness Gmc. Hence, it is very important to take the thermal stress effect into account. From experimental observati ons by Yokozeki and his group (Yokozeki, et al., 2005; Yokozeki & Ishikawa, 2005 and Yokozeki, Aoki, & Ishikawa, 2005), for the same thickness of o-ply the stitch crack is more likely to develop in the laminates with greater values of It can be easily seen from Figure 5-3 that in lay-up [0/60/90]s the SERR drops down significantly after it reaches the maximum, while in [0/30/90]s the SERR drops down a little from the peak and then reaches a plateau. Hence, if the SERR decreased below the microcracking fracture toughness in [0/60/90]s, then the stitch crack would stop propagating. In [0/30/90]s laminate, if the plateau value is higher than the microcracking toughne ss, then the stitch crack would continue to propagate to form a developed crack. Further comparison is shown in Figure 5-4, where the [0/60/90]s laminate is under 1% strain and the [0/30/90]s is under 1.2% strain. It can be seen that PAGE 65 65 the SERR in [0/30/90]s requires more strain to reach th e same peak value as in [0/60/90]s, i.e., the onset strain of st itch crack in [0/30/90]s is greater than in [0/60/90]s. And once it reaches same value the SERR in [0/30/90]s remains higher than in [0/60/90]s. A typical value of microcracking fracture toughness Gmc (240 J/m2) for epoxy material has been taken (Nairn, 2000). It can be seen that under constant strain loading, the SERR in [0/60/90]s drops below Gmc at ac, which gives the stitch crack length at cu rrent strain state. Meanwhile, the SERR in [0/30/90]s is always greater than Gmc as a increases, which means that the stitch crack will continue to propagate to become a developed crack at this strain state. This confirms the results of experiments and analysis performed by Yokozekis group (Yokozeki, et al., 2005; Yokozeki & Ishikawa, 2005 and Yokozeki, Aoki, & Ishikawa, 2005). From Figure 5-4, there is a unique stitch crack length ac for a given loading state. Hence the stitch crack length can be represented in terms of the loading strain. Figure 5-5 shows the variation of stitch crack in [0//90]s laminate with respect to strain. It can be seen that the onset strain of stitch crack in [0/60/90]s is much less than in [0/30/90]s. Also, stitch crack propagation in [0/60/90]s is slower and st eady. In [0/30/90]s the stitch crack appears at high stain levels and then propagates very quickly to become a developed crack. Another phenomenon observed by Yokozeki et al. (2005) is th at the thickness of the o-ply affects the formation of stit ch cracks. The thicker the o-ply, the easier it is to form developed cracks instead of stitch cracks. Figures 5-6 and 5-7 show the effect of thickness of the o-ply on the SERR of stitch cracks. In Figure 5-6, the [0/60n/90]s laminate is investigated. In this figure, the line with diamond symbols repr esents the SERR variation with respect to stitch crack length in [0/60/90]s under x = 1% and the line with square sym bols represents the SERR variation in [0/602/90]s under x = 0.8%. It can be seen that even under lower strain le vels, the SERR in PAGE 66 66 [0/602/90]s is much greater than in [0/60/90]s and the stitch crack length ac is much greater in [0/602/90]s. For the laminates [0/30n/90]s, the same tendency can be found. As shown in Figure 5-7, the SERR in [0/302/90]s under 0.8% strain is hi gher than in [0/30/90]s under 1.2% strain. However, the fully developed crack s are present in both laminates. Effects of Spacing lx After the dependences of the SERR vari ation on the stitch crack length, angle and thickness of angle ply are inves tigated, the effects of spacing lx and ly will be investigated. First the variation of SERR with respect to lx in [0/60/90]s is computed for fixed a (0.6mm) and ly (1.0mm) as shown in Figure 5-8. J-integral method is used to compute the SERR. It can be seen that the SERR almost remains at the sa me level, which means that the spacing lx does not affect the SERR of stitch cracks so much. In this case, the prediction of the crack density x (=1/lx) in the 90o-ply and the stitch crack density in the angl e ply could be two i ndependent processes. Effects of Spacing ly Next, the SERR for generating new stitch cr acks is computed, i.e., the effect of ly on SERR is studied for a fixed a (0.6mm) and lx (2.5mm). Since we assume that the stitch cracks are evenly distributed, the new stitch crack has to ap pear at the middle of two existing cracks. Hence one unit cell becomes two, and the spacing ly reduces to half, as shown by Figure 5-9A. In this case, the old cracks will not propaga te and the released strain en ergy is due to newly generated cracks, hence the J-integral method is not vali d any more. So the SERR is computed using the strain energy difference method. Figure 5-9B shows the variation of the computed SERR, which tells us that the SERR remains almost as a constant for large ly and then drops down significantly as ly goes to zero. This means that once stitch cracks appear they will accumulate very fast until saturation occurs. PAGE 67 67 Permeability Prediction Model In Chapter 2, the gas permeability was derived for cross-ply laminates. Based on our assumption that the gas permeability is mainly de termined by the total intersection areas between layers, the permeation model can be easily ex tended to multi-directional laminates. The difference between multi-directional laminates and cross-ply laminates is in the method to compute the number of intersections and the intersection areas. As shown in Figure 5-10, instead of a rectangle, the intersection area formed by the crack openings of stitch cracks and transverse cracks in the 90o-ply is now a parallelogram Hence the intersection area A is given by d1d2/cos Under biaxial loading, cracks are generated in every layer. Compared to cross-ply laminates, there are more scenarios in multi-dire ctional laminates because of the crack formation in angle ply. In this work, we consider a special scenario as shown in Figure 5-11. In this scenario, we assume that the cracks in 0o-ply, 90o-ply and o-ply are all evenly distributed and characterized by respective crack densities. Stitch cracks are formed in angle ply and one stitch crack interact only once with the transverse cracks in 0o-ply and/or 90o-ply. We assume that there is always a stitch crack at the cro ss point between transverse cracks in 0o-ply and 90o-ply. In this case, only the stitch cracks at the cross poi nt will connect the tr ansverse cracks in 0o-ply and 90o-ply to form the contiguous gas leak pathway. In another word, the numbe r of intersections is only related to the cr ack densities of 0o-ply and 90o-ply. The crack densities in 0o-ply and 90o-ply can be predicted using the method given in Chapter 3. The only difference is that the angle ply with stitch cracks should be included into the unit cell. Based on the assumption that the crack density does not have much influence on the intersection area, for the comparison with the results in cross-ply laminate, the crack densities in PAGE 68 68 0o-ply and 90o-ply is taken as the same crack densities predicted for cross-ply laminates in Chapter 3. The stress ratio = 0.5 is considered, so the crack densities shown in Figure 3-8 are taken in the permeability computation. Laminate [0/60/90]s is taken as an example to demonstrate the prediction process. To compute the intersection area, a unit cell shown in Figure 5-12 is considered, whose dimensions are determined by the crack densities in the 0o-ply and 90o-ply. All the shaded areas are crack surfaces. Delaminations w ith a length of 0.14 mm from both 0o-plies and 90o-plies are introduced in this model (not shown in Figure 5-12). The crack density in the 0o-ply and 90o-ply are taken as 2/cm. From the aforementioned discu ssion, only the crack at the cross-point of the transverse cracks in 0o-ply and 90o-ply will be part of the contiguous leak pathway. So it is reasonable to only model the stitch crack at the cross-point. The results shown in Table 5-1 support this argument. In this tabl e, the intersection areas between 90o-ply and 60o-ply are computed for unit cells with different numbers of stitch cracks along transverse crack in 90o-ply. The unit cells are loaded under the same condition (x = 1%, y = 0). As one can see, as the number of the stitch cracks increase from 1 th rough 4, the intersection areas does not change so much. The difference is less than 1%. In Figure 5-13, the intersection areas are computed as a function of stitch crack length under the same loading condition (Case I: x = 1%, y = 0; Case II: x = 0, y = 1%). The stitch crack length (2a) increases from 1.2 mm to 2.8 mm. Un like the delamination length, the length of a stitch crack does not have much influence on the intersection areas. The intersection areas increase a little bit as the stitch crack length increases. But, when the stitch crack length increases by 133%, the increment of inters ection area is less than 5%. Under both loading cases, the PAGE 69 69 intersection areas between the 90o-ply and 60o-ply are larger than the areas between the 60o-ply and 0o-ply. Following the procedure to predict the progressi ve permeability in cross-ply laminates, the progressive normalized permeability in laminate [0/60/90]s is predicted as a function of applied stress in Figure 5-14. In this figure, the line with cro ss marks represents the permeability in [0/60/90]s laminate. For comparison, the permeability in [0/904/0] is re-plotted in Figure 5-14. It can be seen that the permeability in [0/60/90]s is much less than that in [0/904/0] even though they have same number of layers and same crack densities. At the same time, the permeability is plotted as a function of st rain as shown in Figure 5-15. In this figure, the permeabilities in both laminates are about the same for a given strain (< 0.55%). But, the permeability in [0/904/0] still has a greater increasing trend unde r higher strains. There are two reasons that can explain this. One is that the lay-up in [0/60/90]s is dispersed and in [0/904/0] is grouped. From the investigation in Chapter 4, dispersed lay-up ha s lower permeability than grouped lay-up. Another reason is that the intersection areas in [0/60/90]s is less than the intersection areas in [0/904/0] comparing the results in Table 5-1 and Table 4-2. In this section, a very idealized unit ce ll is taken to predict the permeability in multi-directional laminate [0/60/90]s. In reality, the case would be much more complicated and random. The stitch cracks may intera ct with the transverse cracks more than once, the position of the stitch cracks that connect the transverse cracks may be random, and the number of the intersection areas may be greater than in current model. Hence, some statistical tools may be needed to obtain a more realistic prediction. PAGE 70 70 Table 5-1. Intersection areas for unit cells with different numbers of stitch cracks. Number of stitch crac ks Intersection area (m2) 1 6.1920E-11 2 6.1769E-11 4 6.1505E-11 PAGE 71 71 Figure 5-1. Unit cell taken from laminate [0//90]s. A) Top view of lamina te. B) Top view of the unit cell. C) Three-dimensional view of the half of the unit cell. 0 20 40 60 80 100 120 140 160 180 00.20.40.60.811.2Stitch crack length a mmStrain energy release rate, J/m2 Method I: Definition of SERR Method II: J-integral at crack ft Figure 5-2. Comparison of strain energy releas e rate computed from two different methods. Stitch crack in o-ply Transverse crack in 90o-ply l y l x 2aA B C PAGE 72 72 [0/60/90]s0 50 100 150 200 250 300 00.20.40.60.811.2Stitch crack length a mmSERR, J/m2 Mechanical load only Mechanical and thermal load A [0/30/90]s0 50 100 150 200 250 300 35000.20.40.60.811.2Stitch crack length a mmSERR, J/m2 Mechanical load only Mechanical and thermal load B Figure 5-3. Thermal stress effect due to differenc e in curing temperature and service temperature. A) [0/60/90]s. B) [0/30/90]s. PAGE 73 73 0 50 100 150 200 250 300 350 00.20.40.60.811.2Stitch crack length a mmSERR, J/m2 [0/30/90]s [0/60/90]s ac Gmc Figure 5-4. Effect of angle on the formation of stitch crack. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.50.70.91.11.31.5x, %Stitch crack length ac, mm Figure 5-5. Progressive s titch crack growth in terms of strain in [0//90]s laminates. PAGE 74 74 0 50 100 150 200 250 300 350 400 450 00.20.40.60.811.2Stitch crack length a mmSERR, J/m2 n=2, strain_x=0.8% n=1, strain_x=1.0% a a Figure 5-6. Effect of nu mber of angle plies, n, of lay-up [0/60n/90]s on the strain energy release rate. 0 50 100 150 200 250 300 350 400 00.20.40.60.811.2Stitch crack length a mmSERR, J/m2 n = 1, strain_x = 1.2 % n = 2, strain_x = 0.8% Gmc Figure 5-7. Effect of nu mber of angle plies, n, of lay up [0/30n/90]s on the strain energy release rate. PAGE 75 75 0 50 100 150 200 250 300 1.51.71.92.12.32.52.72.9Spacing lx, mmSERR, J/m2 Figure 5-8. Effect of spacing lx on the strain en ergy release rate of stitch crack. A B Figure 5-9. Strain energy release rate rela ted to the formation of new stitch crack. 0 50 100 150 200 250 300 012345Spacing ly, mmSERR, J/m2 PAGE 76 76 Figure 5-10. Interlaminar cross-sectional area formed by crack openings of stitch crack and transverse crack. Figure 5-11. Special cracking scenario s in multi-directional laminates. d2 d1 PAGE 77 77 Figure 5-12. Unit cell for laminate with stitch cracks. 0.0E+00 2.0E-11 4.0E-11 6.0E-11 8.0E-11 1.0E-10 1.2E-10 1.4E-10 0.50.70.91.11.31.5Stitch Crack Length a, mmIntersection Area, m2 90/60, Case I 60/0, Case I 90/60, Case II 60/0, Case II Figure 5-13. Effect of the stitch crack length on the intersection areas. PAGE 78 78 Figure 5-14. Permeability in [0/60/90]s laminate with stitch crac ks versus average stress compared to the permeability in [0/904/0] laminate. Figure 5-15. Permeability in [0/60/90]s laminate with stitch cracks versus applied strain compared to the permeability in [0/904/0] laminate. PAGE 79 79 CHAPTER 6 CONCLUSION AND DISCUSSION Conclusions In this work, a finite element based method to predict gas permeability was used in conjunction with Darcys Law for porous materials. In Chapter 2, the permeability of cross-ply laminate was derived as a function of crack de nsity and the intersection area formed by crack opening displacements. Based on preliminary 2-D FE M modeling, the problem is decoupled into two parts: crack density predic tion and intersection area computat ion. In Chapter 3, the crack densities in both 0o-ply and 90o-ply in cross-ply laminate are predicted using a strain energy release rate method based on 3-D finite elemen t modeling. Then the intersection areas are investigated parametrically in Chapter 4 and the progressive permeability is predicted in terms of applied load. In Chapter 5, the damage evolutio n in multi-directional laminates is investigated and the permeability model is extended to multi-di rectional laminates. From this study, the following conclusions can be drawn: For cross-ply laminates [0/90n/0], the microcrack initiation stress and the crack density evolution largely depend on the thickness of the 90o-ply. Usually, the thinner the 90o-ply is, the greater the initiation stress and th e greater the crack ac cumulating rate after initiation. Local effects on the intersection area in the cross-ply are not negligible; The intersection area is str ongly influenced by the delami nation shape and delamination length; Leakage paths can still be formed in [0/90/0] laminates if delamination only occurs from the surface ply; For a given laminate thickness, the lamina tes with dispersed-layers show lower permeability than the laminates with grouped layers; Although some of the composite properti es depend strongly on the temperature, temperature independent materi al properties give the uppe r bound for intersection area. The permeability in cross-ply laminate is mainly dependent on crack densities. PAGE 80 80 In multi-directional laminates such as [0/n/90]s, whether stitch cracks or developed cracks will form in the o-ply depends on the angle and the thickness of the angle ply. The thicker angle ply and the smaller angle lead to longer cracks. Th e stitch crack length can be predicted using strain ener gy release rate criterion. Thermal effects can significantly increase the total strain energy release rate of stitch cracks. Crack density in the 90o-ply does not have much effect on the evolution of stitch cracks. The intersection area computed in the unit cell of a multi-directional laminate is not largely affected by the number of stitch cracks in the unit cell and the stitch crack length. The permeability predicted in multi-direction laminates is much lower than that in cross-ply laminate with same crack densities and number of layers. General Recommendations From the modeling results in the current study, we can make the following recommendations for the future desi gn of the composite cryogenic tank: Thin laminates are preferable to prevent microcrack initiation. For a given thickness, disp ersed lay-up performs better than grouped lay-up. Multi-directional laminate may have lower perm eability than cross-ply laminate with the same thickness. Future Work Based on the work done in this study, some future work can be suggested. First, the material constant involved in permeability calc ulation needs to be determined using some experimental techniques. Second, from what we found, the delaminati on length and the delamination shape affect the intersection area significantly. Hence, our model could predict the permeability in composite more accurately if the delamination length a nd the delamination shape can be measured experimentally. PAGE 81 81 Third, currently, an ideal unit cell is taken in the multi-direct ional laminates to simplify the prediction procedure, which could be improved in the future. Some statisti cal tools may be used to obtain the more realistic unit cell dimensions and the numbers of in tersection areas per unit area. Finally, a linear material consti tutive relationship was adopted in the FE analysis and a constant value of the microcrack ing fracture toughness is used in crack density prediction, which may not accurately characterize gr aphite/epoxy composite laminates. Hence, in the future, more realistic constitutive relations may be used, and the variability of the microcracking fracture toughness, e.g., as a function of te mperature, could be considered. PAGE 82 82 APPENDIX A SOME FLUID FLOW LAWS RE LATED TO PERMEABILITY Darcys Law Darcys Law is a phenomenological equation that describes the fl ow in porous materials. L P A Q (A-1) where Q is the volume flow rate, is the viscosity, is the permeability, A is the cross-section area, and P is the pressure drop along length L. Dividing both sides of Equation A-1 by cross-section area A yields more general form P q (A-2) where q is volume flow rate per unit area and P is pressure gradient. Poiseuille Law Poiseuilles Law (or Hagen-Poiseuille law) fo r laminar Newtonian fluid in a capillary is given as (Kulinchenko, 2005) 48r L P Q (A-3) where r is the diameter of the capillary. For the rectangular cross s ection with edge length a and b (a<=b), Poiseuilles Law take following form b a L P Q312 1 (A-4) It can be seen that volume flow rate is la rgely dependent on the si ze of the cross-section areas. Also, comparison between Darcys Law and Poiseuille Law implie s that the permeability is related to the square of th e size of the cross-section area. Darcy-Weisbach Equation Darcy-Weisbach Equation for pressure loss in pipe is (Kulinchenko, 2005) PAGE 83 83 hd L v P22 (A-5) where is density of fluid, v is the velocity of fluid, dh is hydraulic diameter of pipe, and is called friction factor which depe nds on the Reinolds Number Re. For the laminar flow (Re < 2300), the friction factor is given by Re 64 (A-6) Substituting Equation A-6 back to Equation A-5 leads to Poiseuilles Law (Equation A-3). PAGE 84 84 APPENDIX B MATLAB CODES FOR CRACK DENSITY PREDICTION Matlab Codes %*********************************************************************** % PROGRE SSIVE BIDIRECTIONA L MICROCRACKING. % This program predict crack densities for laminate [0/90/0] % The program first fits a complete cubic po lynomial to each of A11, A12 and A22 using % least squares approximation from data computed from FEM. Then it calculates % progressive bidirectional microcracking. % Laminate with ply thickness = 0.33mm and material properties from Table 2-1 % For other laminate, change the crack densities in the 'x' and 'y' matrices and % then change the three 'Z' matrices. % % Developed by: Satish Bapanapalli % Modified by: Jianlong Xu % % The unit used in this program is S.I. unit %*********************************************************************** clear all; % crack densities in x and y -direction x = [0 0.5 1 2 4 6 9.5]; y = [0 0.5 1 2 4 8 10 12.6]; % Cubic equation for A11. % Matrix Z is the A11 values obtained from FEM simulation for different combination of x and y value. % Then the A11 is fitted using least squire approximation as a function of crack density lambda_x and lambda_y Z=[115124166.300 115051076.169 114983796.644 114830976.957 114538885.788 114252612.330 113788986.829 115072076.088 115001789.291 114931532.927 114779254.329 114487140.722 114200061.974 113736878.189 115020420.990 114947917.088 114879789.965 114727684.948 114435504.389 114149182.967 113685477.863 114921525.458 114849519.101 114780990.642 114640564.895 114337532.488 114049833.473 113587040.365 114766384.882 114693253.386 114619226.459 114472602.251 114204007.518 113893268.851 113429330.391 114606404.382 114531177.475 114459259.306 114312180.932 114022654.509 113733427.377 113270133.621 114568529.731 114494422.741 114421387.135 114274040.267 113979649.641 113691711.390 113232333.069 114535304.276 114461680.841 114388158.047 114241385.079 113947459.433 113659520.466 113200147.736]; % Cubic_coef is a function to obtain the ten coefficients in the cubic polynomial using least square approximation. A11 = Cubic_coef(x, y, Z); % Cubic equation for A12 Z=[3043009.8959 3017132.5585 2993389.7655 2939459.3765 2836380.4873 2735356.0127 2571741.4473 2895607.0959 2870814.1059 2845982.6807 2793143.2410 2689950.8282 2588589.3957 2424984.7627 2750099.9526 2726547.0760 2700418.8413 2648510.5360 2545191.4487 2444046.3998 2280224.1904 2471521.6021 2449373.3120 2421934.0958 2372175.9896 2268559.7048 2166509.8680 2002994.4629 2034505.9375 2009605.6089 1982448.8061 1930523.6023 1831004.3452 1725733.2698 1563755.5200 1583854.7372 1553068.4196 1531865.5126 1479536.2481 1379434.9310 1274800.2837 1115510.3115 PAGE 85 85 1477165.5944 1449819.6792 1425196.6639 1377075.4985 1273097.6661 1171398.9266 1009129.8300 1383570.7846 1357575.4708 1331615.1250 1279790.8474 1182483.2131 1080813.3098 918597.3278 ]; A12 = Cubic_coef(x, y, Z); % Cubic equation for A22 Z=[61852537.9986 61843406.6387 61835027.8469 61815995.8661 61779618.8095 61743968.0339 61686227.2949 61437319.2629 61428601.5376 61419731.9961 61403619.5862 61366906.4832 61330987.6387 61272841.4494 61027439.9594 61024857.4708 61009760.9249 60996427.7670 60959380.1429 60923394.7747 60865077.7963 60252963.0372 60244087.9876 60225235.3082 60207108.2724 60179906.4239 60142423.7126 60084173.7244 59011682.3750 59005453.5207 58992962.0937 58974127.9208 58925680.0016 58900890.0190 58848651.0728 57742242.2748 57719438.2380 57723723.0053 57704059.4061 57674572.4809 57630416.7839 57586055.8405 57441705.2035 57428693.0542 57423247.8372 57417188.6166 57380249.8881 57344116.7740 57286422.5137 57178059.2794 57168846.8293 57159646.4955 57141279.8226 57125021.0986 57088979.9209 57031447.9881]; A22 = Cubic_coef(x, y, Z); % Input the initial crack densities and the ratio between Ny and Nx lambda_x = input('\nInitial crack density in the x-direction (cracks/cm) = '); lambda_y = input('\nInitial crack density in the y-direction (cracks/cm)= '); Ny_limit = input('\nUpper limit for Ny (N*m) = '); fid = fopen('results.txt','a'); alpha = input('\nALPHA: Ratio of Nx to Ny ='); % Microcracking fracture t oughness and ply toughness Gmc = 240; t1 = 0.00033; % % Create the current [ A ] matr ix from the cubic functions A1(1,1) = A11(1)*lambda_x^3 + A11(2)*lambda_x^2*lambda_y + A11(3)*lambda_x*lambda_y^2 + A11(4)*lambda_y^3 + A11(5)*lambda_x^2 + A11(6)*lambda_x*lambda_y + A11(7)*lambda_y^2 + A11(8)*lambda_x + A11(9)*lambda_y + A11(10); A1(1,2) = A12(1)*lambda_x^3 + A12(2)*lambda_x^2*lambda_y + A12(3)*lambda_x*lambda_y^2 + A12(4)*lambda_y^3 + A12(5)*lambda_x^2 + A12(6)*lambda_x*lambda_y + A12(7)*lambda_y^2 + A12(8)*lambda_x + A12(9)*lambda_y + A12(10); A1(2,1) = A1(1,2); A1(2,2) = A22(1)*lambda_x^3 + A22(2)*lambda_x^2*lambda_y + A22(3)*lambda_x*lambda_y^2 + A22(4)*lambda_y^3 + A22(5)*lambda_x^2 + A22(6)*lambda_x*lambda_y + A22(7)*lambda_y^2 + A22(8)*lambda_x + A22(9)*lambda_y + A22(10); A1bar = inv(A1); % fprintf('\n\nSTEP lambda_x lambda_y Ny_1 Ny_2 Ny_3\n'); fprintf(fid, '\n\n lambda_x0 = %e, lambda_y0 = %e', lambda_x, lambda_y); % Start the loop step = 1; Ny = 0.0; for i = 1:1:20 lambda_x2 = 2*lambda_x; lambda_y2 = 3*lambda_y; len_x = 0.01 / lambda_x; len_y = 0.01 / lambda_y; PAGE 86 86 % lambda is in 1/cm, len is in m if (lambda_x2 > 9.5) fprintf('\nCRACK DENSITY LIMIT: x-DIRECTION'); break; end if (lambda_y2 > 12.6) fprintf('\nCRACK DENSITY LIMIT: y-DIRECTION'); break; end % % Case I: new crack in the surface ply. % Create the [A] matrix for Case I. A2_1(1,1) = A11(1)*lambda_x^3 + A11(2)*lambda_x^2*lambda_y2 + A11(3)*lambda_x*lambda_y2^2 + A11(4)*lambda_y2^3 + A11(5)*lambda_x^2 + A11(6)*lambda_x*lambda_y2 + A11(7)*lambda_y2^2 + A11(8)*lambda_x + A11(9)*lambda_y2 + A11(10); A2_1(1,2) = A12(1)*lambda_x^3 + A12(2)*lambda_x^2*lambda_y2 + A12(3)*lambda_x*lambda_y2^2 + A12(4)*lambda_y2^3 + A12(5)*lambda_x^2 + A12(6)*lambda_x*lambda_y2 + A12(7)*lambda_y2^2 + A12(8)*lambda_x + A12(9)*lambda_y2 + A12(10); A2_1(2,1) = A2_1(1,2); A2_1(2,2) = A22(1)*lambda_x^3 + A22(2)*lambda_x^2*lambda_y2 + A22(3)*lambda_x*lambda_y2^2 + A22(4)*lambda_y2^3 + A22(5)*lambda_x^2 + A22(6)*lambda_x*lambda_y2 + A22(7)*lambda_y2^2 + A22(8)*lambda_x + A22(9)*lambda_y2 + A22(10); A2_1bar = inv(A2_1); % %Calculate Ny Ny_1 = sqrt(8*Gmc*t1/((alpha^2*(A2_1bar(1,1)-A1bar(1,1))+2*alpha*(A2_1bar(1,2)-A1bar(1,2))+(A2_1bar(2,2)-A1bar(2, 2)))*len_y)); % % Case II: new crack in the middle ply. % Create the [A] matrix for Case II. A2_2(1,1) = A11(1)*lambda_x2^3 + A11(2)*lambda_x2^2*lambda_y + A11(3)*lambda_x2*lambda_y^2 + A11(4)*lambda_y^3 + A11(5)*lambda_x2^2 + A11(6)*lambda_x2*lambda_y + A11(7)*lambda_y^2 + A11(8)*lambda_x2 + A11(9)*lambda_y + A11(10); A2_2(1,2) = A12(1)*lambda_x2^3 + A12(2)*lambda_x2^2*lambda_y + A12(3)*lambda_x2*lambda_y^2 + A12(4)*lambda_y^3 + A12(5)*lambda_x2^2 + A12(6)*lambda_x2*lambda_y + A12(7)*lambda_y^2 + A12(8)*lambda_x2 + A12(9)*lambda_y + A12(10); A2_2(2,1) = A2_2(1,2); A2_2(2,2) = A22(1)*lambda_x2^3 + A22(2)*lambda_x2^2*lambda_y + A22(3)*lambda_x2*lambda_y^2 + A22(4)*lambda_y^3 + A22(5)*lambda_x2^2 + A22(6)*lambda_x2*lambda_y + A22(7)*lambda_y^2 + A22(8)*lambda_x2 + A22(9)*lambda_y + A22(10); A2_2bar = inv(A2_2); Ny_2 = sqrt(2*Gmc*t1/((alpha^2*(A2_2bar(1,1)-A1bar(1,1))+2*alpha*(A2_2bar(1,2)-A1bar(1,2))+(A2_2bar(2,2)-A1bar(2, 2)))*len_x)); % PAGE 87 87 % Case III: new crack in both the middle and the surface ply. % Create the [A] matrix for Case III. A2_3(1,1) = A11(1)*lambda_x2^3+A11(2)*lambda_x2^2*lambda_y2 + A11(3)*lambda_x2*lambda_y2^2 + A11(4)*lambda_y2^3 + A11(5)*lambda_x2^2 + A11(6)*lambda_x2*lambda_y2 + A11(7)*lambda_y2^2 + A11(8)*lambda_x2 + A11(9)*lambda_y2 + A11(10); A2_3(1,2) = A12(1)*lambda_x2^3+A12(2)*lambda_x2^2*lambda_y2 + A12(3)*lambda_x2*lambda_y2^2 + A12(4)*lambda_y2^3 + A12(5)*lambda_x2^2 + A12(6)*lambda_x2*lambda_y2 + A12(7)*lambda_y2^2 + A12(8)*lambda_x2 + A12(9)*lambda_y2 + A12(10); A2_3(2,1) = A2_3(1,2); A2_3(2,2) = A22(1)*lambda_x2^3+A22(2)*lambda_x2^2*lambda_y2 + A22(3)*lambda_x2*lambda_y2^2 + A22(4)*lambda_y2^3 + A22(5)*lambda_x2^2 + A22(6)*lambda_x2*lambda_y2 + A22(7)*lambda_y2^2 + A22(8)*lambda_x2 + A22(9)*lambda_y2 + A22(10); A2_3bar = inv(A2_3); Ny_3 = sqrt((2*Gmc*len_y*t1 + 8*Gmc*len_x*t1)/((alpha^2*(A2_3bar(1,1)-A1bar(1,1))+2*alpha*(A2_3bar(1,2)-A1bar(1,2))+(A2_3bar(2,2)-A1ba r(2,2)))*len_x*len_y)); % %fprintf('\n%-8d %-10.2f %-10.2f %-15.2f %-15.2f %-15.2f', step, lambda_x, lambda_y, Ny_1, Ny_2, Ny_3); % % Compare the 3 values of Ny and find the least amongst them. if (Ny_1 < Ny_2) if (Ny_1 < Ny_3) least = 1; else least = 3; end; elseif (Ny_2 < Ny_3) least = 2; else least = 3; end; % % Prepare for next step if (least == 1) Ny = Ny_1; A1 = A2_1; A1bar = inv(A1); lambda_y = lambda_y2; elseif (least == 2) Ny = Ny_2; A1 = A2_2; A1bar = inv(A1); lambda_x = lambda_x2; else Ny = Ny_3; A1 = A2_3; A1bar = inv(A1); lambda_x = lambda_x2; lambda_y = lambda_y2; end; PAGE 88 88 % % Print this step summary fprintf('\n%-8d %-10.2f %-10.2f %-15.2f %-15.2f %-15.2f %-15.2f', step, lambda_x, lambda_y, Ny_1, Ny_2, Ny_3, Ny); fprintf(fid, '\n%-8d %-15.2f %-10.3f %-10.3f', step, Ny, lambda_x, lambda_y); step = step + 1; if (Ny >= Ny_limit) fprintf('\n LOAD LIMIT REACHED.\n'); break; end; end; % for loop ends. fprintf(fid,'\n'); fclose(fid); fprintf('\n END PROGRAM.\n'); function [A] = Cubic_coef(x, y, Z) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Get the ten coefficients for a complete cubic polynomial % % Input: x first variable (N by 1) % y second variable (M by 1) % Z values of function (N by M) % % Output: A matrix of ten coefficients %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C = zeros(10,10); R = zeros(10,1); for j = 1:1:7 for i = 1:1:8 C(1,1)=C(1,1)+x(j)^6; C(1,2)=C(1,2)+x(j)^5*y(i); C(1,3)=C(1,3)+x(j)^4*y(i)^2; C(1,4)=C(1,4)+x(j)^3*y(i)^3; C(1,5)=C(1,5)+x(j)^5; C(1,6)=C(1,6)+x(j)^4*y(i); C(1,7)=C(1,7)+x(j)^3*y(i)^2; C(1,8)=C(1,8)+x(j)^4; C(1,9)=C(1,9)+x(j)^3*y(i); C(1,10)=C(1,10)+x(j)^3; C(2,1)=C(1,2); C(2,2)=C(1,3); C(2,3)=C(1,4); C(2,4)=C(2,4)+x(j)^2*y(i)^4; C(2,5)=C(1,6); C(2,6)=C(1,7); C(2,7)=C(2,7)+x(j)^2*y(i)^3; C(2,8)=C(1,9); C(2,9)=C(2,9)+x(j)^2*y(i)^2;C(2,10)=C(2,10)+x(j)^2*y(i); C(3,1)=C(2,2); C(3,2)=C(2,3); C(3,3)=C(2,4); C(3,4)=C(3,4)+x(j)*y(i)^5; C(3,5)=C(2,6); C(3,6)=C(2,7); C(3,7)=C(3,7)+x(j)*y(i)^4; C(3,8)=C(2,9); C(3,9)=C(3,9)+x(j)*y(i)^3; C(3,10)=C(3,10)+x(j)*y(i)^2; C(4,1)=C(3,2); C(4,2)=C(3,3); C(4,3)=C(3,4); C(4,4)=C(4,4)+y(i)^6; C(4,5)=C(3,6); C(4,6)=C(3,7); C(4,7)=C(4,7)+y(i)^5; C(4,8)=C(3,9); C(4,9)=C(4,9)+y(i)^4; C(4,10)=C(4,10)+y(i)^3; C(5,1)=C(1,5); C(5,2)=C(2 ,5); C(5,3)=C(3,5); C(5,4)=C(4,5); C(5,5)=C(1,8); C(5,6)=C(2,8); C(5,7)=C(3,8); C(5,8)=C(1,10); C(5,9)=C(2,10); C(5,10)=C(5,10)+x(j)^2; C(6,1)=C(5,2); C(6,2)=C(5 ,3); C(6,3)=C(5,4); C(6,4)=C(4,6); C(6,5)=C(5,6); C(6,6)=C(5,7); C(6,7)=C(4,8); C(6,8)=C(5,9); C(6,9)=C(3,10); C(6,10)=C(6,10)+x(j)*y(i); C(7,1)=C(6,2); C(7,2)=C(6 ,3); C(7,3)=C(6,4); C(7,4)=C(4,7); C(7,5)=C(6,6); C(7,6)=C(6,7); C(7,7)=C(4,9); C(7,8)=C(6,9); C(7,9)=C(4,10); C(7,10)=C(7,10)+y(i)^2; C(8,1)=C(1,8); C(8,2)=C(2 ,8); C(8,3)=C(3,8); C(8,4)=C(3,9); C(8,5)=C(5,8); C(8,6)=C(6,8); C(8,7)=C(7,8); C(8,8)=C(5,10); C(8,9)=C(6,10); C(8,10)=C(8,10)+x(j); C(9,1)=C(1,9); C(9,2)=C(2 ,9); C(9,3)=C(3,9); C(9,4)=C(4,9); C(9,5)=C(5,9); C(9,6)=C(6,9); C(9,7)=C(7,9); C(9,8)=C(8,9); C(9,9)=C(7,10); C(9,10)=C(9,10)+y(i); C(10,1)=C(1,10); C(10,2)=C(2,10); C(10,3)= C(3,10); C(10,4)=C(4,10); C(10,5)=C(5,10); C(10,6)=C(6,10); C(10,7)=C(7,10); C(10,8)=C(8,10); C(10,9)=C(9,10); C(10,10)=C(10,10)+1; PAGE 89 89 R(1)=R(1)+Z(i,j)*x(j)^3; R(2)=R(2)+Z(i,j)*x(j)^2*y(i); R(3)=R(3)+Z(i,j)*x(j)*y(i)^2; R(4)=R(4)+Z(i,j)*y(i)^3; R(5)=R(5)+Z(i,j)*x(j)^2; R(6)=R(6)+Z(i,j)*x(j)*y(i); R(7)=R(7)+Z(i,j)*y(i)^2; R(8)=R(8)+Z(i,j)*x(j); R(9)=R(9)+Z(i,j)*y(i); R(10)=R(10)+Z(i,j); end; end; A = inv(C) R; PAGE 90 90 Least Square Approximation for [A] Matrix Suppose the computed values for [A] matrix are available for the crack density combinations of m xs and n ys, then for each entry in [A] matrix, there are mn values. Take A11 as an example m x x x x 2 1 (B-1) n y y y y 2 1 (B-2) The computed values of A11 from finite element analysis are mn m m n nZ Z Z Z Z Z Z Z Z Z 2 1 2 22 21 1 12 11 (B-3) Suppose A11 is a complete cubic polynomial of x and y: 10 9 8 2 7 6 2 5 3 4 2 3 2 2 3 1 11a a a a a a a a a a Ay x y y x x y y x y x x (B-4) where a1 to a10 are the coefficients need to be determined. Then the error can be computed as m i n j j y j y i x j y i x i x ija a a a a Z E1 2 10 3 4 2 3 2 2 3 1 (B-5) Want the error to be minimum, so m i n j j y j y i x j y i x i x ij m i n j j y i x j y j y i x j y i x i x ij m i n j i x j y j y i x j y i x i x ij ia a a a a Z a a a a a Z a a a a a Z a E11 10 3 4 2 3 2 2 3 1 11 2 10 3 4 2 3 2 2 3 1 11 3 10 3 4 2 3 2 2 3 10 0 0 0 (B-6) Writing this system of equation into matrix format yield 1 10 1 10 10 10 R a C (B-7) where PAGE 91 91 10 2 1a a a a ij ij ij j y i x ij ij i x ijZ Z Z R 2 3, and ij ij j y i x ij j y i x ij j y i x ij j y i x ij j y i x ij i x ij j y i x ij j y i x ij i xSym C1 .2 4 2 2 3 3 2 4 3 2 4 5 6 Solving the linear system given by Equation B-7 will give all ten coefficients needed in the cubic polynomial for A11. The same procedures can be applied to A12 and A22 as well. PAGE 92 92 APPENDIX C MATLAB CODES TO PREDICT PR OGRESSIVE PERMEABILITY Functions Amatrix_n1, Amatrix_n2, and Amatrix_n4 are used in the main program later. function [A_matrix] = Amatrix_n1(lambda_x, lambda_y) %************************************************************************** % This function is to compute the A matrix for different combination of % lambda_x and lambda_y for the laminate [0/90/0]. % % Input: lambda_x, lambda_y (/cm) % Output: A matrix (4X4) % Function Cubic_coef is listed in Appendix B. %************************************************************************** x = [0 0.5 1 2 4 6 9.5]; y = [0 0.5 1 2 4 8 10 12.6]; % Cubic equation for A11 Z = [115124166.300 115051076.169 114983796.644 114830976.957 114538885.788 114252612.330 113788986.829 115072076.088 115001789.291 114931532.927 114779254.329 114487140.722 114200061.974 113736878.189 115020420.990 114947917.088 114879789.965 114727684.948 114435504.389 114149182.967 113685477.863 114921525.458 114849519.101 114780990.642 114640564.895 114337532.488 114049833.473 113587040.365 114766384.882 114693253.386 114619226.459 114472602.251 114204007.518 113893268.851 113429330.391 114606404.382 114531177.475 114459259.306 114312180.932 114022654.509 113733427.377 113270133.621 114568529.731 114494422.741 114421387.135 114274040.267 113979649.641 113691711.390 113232333.069 114535304.276 114461680.841 114388158.047 114241385.079 113947459.433 113659520.466 113200147.736]; A11 = Cubic_coef(x, y, Z); % % Cubic equation for A12 Z=[3043009.8959 3017132.5585 2993389.7655 2939459.3765 2836380.4873 2735356.0127 2571741.4473 2895607.0959 2870814.1059 2845982.6807 2793143.2410 2689950.8282 2588589.3957 2424984.7627 2750099.9526 2726547.0760 2700418.8413 2648510.5360 2545191.4487 2444046.3998 2280224.1904 2471521.6021 2449373.3120 2421934.0958 2372175.9896 2268559.7048 2166509.8680 2002994.4629 2034505.9375 2009605.6089 1982448.8061 1930523.6023 1831004.3452 1725733.2698 1563755.5200 1583854.7372 1553068.4196 1531865.5126 1479536.2481 1379434.9310 1274800.2837 1115510.3115 1477165.5944 1449819.6792 1425196.6639 1377075.4985 1273097.6661 1171398.9266 1009129.8300 1383570.7846 1357575.4708 1331615.1250 1279790.8474 1182483.2131 1080813.3098 918597.3278 ]; A12 = Cubic_coef(x, y, Z); % % Cubic equation for A22 Z=[61852537.9986 61843406.6387 61835027.8469 61815995.8661 61779618.8095 61743968.0339 61686227.2949 61437319.2629 61428601.5376 61419731.9961 61403619.5862 61366906.4832 61330987.6387 61272841.4494 61027439.9594 61024857.4708 61009760.9249 60996427.7670 60959380.1429 60923394.7747 60865077.7963 60252963.0372 60244087.9876 60225235.3082 60207108.2724 60179906.4239 60142423.7126 60084173.7244 PAGE 93 93 59011682.3750 59005453.5207 58992962.0937 58974127.9208 58925680.0016 58900890.0190 58848651.0728 57742242.2748 57719438.2380 57723723.0053 57704059.4061 57674572.4809 57630416.7839 57586055.8405 57441705.2035 57428693.0542 57423247.8372 57417188.6166 57380249.8881 57344116.7740 57286422.5137 57178059.2794 57168846.8293 57159646.4955 57141279.8226 57125021.0986 57088979.9209 57031447.9881]; A22 = Cubic_coef(x, y, Z); % % Create the current [ A ] matrix A_matrix = zeros(2,2); A_matrix(1,1) = A11(1)*lambda_x^3 + A11(2)*lambda_x^2*lambda_y + A11(3)*lambda_x*lambda_y^2 + A11(4)*lambda_y^3 + A11(5)*lambda_x^2 + A11(6)*lambda_x*lambda_y + A11(7)*lambda_y^2 + A11(8)*lambda_x + A11(9)*lambda_y + A11(10); A_matrix(1,2) = A12(1)*lambda_x^3 + A12(2)*lambda_x^2*lambda_y + A12(3)*lambda_x*lambda_y^2 + A12(4)*lambda_y^3 + A12(5)*lambda_x^2 + A12(6)*lambda_x*lambda_y + A12(7)*lambda_y^2 + A12(8)*lambda_x + A12(9)*lambda_y + A12(10); A_matrix(2,1) = A_matrix(1,2); A_matrix(2,2) = A22(1)*lambda_x^3 + A22(2)*lambda_x^2*lambda_y + A22(3)*lambda_x*lambda_y^2 + A22(4)*lambda_y^3 + A22(5)*lambda_x^2 + A22(6)*lambda_x*lambda_y + A22(7)*lambda_y^2 + A22(8)*lambda_x + A22(9)*lambda_y + A22(10); function A_matrix = Amatrix_n2(lambda_x, lambda_y) %************************************************************************** % This function is to compute the A matrix for different combination of % lambda_x and lambda_y for the laminate [0/90_2/0]. % % Input: lambda_x, lambda_y (/cm) % Output: A matrix (4X4) %************************************************************************** x = [0 0.5 1 2 4 6 9.5]; y = [0 0.5 1 2 4 8 10 12.6]; % Cubic equation for A11 Z = [117987202.37 117687685.17 117389249.80 116795817.32 115664529.51 114765880.45 113848508.31 117941618.09 117642126.80 117343636.32 116751802.12 115617652.23 114718967.31 113781643.83 117896306.71 117596545.42 117298299.63 116704668.88 115573290.73 114674561.91 113736077.79 117807051.38 117507160.86 117208181.20 116615470.79 115484007.34 114585203.44 113646618.19 117650638.19 117350412.42 117051715.53 116458961.59 115327417.07 114428548.97 113505654.61 117477822.48 117177546.83 116878843.86 116285986.90 115156132.08 114262635.21 113339153.63 117434051.83 117133789.81 116835059.39 116244529.51 115114618.48 114220867.69 113297837.96 117401559.94 117101227.67 116802427.27 116209509.48 115079453.48 114186189.03 113263112.80 ]; A11 = Cubic_coef(x, y, Z); % % Cubic equation for A12 Z = [4053287.4987 3947586.8143 3842268.7924 3632846.3491 3233614.7176 2916482.5496 2592740.8214 3924871.3483 3819266.3069 3713804.1817 3504906.9703 3104587.7396 2787215.2078 2456267.5829 3797234.0365 3691409.3130 3586119.4655 3376178.4570 2976760.3371 2659451.6819 2328029.6941 3545811.4744 3440152.5682 3334320.0401 3125004.0099 2725414.9172 2407941.0838 2076263.1169 3105210.1150 2999169.6788 2893669.0654 2684306.8171 2284627.6703 1967063.5446 1641732.8102 2618406.8082 2512385.9593 2406920.4852 2197595.5123 1798659.9210 1483137.2799 1159552.7446 2495109.0042 2389106.5224 2283644.7407 2081399.9765 1682500.0854 1366943.1636 1043567.2000 2403579.8806 2297567.5454 2192095.4843 1982804.7475 1583908.4180 1268578.1089 946272.8516 ]; A12 = Cubic_coef(x, y, Z); % % Cubic equation for A22 PAGE 94 94 Z = [117979056.58 117941755.98 117904588.47 117830683.21 117689796.28 117577879.52 117463630.49 117617321.86 117580294.89 117542726.62 117468898.10 117327408.12 117214766.65 117097511.80 117257780.04 117220326.26 117183061.64 117107711.46 116966323.50 116853929.04 116736325.35 116549547.72 116512749.94 116474495.85 116400212.88 116258368.98 116145527.28 116027214.60 115308416.29 115270744.37 115233262.23 115158877.27 115016829.15 114903767.56 114789971.37 113937140.24 113899578.02 113862212.51 113788047.45 113646675.82 113534743.97 113427066.78 113589813.46 113552296.52 113514969.80 113460932.98 113319730.54 113207942.00 113100476.92 113331986.95 113294507.36 113257217.03 113183219.80 113042172.88 112930618.35 112826593.40 ]; A22 = Cubic_coef(x, y, Z); % % Create the current [ A ] matrix A_matrix = zeros(2,2); A_matrix(1,1) = A11(1)*lambda_x^3 + A11(2)*lambda_x^2*lambda_y + A11(3)*lambda_x*lambda_y^2 + A11(4)*lambda_y^3 + A11(5)*lambda_x^2 + A11(6)*lambda_x*lambda_y + A11(7)*lambda_y^2 + A11(8)*lambda_x + A11(9)*lambda_y + A11(10); A_matrix(1,2) = A12(1)*lambda_x^3 + A12(2)*lambda_x^2*lambda_y + A12(3)*lambda_x*lambda_y^2 + A12(4)*lambda_y^3 + A12(5)*lambda_x^2 + A12(6)*lambda_x*lambda_y + A12(7)*lambda_y^2 + A12(8)*lambda_x + A12(9)*lambda_y + A12(10); A_matrix(2,1) = A_matrix(1,2); A_matrix(2,2) = A22(1)*lambda_x^3 + A22(2)*lambda_x^2*lambda_y + A22(3)*lambda_x*lambda_y^2 + A22(4)*lambda_y^3 + A22(5)*lambda_x^2 + A22(6)*lambda_x*lambda_y + A22(7)*lambda_y^2 + A22(8)*lambda_x + A22(9)*lambda_y + A22(10); function A_matrix = Amatrix_n4(lambda_x, lambda_y) %************************************************************************** % This function is to compute the A matrix for different combination of % lambda_x and lambda_y for the laminate [0/90_4/0]. % % Input: lambda_x, lambda_y (/cm) % Output: A matrix (4X4) %************************************************************************** x = [0 0.5 1 2 4 6 9.5]; y = [0 0.5 1 2 4 8 10 12.6]; % Cubic equation for A11 Z = [123712759.73 122494744.24 121297450.60 119038610.88 116023569.18 114705419.21 113785266.84 123669288.74 122451361.56 121254657.92 118996636.75 115973574.90 114661655.69 113718223.70 123625952.46 122406677.38 121210501.44 118953188.95 115936262.77 114618039.67 113674574.27 123539670.48 122320030.72 121123178.60 118865082.87 115849801.05 114531509.07 113602596.01 123384655.15 122164966.14 120967580.98 118709754.04 115694346.57 114376010.47 113447311.49 123205464.39 121984694.65 120788096.83 118530103.50 115526392.29 114211221.65 113278207.24 123161668.09 121940828.47 120744172.14 118486089.06 115483480.72 114168163.28 113235476.85 123127912.91 121907512.26 120711288.31 11 8455106.86 115447337.37 114132440.23 113190438.82 ]; A11 = Cubic_coef(x, y, Z); % % Cubic equation for A12 Z = [6073842.5182 5644003.6413 5221478.4213 4424333.0737 3360324.0431 2895150.1467 2570427.5831 5951378.9434 5521864.8269 5099295.8205 4302290.0779 3235347.5148 2772095.4604 2439032.2755 5829302.9700 5399231.3841 4976779.3288 4180066.3604 3114834.4295 2649453.5407 2316254.3357 5586256.0118 5156190.0921 4733449.2430 3936199.3370 2871721.7207 2406143.5528 2077463.2805 5149592.8446 4719652.1564 4296355.7369 3499325.2493 2434703.1371 1968993.7253 1640732.5978 4644832.4978 4213938.5267 3791576.2978 2994568.9492 1936441.4063 1472009.8981 1144726.7438 4521461.0644 4090571.6296 3668217.2933 2871234.1053 1816189.9701 1351810.6800 1025007.0641 4426377.3502 3995666.8782 3573488.8468 2777223.0115 1715663.9635 1251522.1743 920068.8042 ]; A12 = Cubic_coef(x, y, Z); % % Cubic equation for A22 PAGE 95 95 Z = [230232128.00 230080439.85 229931330.67 229650017.85 229274529.80 229110369.49 228995775.35 229887161.45 229736413.01 229586580.22 229304895.60 228928080.98 228763820.03 228645930.72 229543286.64 229392112.42 229242127.24 228960666.31 228583176.00 228418433.73 228300146.88 228858646.85 228707858.77 228557633.20 228275249.25 227898515.39 227733217.55 227614786.31 227628610.90 227478272.57 227326817.56 227044863.59 226667797.04 226502129.16 226384812.14 226206750.23 226054463.16 225905187.18 225623475.89 225255360.90 225090596.43 224980665.49 225859208.72 225707000.42 225557805.12 225276257.61 224916825.38 224752372.53 224643635.55 225591371.96 225439292.89 225290224.48 225009059.31 224634109.76 224469991.25 224355784.54 ]; A22 = Cubic_coef(x, y, Z); % % Create the current [ A ] matrix A_matrix = zeros(2,2); A_matrix(1,1) = A11(1)*lambda_x^3 + A11(2)*lambda_x^2*lambda_y + A11(3)*lambda_x*lambda_y^2 + A11(4)*lambda_y^3 + A11(5)*lambda_x^2 + A11(6)*lambda_x*lambda_y + A11(7)*lambda_y^2 + .. A11(8)*lambda_x + A11(9)*lambda_y + A11(10); A_matrix(1,2) = A12(1)*lambda_x^3 + A12(2)*lambda_x^2*lambda_y + A12(3)*lambda_x*lambda_y^2 + A12(4)*lambda_y^3 + A12(5)*lambda_x^2 + A12(6)*lambda_x*lambda_y + A12(7)*lambda_y^2 + A12(8)*lambda_x + A12(9)*lambda_y + A12(10); A_matrix(2,1) = A_matrix(1,2); A_matrix(2,2) = A22(1)*lambda_x^3 + A22(2)*lambda_x^2*lambda_y + A22(3)*lambda_x*lambda_y^2 + A22(4)*lambda_y^3 + A22(5)*lambda_x^2 + A22(6)*lambda_x*lambda_y + A22(7)*lambda_y^2 + A22(8)*lambda_x + A22(9)*lambda_y + A22(10); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % This prgram is to plot the norma lized permeability with respect to the % force resultant Ny(or average stress sigma_y) with Nx be proportional to % Ny and neglecting the thermal stress effect. % case I: epsilon_x = 0, epsilon_y = 1%, dT = 0 % case II: epsilon_x = 1%, epsilon_y = 0, dT = 0 % The basic idea is to get epsilon_x, epsilon_y from given N_x, N_y % and then superpose above two cases to get u v and the intersection area. % The normalized permeability is co mputed according to Equation 2-11. % In this program, the crack density pr edicted from FE analysis are used. % The ratio between Nx and Ny: ratio is input and the permeability for all three laminates % are plotted in the same figure as a function of Ny. % In this program, dT = 0. % % Developed by: Jianlong Xu, May, 2007. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all; % ply thickness t_ply = 0.00033; % The ratio between N_x and N_y ratio = 0.5; % Make a choice for three different lay-ups display('what kind of lay-up?'); display('[0/90/0], [0/90_2/0] or [0/90_4/0]?'); n = input('make your choice (1, 2 or 3):'); % Laminate thickness for different lay-up switch n case 1 t = 3*t_ply; case 2 PAGE 96 96 t = 4*t_ply; case 3 t = 6*t_ply; otherwise display('wrong choice'); end % Typical crack opening displacements for different lay-ups. % The crack opening displacement is obtained from unit cell where lx=ly=2.0/cm % % Case I: epsilon_x = 1%, epsilon_y = 0, dT = 0 epsilon_x1 = 0.01 switch n case 1 dU1 = 6.98E-06; dV1 = 2.98E-06; case 2 dU1 = 8.96E-06; dV1 = 2.94E-06; case 3 dU1 = 1.06E-05; dV1 = 2.97E-06; otherwise display('wrong choice'); end % Case II: epsilon_x = 0, epsilon_y = 1%, dT = 0 epsilon_y2 = 0.01 switch n case 1 dU2 = 2.41E-06; dV2 = 8.01E-06; case 2 dU2 = 3.28E-06; dV2 = 7.68E-06; case 3 dU2 = 3.91E-06; dV2 = 7.55E-06; otherwise display('wrong choice'); end % Stress range is determined from the crack density prediction. stress_y_n1 = 34e7:1e7:65e7; stress_y_n2 = 47e7:1e7:85e7; stress_y_n4 = 44e7:1e7:76e7; stress_x_n1=stress_y_n1*ratio; stress_x_n2=stress_y_n2*ratio; stress_x_n4=stress_y_n4*ratio; % Compute the corresponding force resultant, which will be used to compute % the A matrix for different lx, ly combination N_x_n1 = stress_x_n1*t; N_x_n2 = stress_x_n2*t; N_x_n4 = stress_x_n4*t; N_y_n1 = stress_y_n1*t; N_y_n2 = stress_y_n2*t; N_y_n4 = stress_y_n4*t; % Obtain crack density values for given loads. Crack densities are fitted as cubic polynomials % based on the predicted values using Trendline function in EXCEL. % If the stress level is less than initiation stress, crack density is considered to be 0.01 % (close to ZERO but not ZERO). % And if the stress level is greater than stress limit, crack density is considered to remain PAGE 97 97 % same as the maximum value. switch n case 1 % For layup [0/90/0] for i = 1:length(N_x_n1) lambda_x(i) = 0.01; lambda_y(i) = 2.6206450E-25*stress_y_n1(i)^3 4.5320139E-16*stress_y_n1(i)^2 + ... 2.8746996E-07*stress_y_n1(i) 5.5172903E+01; end case 2 % For layup [0/90_2/0] for i = 1:length(N_x_n2) if stress_y_n2(i)<77e7 lambda_x(i) = 0.01; else lambda_x(i) = 2.1493319 E-23*stress_y_n2(i)^3 5.30 50630E-14*stress_y_n2(i)^2+... 4.3677142E-05*stress_y_n2(i) 1.1989756E+04; end lambda_y(i) = 1.1562576E-25*stress_y_n2(i)^3 2.7728675E-16*stress_y_n2(i)^2 + ... 2.4313246E-07*stress_y_n2(i) 6.5196048E+01; end case 3 % For layup [0/90_4/0] for i = 1:length(N_x_n4) lambda_x(i) = 1.0773194E-25*stress_y_n4(i)^3 2.3211953E-16*stress_y_n4(i)^2 + ... 1.7951598E-07*stress_y_n4(i) 4.3240891E+01; if stress_y_n4(i) < 55e7 lambda_y(i) = 0.01; else lambda_y(i) = 3.6124431 E-26*stress_y_n4(i)^3 1.35 45644E-16*stress_y_n4(i)^2+... 1.6900304E-07*stress_y_n4(i) 5.7924787E+01; end end otherwise display ('wrong choice') end % % Compute cross area from two basic cases % Amatrix_n1, Amatrix_n2, and Amatrix_n4 are functions to obtain [A] matrix for given (lx, ly) switch n case 1 for i = 1:length(N_x_n1) A = Amatrix_n1(lambda_x(i), lambda_y(i)); N = [N_x_n1(i); N_y_n1(i)]; epsilon = inv(A)*N; dU = dU1*epsilon(1)/epsilon_x1+dU2*epsilon(2)/epsilon_y2; dV = dV1*epsilon(1)/epsilon_x1+dV2*epsilon(2)/epsilon_y2; area(i) = 4*dU*dV; end case 2 for i = 1:length(N_x_n2) A = Amatrix_n2(lambda_x(i), lambda_y(i)); N = [N_x_n2(i); N_y_n2(i)]; epsilon = inv(A)*N; dU = dU1*epsilon(1)/epsilon_x1+dU2*epsilon(2)/epsilon_y2; PAGE 98 98 dV = dV1*epsilon(1)/epsilon_x1+dV2*epsilon(2)/epsilon_y2; area(i) = 4*dU*dV; end case 3 for i = 1:length(N_x_n4) A = Amatrix_n4(lambda_x(i), lambda_y(i)); N = [N_x_n4(i); N_y_n4(i)]; epsilon = inv(A)*N; dU = dU1*epsilon(1)/epsilon_x1+dU2*epsilon(2)/epsilon_y2; dV = dV1*epsilon(1)/epsilon_x1+dV2*epsilon(2)/epsilon_y2; area(i) = 4*dU*dV; end otherwise display('wrong choice of n'); end % Compute normalized permeability k/C using Equation 2-11. % Note that since crack densities in every plies are assumed to be same, the equation for % permeability is simplified to following form. In generally case, should following % Equation (2.11) % Also Note that the unit for lambda is 1/cm, hence there's a factor of 1e-4 appeared in % the formula. k = lambda_x.*lambda_y.*area*1E4; % Plot the computed permeability in terms of stress_y. switch n case 1 plot(stress_y_n1, k); case 2 plot(stress_y_n2, k); case 3 plot(stress_y_n4, k); end hold on; PAGE 99 99 APPENDIX D SUPERPOSITION OF STRAIN ENERGY RELEASE RATE In linear elastic fracture mechanics, the stre ss intensity factors for angled crack subjected to biaxial loading (Figure D-1) is given by (Anderson, 1995) B K K B K KI II I I 1 cos sin sin cos0 2 2 0 (D-1) where is the angle between the crack and the 1 plane, B is the ratio between 2 and 1 (2<1), KI(0) is the Mode I stress intensity factor when = 0. In 2-D problem, stress intensity factor is related to strain energy release rate by E K E K GII I 2 2 (D-2) where E E for plane stress and 21 E E for plane strain. Here we only consider Mode I and Mode II. For different loading case, the stress intensity factors for same Mode can be superposed. Then the total strain energy release rate can be ob tained from the superposed total stress intensity factors for each mode using Equation D-2. In our case, two kinds of load ing are considered: uniaxial loading and thermal loading. The strain energy release rates under mechanical lo ading and under thermal loading are denoted as GM and GT. In the uniaxial loading, B = 0, = 90 where is the fiber direction in angle ply. From Equation D-1, the stress intensity factors are cos sin sin0 2 0 I II I IK K K K (D-3) Substituting Equation D-3 back into Equation D-2 yields sin0E G KM I (D-4) and PAGE 100 100 cos sin E G K E G KM II M I (D-5) For the thermal loading case, we assume B = 1. So only Mode I stress intensity factor left, which is E G KT I (D-6) By superposing, the total stress intens ity factors for Mode I and Mode II are cos sin E G K E G E G KM Total II T M Total I (D-7) Subsequently, the total strain energy release rate is given by sin 2T M T M TotalG G G G G (D-8) These equations are derived based on linea r elastic fracture mechanics for isotropic homogeneous material. In this work, we use th ese equations for composite material as an approximation. PAGE 101 101 Figure D-1. Angled crack subj ected to biaxial loading. PAGE 102 102 LIST OF REFERENCES Anderson, T. L. (1995). Fracture Mechanics: Fundamenta ls and Applications (2nd ed.) Boca Raton: CRC Press. Aoki, T., Kumazawa, H., & Susuki, I. (2002, Oct.). Propellant Leakage through Laminated Structures. Proceeding of American Socie ty for Composites (ASC) 17th Annual Technical Conference West Lafayette, IN. ASTM D1434-82 (1992). Standard Test Method for Dete rmining Gas permeability Characteristics of Plastic Film and Sheeting ASTM, 203-213. Bapanapalli, S. K., Sankar, B. V., & Primas, R. J. (2006). Microcracking in Cross-ply Laminates due to Biaxial Mechanical and Thermal Loading. AIAA Journal, 44 (12), 2949-2957. Bechel, V. T., & Arnold, F. (2006, May). Pe rmeability of Polymer Composites for Cryogenic Applications. Proceeding of 47th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics & Materials conference Newport, RI, AIAA 2006-2015. Bechel, V. T., Camping, J. D. & Kim, R. Y. (2005). Cryogenic/elevat ed Temperature Cycling Induced Leakage Paths in PMCs. Composites: Part B, 36 (2), 171-182. Berthelot, J. M. (2003). Transv erse Cracking and Delamination in Cross-ply Glass-fiber and Carbon-fiber Reinforced Plastic Lamina tes: Static and Fatigue Loading. Applied Mechanics Review, 56 (1), 111-147. Choi, S. (2005). Micromechanics, Fracture Mechanics and Gas Permeability of Composite Laminates for Cryogenic Storage Systems Doctoral Dissertation. University of Florida, Gainesville, FL. Grenoble, R. W., & Gates, T. S. (2006, Ap ril). Hydrogen Permeability of Polymer Matrix Composites at Cryogenic Temperatures. Proceeding of 46th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics & Materials conference Austin, TX, AIAA 2005-2086. Hutchinson, J.W., & Suo, Z. (1992). Mixe d Mode Cracking in Layered Materials. Advances in Applied Mechanics, 29 63-191. Jamison, R. D., Schulte, K., Reifsnider, K. L ., & Stinchcomb, W. W. (1984). Characterization and Analysis of Damage Mechanism in Tension-Tension Fatigue of Graphite/Epoxy Laminates. In D. J. Wilkins (Ed.), Effects of Defects in Composites Materials (pp. 21-55). West Conshohocken: ASTM International. Kessler, S. S., Matuszeski, T., & McManus, H. ( 2001, Sept.). Cryocycling and Mechanical Testing of CFRP for theX-33 Liquid H2 Fuel Tank Structure. Proceedings of American Society for Composites (ASC) 16th Technical Conference Blacksburg, VA, ASC-2001-093. Kulichenko, A.V. (2005). Theoretical Analysis Calculation, and Pr ediction of the Air Permeability of Textiles. Fibre Chemistry, 37 (5), 371-380. PAGE 103 103 Kumazawa, H., Aoki, T., & Susuki, I. (2003). An alysis and Experiment of Gas Leakage through Composite Laminates for Propellant Tanks. AIAA Journal, 41 (10), 2037-2044. Kumazawa, H., Susuki, I., & Aoki, T. (2006). Gas Leakage Evaluation of CFRP Cross-ply Laminates under Biaxial Loadings. Journal of Composite Materials, 40 (10), 853-871. Kumazawa, H., Hayashi, H., Susuki, I., & Uts unomiya, T. (2006). Damage and Permeability Evolution in CFRP Cross-ply Laminates. Composite Structures, 76 73-81. Lavoie, J. A., & Adolfsson, E. (2001). Stitch Cr acks in Constraint Plies Adjacent to a Cracked Ply. Journal of Composite Materials, 35 (23), 2077-2097. McManus, H. L., Faust, A., & Uebelhart, S. (2001, Sept.). Gas Permeability of Thermally Cycled Graphite-epoxy Composites. Proceeding of American Socie ty for Composites (ASC) 16th Technical Conference Blacksburg, VA, ASC-2001-092. Nairn, J. A. (2000). Matrix Microcracking in Co mposites. In R. Talreja & J.A. Manson (Eds.), Polymer Matrix Composites, Chapter 13 (pp. 403-432), Volume 2 of Kelly, A and Zweben (Eds.), Comprehensive Composite Materials Elsevier Science. Nairn, J. A., Hu, S., & Bark, J. S. (1993). A critical Evaluation of Theories for Predicting Microcracking in Composite Laminates. Journal of Materials Science, 28 (18), 5099-5111. Nettles, A. T. (2001). Permeability Testing of Impacted Composite Laminates for Use on Reusable Launch Vehicles. NASA/TM-2001-210799 Noh, J., Whitcomb, J., Peddiraju, P., & Lagoudas, D. (2004, April). Prediction of Leakage Rate through Damage Network in Cryogenic Composite Laminates. Proceeding of 45th AIAA/ASME/ASCE/AHS/ASC Structures, Stru ctural Dynamics & Materials Conference Palm Springs, CA, AIAA-2004-1861. Peddiraju, P., Grenoble, R., Fried, K., Gates, T., & Lagoudas, D. (2005, April). Analytical Predictions and Experimental Measurements of Hydrogen Permeability in a Microcrack Damaged Composite. Proceeding of 46th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics & Materials Conference Austin TX, AIAA-2005-2087. Peddiraju, P., Lagoudas, D., Noh, J., & Whitcom b, J. (2004, April). Numerical Modeling of Cryogen Leakage through Composite Laminates. Proceeding of 45th AIAA/ASME/ASCE/AHS/ASC Structures, Stru ctural Dynamics & Materials Conference Palm Springs, CA, AIAA-2004-1862. Peddiraju, P., Noh, J., Whitcomb, J., & Lagoudas, D.C. (2007). Prediction of Cryogen Leak Rate through Damaged Composite Laminates. Journal of Composite Materials, 41 (1), 41-71. Roy, S., & Benjamin, M. (2004). Modeling of Permeation and Damage in Graphite-epoxy Laminates for Cryogenic Fuel Storage. Composites Science and Technology, 64 (13-14), 2051-2065. PAGE 104 104 Roy, S., & Nair A. (2006). Modeling of Permeability in Composite Laminates with Delaminations and Stitch Cracks. Proceeding of 47th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics & Materials Conference Newport, RI, AIAA-2006-2094. Schulz, W. A. (2005). Determination of Residual Stress and Thermal Behavior for Composite Laminates Master Thesis, University of Florida, Gainesville, FL. Speriatu, L. M. (2005). Temperature Dependent Mechanical Pr operties of Composite Materials and Uncertainties in Experimental Measurements Doctoral Dissertation, University of Florida, Gainesville, FL. Stokes, E. (2002, Sept.). Hydrogen Permeability of Polymer Based Com posite Tank Material under Tetra-Axial Strain. Proceedings of the Fifth Conference on Aerospace Materials, Processes and Environment Technology (AMPET) Huntsville, Alabama. Sun, C.T., & Jih, C. J. (1987). On Strain Energy Release Rate for Interfacial Cracks in Bimaterial media. Engineering Fracture Mechanics, 28 (1), 13-20. Xu, J. L., & Sankar, B. V. (2007). Parametric Investigation of Gas Pe rmeability in Cross-Ply Laminates Using Finite Elements. AIAA Journal, 45 (4), 934-941. Yokozeki, T., Aoki, T., Ogasawara, T., & Ishi kawa, T. (2005). Effects of layup angle and ply thickness on matrix crack interaction in contiguous plies of composite laminates. Composites: Part A, 36 (9), 1229-1235. Yokozeki, T., & Ishikawa, T. (2005). ThroughThickness Connection of Matrix Cracks in Laminate Composites for Propellant Tank. Journal of Spacecraft and Rockets, 42 (4), 647-653. Yokozeki, T., Aoki, T., & Ishikawa, T. (2005). Consecutive Matrix Cracking in Contiguous Plies of Composite Laminates. International Journal of Solids and Structures, 42 (9-10), 2785-2802. PAGE 105 105 BIOGRAPHICAL SKETCH Jianlong Xu was born in Gansu, China in 1978. He received his Bachelor of Science degree in mechanical engineer ing from the University of Science and Technology of China (USTC) in July 2000. At the same time, he obtai ned a secondary B.S. degree in computer science from USTC. He continued his study in USTC pursu ing his masters degree. In July 2003, he received his Master of Science degree in so lid mechanics. In August 2003, Jianlong Xu was admitted to the Department of Mechanical and Aerospace Engineering at University of Florida to pursue his doctoral degree. |