A ROBUST CRASH SIMULATION MODEL FOR COMPOSITE STRUCTURES By Danghe Shi               A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of  Mechanical Engineering - Doctor of Philosophy 2016    ABSTRACT A ROBUST CRASH SIMULATION MODEL FOR COMPOSITE STRUCTURES By Danghe Shi Fiber reinforced composites are widely used in aerospace, automotive and other industries due to their high stiffness & strength-to-weight ratio, corrosion resistant and energy absorption ability. The use of composites in primary energy absorbing components in vehicles, however, is very limited. The lack of robust and accurate models for the prediction of crashworthiness performance of composite structures is a critical factor. The capability of crash simulations is often examined with axial impact of tubes, a benchmark problem. Among the large amount of published works on modeling composite tubes under axial impact, very few were able to predict both the failure morphologies and force-displacement responses of a crushing tube, especially for the ones without a plug initiator at the crash front. One of the problems is related to the limitations of material constitutive models. The existing composite material models used in crash simulations generally contain a few non-measurable parameters while ignoring the irreversible strains which can be important to the prediction of energy absorption. The other problem is related to the robustness of the finite element (FE) models. Traditional elements tend to have instability issues under in-plane compression.  In order to solve these problems, this study proposes a new crash model which composes of an Enhanced Continuum Damage Mechanics (ECDM) model and a Shell-Beam (SB) element modeling method. The ECDM model consists of a pre-failure sub-model, based on   the modified Ladevèze model, and a post-failure sub-model. It considers both the matrix plasticity and irreversible strains due to damage. The SB element method enhances the shell element with real out-of-plane properties by introducing beam elements in the out-of-plane direction. As a result, the SB element is stable under in-plane compression while retaining the efficiency of the shells. The ECDM has been implemented in LS-DYNA as a user-defined material model.  To evaluate the predictive capability of this new crash model, simulations are carried out for both quasi-static tests and dynamic tube crash tests of 2D triaxial braided composites (2D3A). For dynamic tube crash tests, two different groups of tests are done on composite tubes, with and without a plug initiator. The results show that the ECDM model is able to simulate the force and energy responses as well as the failure morphologies more accurately than a widely used CDM model available in LS-DYNA. On the other hand, the SB element method allows one to capture most of the damage modes of a composite tube under axial crash without the numerical difficulties experienced by using traditional elements, like negative element volume, infinite small timesteps, and instability. Both ECDM model and SB element method are also robust enough to be used on different geometries.  To have a better understanding of ECDM model and shell-beam method, a sensitivity study is carried out at the end of this study for key model parameters, such as the tiebreak contact strength, friction coefficient, element deletion strain, etc.     iv                      To my grandmothers and my grandfather    v ACKNOWLEDGMENTS  The PhD study is a long fruitful journey, full of frustration, hope, happiness and surprise. I am so glad I made the choice to purse a PhDs degree four years ago since I just had one of the best four years of my life. I learned how to face and deal with problems and enjoyed every part of the process. I would like to take this opportunity to acknowledge all the people involved in this wonderful journey.  I want to express my deepest gratitude to my supervisor, Dr. Xinran Xiao. I also got my Masters degree under her supervision. She guided me onto the road of research, shared her life experience with me and instructed me to be an independent thinker. Without her guidance, encouragement, understanding and patience, it would have been impossible for me to finish this dissertation. I appreciate everything she did for me.   Many thanks also go to Dr. Ronald Averill, Dr. Rigoberto Burgueño and Dr. Dahsin Liu. Thank them all for being willing to join in my dissertation committee. I appreciate their insights and helpful suggestions.  I also want to thank DOE (DE-EE0005397), Michigan State University (college of engineering, graduate school, OISS) and MGA Research Corporation for supporting this work. It is a big encouragement for me to have their faith in me.    vi I am also grateful for those people who kindly helped me throughout this research. Firstly, I would like to thank Christopher Cater for his willingness to discuss on problems and share his valuable knowledge with me all the time. I puzzled over many problems with him. Secondly, I want to thank my other fellow group members, Miao Wang, Shutian Yan, Yalla Abushawashi and Vahid Khademi. Even though we work on different projects, we support each other with positive feedbacks and suggestions. I enjoy working, studying and brainstorming with all of them every day. Thirdly, I want to express my appreciations to my friends, Fang Hou, Xing Xing, Wei Zhang, Huiyi Zhou, Mengtian Jiang, Wu Zhou, Shupeng Zhang and my girlfriend Yanjie Wang. Thank them for participating in my life. Without them, the past four years cannot be so colorful and memorable.   Last but not least, I would like to thank my grandmothers, my grandfather, my dad, my mom and my sister. Their endless love and support are always the most powerful driving force of my life. I love them all so much. I hope I have lived up to their expectations.      vii TABLE OF CONTENTS LIST OF TABLES ................................................................................................................. x LIST OF FIGURES .............................................................................................................. xi CHAPTER 1:  INTRODUCTION ............................................................................................ 1 1.1. Fiber Reinforced Composites...................................................................................... 1 1.2. Crashworthiness of Composites .................................................................................. 3 1.2.1. Experimental Work on Crashworthiness of Composites ................................ 5 1.2.2. Modeling Work on Crashworthiness of Composites ...................................... 9 1.2.2.1. Composite constitutive material models ........................................... 10 1.2.2.2. FE model ........................................................................................... 12 1.2.3. Summary ....................................................................................................... 14 1.3. Research Objectives .................................................................................................. 15 1.4. Thesis Structure......................................................................................................... 17 CHAPTER 2:  COMPOSITE MATERIAL MODELS ............................................................. 19 2.1. Considerations of a Composite Material Model for Crash Simulation ..................... 19 2.2. Literature Review on CDM Models.......................................................................... 22 2.2.1. Concept and Characters of CDM Models ..................................................... 22 2.2.2. Popular CDM Models for Composites.......................................................... 23 2.3. Literature Review on Post-Failure Models ............................................................... 28 2.3.1. Experimental work on strain-softening ......................................................... 29 2.3.2. Modeling work on strain-softening ............................................................... 30 CHAPTER 3:  ECDM MODEL .......................................................................................... 36 3.1. Proposed ECDM Material Model ............................................................................. 36 3.2. Pre-failure Material Model ........................................................................................ 37 3.2.1. Damage kinematics of the elementary ply .................................................... 38 3.2.2. Plasticity modeling and damage-plasticity coupling .................................... 41 3.3. Initial Failure Criteria ................................................................................................ 43 3.4. Strain Softening Laws ............................................................................................... 44 3.5. Post-Failure Residual Model ..................................................................................... 44 3.5.1. Post-failure residual state model ................................................................... 44 3.5.2. Post-failure irreversible strains ..................................................................... 46 3.6. Final Failure Criteria ................................................................................................. 48 3.7. Implementation of ECDM Model in LS-DYNA ...................................................... 49 CHAPTER 4: VALIDATION OF ECDM MODEL ................................................................ 51 4.1. Material and Material Properties ............................................................................... 51 4.1.1. 2D3A for quasi-static tests ............................................................................ 51 4.1.2. 2D3A for dynamic tests ................................................................................ 52 4.2. Validation of ECDM Model on Quasi-Static One Element Model .......................... 56 4.2.1. Growth of damage and plasticity .................................................................. 57 4.2.2. Tensile responses of 2D3A ........................................................................... 57   viii 4.3. Validation of ECDM Model on Dynamic Crash Model with Plug ........................... 61 4.3.1. Simulation model for composite tubes with plug initiator ............................ 61 4.3.1.1. Dynamic crash model........................................................................ 61 4.3.1.2. Delamination model .......................................................................... 62 4.3.1.3. Boundary conditions and contact definitions .................................... 63 4.3.2. Crash simulation results of thin-walled composite tubes with plug initiator 65 4.3.2.1. Compare different ECDM models .................................................... 65 4.3.2.2. Compare ECDM with MAT58 ......................................................... 73 CHAPTER 5:  SHELL-BEAM ELEMENT METHOD ............................................................ 80 5.1. Simulation Difficulties of Composite Tubes without Plug Initiator ......................... 80 5.2. Investigation of Different Element Formulations in LS-DYNA .............................. 81 5.2.1. Different element formulations in LS-DYNA .............................................. 82 5.2.2. Three test cases ............................................................................................. 85 5.2.2.1. Test I: One layer of triaxial braided composites without trigger ...... 86 5.2.2.2. Test II: One layer of triaxial braided composites with trigger .......... 88 5.2.2.3. Test III: Two layers of triaxial braided composites with trigger ...... 91 5.2.2.4. Test case summary ............................................................................ 93 5.3. Proposed New Modeling Method ............................................................................. 94 5.3.1. Shell-Beam element method ......................................................................... 94 5.3.2. Evaluation of shell-beam element method .................................................... 96 5.3.2.1. One element tests .............................................................................. 98 5.3.2.2. Cantilever Beam Tests .................................................................... 100 5.3.2.3. Plate Crash Tests ............................................................................. 103 CHAPTER 6: CRASH SIMULATION OF COMPOSITE TUBES WITH ECDM AND SHELL-BEAM ELEMENT METHOD ................................................................................... 105 6.1. Crash Simulation Models ........................................................................................ 105 6.2. Crash Simulation of 2 Tubes and 4 Tubes .................... 108 6.2.1. Damage/failure morphology ....................................................................... 108 6.2.2. Force-displacement and SEA responses ..................................................... 113 6.3. Crash Simulation of Some Other Triaxial Braided Composite Tubes without Plug Initiator .................................................................................................................. 118 CHAPTER 7: INVESTIGATION ON THE SIMULATION MODELS AND DETAILED RESULTS ............................................................................................................... 123 7.1. Sensitivity Studies ................................................................................................... 123 7.1.1. Trigger ......................................................................................................... 123 7.1.2. Eratios ......................................................................................................... 128 7.1.3. Tiebreak Contact Strength .......................................................................... 131 7.1.4 Friction coefficient ....................................................................................... 135 7.1.5 Mechanical Properties and Residual Strength ............................................. 137 7.1.6 Element deletion strains ............................................................................... 142 7.2. Investigation on detailed simulation results ............................................................ 146 7.2.1. Damage conditions ...................................................................................... 146 7.2.2. Residual stiffness ........................................................................................ 148 7.2.3. Delamination ............................................................................................... 149   ix 7.2.4. Energy breakdown ...................................................................................... 151 CHAPTER 8:  CONCLUSION AND FUTURE WORK .......................................................... 153 8.1. Conclusion .............................................................................................................. 153 8.2. Future Work ............................................................................................................ 155 APPENDICES .............................................................................................................. 160 APPENDIX A: ECDM HISTORY VARIABLES ............................................................. 161 APPENDIX B: MATERIAL CARDS FOR ECDM MODEL IN LS-DYNA ...................... 166 BIBLIOGRAPHY ........................................................................................................ 168     x LIST OF TABLES Table 2-1. Comparison of different CDM models. .......................................................... 35 Table 4-1. Input data of 2D3A for quasi-static simulation. ............................................. 52 Table 4-2. Dimensions of different composite tubes ....................................................... 53 Table 4-3. Input data of 2D3A (GM) for dynamic simulation ........................................ 55 Table 5-1. Comparison between thin-thick shell and thick-thick shell elements. ........... 85 Table 5-2. Material properties of beam elements ............................................................ 97 Table 5-3. Comparison of input moduli and simulated results by one shell-beam element .................................................................................................................. 100 Table 5-4. Summary of specimen configurations. ......................................................... 101 Table 6-1. Comparison of computational times of 2 tube models with plug initiators using the shell element and the shell-beam element. ..................... 107   xi LIST OF FIGURES Figure 1-1. Classification of fabric forms [1]. ................................................................... 2 Figure 1-2. Typical SEA values of different materials. [9] ............................................... 4 Figure 1-3.  Application of fiber-reinforced composite tubes in the front of a vehicle [10]. ........................................................................................................................... 4 Figure 1-4. Schematic of axial crushing with and without a trigger mechanism [18]. ...... 5 Figure 1-5. Various variables that may influence the energy absorption abilities of composite materials [9]. ............................................................................................ 8 Figure 1-6. Simulation of a [0/±30]2-braided composite tube without plug initiator by (a) one layer of shell elements [47] and (b) two layers of shell elements [11]. ...... 12 Figure 1-7. (a) Progressive collapse of a CFRP quarter-model in [26]; (b) use of fixed debris wedge and rigid walls to assist the formation of failure morphologies in [22]; (c) final deformation of the simulation model obtained in [28]. .................... 14 Figure 1-8.  Schematic of damage/failure mechanisms of braided composite tubes in the crash front with and without plug initiator [22]. ............................................... 16 Figure 2-1. Comparison of the stress-strain response predicted by MAT 58 to experiment data of a [0/±30]2 triaxial braided composite. The specimen was loaded in compression, then unloaded and loaded to failure in tension [11]. ......... 21 Figure 2-2. Effect of variations in the exponent, m, on the longitudinal stress-strain behavior predicted by the MLT model. [118] ......................................................... 33 Figure 3-1. Schematic of the uniaxial stress-strain response of the ECDM model. The model is composed of two sub-models: a pre-failure model and post-failure model. The post-failure model considers two stages:  strain softening and residual state with two different property evolution laws. ...................................... 37 Figure 3-2. Elementary ply [82] ....................................................................................... 38 Figure 3-3. Under compression, a) kink bands formed in a unidirectional carbon reinforced composite [104]; b) kink bands formed in a 2D3A material [105]; c) micro-buckling observed in a 2D3A material [103] ............................................... 40 Figure 3-4. Schematic of the phenomenon of increasing residual load carrying capability with deformation. (a) Assume 90% of the fibers fail earlier than the remaining 10% fibers; (b) the response of the material is the summation of the two responses in (a). ............................................................................................... 45   xii Figure 3-5. Energy comparison between (a) CDM model and plastic CDM model; (b) predicted results by MAT58 in LS-DYNA and LLNL compression tension test data [11]. ................................................................................................................. 48 Figure 3-6. The general working flow chart of the ECDM model implemented in LS-DYNA. .................................................................................................................... 50 Figure 4-1. a) Schematic of the preparation of uncut and cut specimens; b) comparison of results for the cut and uncut [±25]s braided specimens. [101] ....... 54 Figure 4-2. Shear damage and plasticity responses and fitted curves. ............................. 56 Figure 4-3. Comparison of simulation results of ECDM model and MAT58 and experimental results on  (a) 0º, (b) 45º and (c) 90º specimens. .............................. 59 Figure 4-4. Comparison of predicted damage evolution in longitudinal direction by MAT58 and ECDM model. .................................................................................... 60 Figure 4-5. Detailed FE crash model and the drop tower test setup at UBC [22]. .......... 62 Figure 4-6. Experiment setup designed to measure the friction coefficient between the plug initiator and the inner surface of a composite tube ......................................... 64 Figure 4-7. Comparison of the crash front morphologies predicted by different ECDM models with the one obtained by experiment. ........................................................ 66 Figure 4-8. Predicted failure morphology is highly affected by residual strains. (a) Prediction with ECDM-ir; (b) Prediction with ECDM-no-ir and ECDM-pre-ir. ... 67 Figure 4-9. Longitudinal strain responses of two adjacent elements from the outer surface and inner surface for cases predicted with ECDM-ir and ECDM-pre-ir; positive and negative strains corresponding to tension and compression load. ...... 68 Figure 4-10. a) Predicted force-displacement and b) predicted velocity responses by different ECDM models. ......................................................................................... 69 Figure 4-11. Detailed comparison of a) peak force, b) average plateau force, c) load ratio, d) crash length, e) SEA and f) average plateau force vs. crash length. ......... 70 Figure 4-12. Comparison of the predicted internal energy densities by ECDMs ............ 72 Figure 4-13. Predicted longitudinal strains and internal energy densities of the same element in the inner layer of a tube model by three ECDM models....................... 72 Figure 4-14. Comparison of the crash front morphologies predicted by MAT58 and ECDM models with the one obtained by experiment. ............................................ 74 Figure 4-15. Comparison of experimental results with the predicted force-displacement responses using the ECDM model and MAT58. .............................. 75   xiii Figure 4-16. Detailed comparison of a) peak force, b) average plateau force, c) load ratio, d) crash length and e) SEA. ........................................................................... 75 Figure 4-17. Comparison of the predicted internal energy densities by MAT58 and ECDM model. ......................................................................................................... 77 Figure 4-18. Predicted longitudinal strains and internal energy densities of the same elements in the inner layer of a tube model by MAT58 and ECDM. ..................... 78 Figure 5-1. Unstable shell element model could predict different failure morphologies under the same loading conditions. ......................................................................... 83 Figure 5-2. Element discretization and nodal degree of freedom in local coordinate system for (a) fully integrated shell element; (b) thin-thick shell element; (c) thick-thick shell element. [139] .............................................................................. 84 Figure 5-3. Comparison of the simulation results of 1 layer of triaxial braided composite tubes without trigger by (a) fully integrated shell elements; (b) thin-thick shell element with ELFORM 2; (c) thick-thick shell element with ELFORM 5. ............................................................................................................ 87 Figure 5-4. Progressive crash of a l-layer composite tube made by fully integrated shell elements without trigger in the crash front. .................................................... 87 Figure 5-5. Schematic of the triggers used in the composite tubes in test case II modeled by (a) fully integrated shell elements; (b) thin-thick shell element with ELFORM 2; (c) thick-thick shell element with ELFORM 5. ................................. 89 Figure 5-6. Comparison of the simulation results of 1 layer of triaxial braided composite tubes with trigger by (a) fully integrated shell elements; (b) thin-thick shell element with ELFORM 2; (c) thick-thick shell element with ELFORM 5.... 89 Figure 5-7. Progressive crash of a l-layer composite tube made by (a) thin-thick shell elements and (b) shell elements with triggers in the crash front. ............................ 90 Figure 5-8. Schematic of the triggers used in the composite tubes in test case III modeled by (a) fully integrated shell elements; (b) thin-thick shell element with ELFORM 2; (c) thick-thick shell element with ELFORM 5. ................................. 92 Figure 5-9. Comparison of the simulation results of 2 layers of triaxial braided composite tubes with trigger by (a) fully integrated shell elements; (b) thin-thick shell element with ELFORM 2; (c) thick-thick shell element with ELFORM 5.... 92 Figure 5-10. Schematic of the proposed shell-beam element .......................................... 95 Figure 5-11. Stress-strain response of beam elements ..................................................... 98   xiv Figure 5-12. One element tests on single shell-beam element. T, C, S refer to tension, compression and shear separately. .......................................................................... 99 Figure 5-13. Cantilever beam finite element models with a length to depth ratio 8:1. . 101 Figure 5-14. Force-deflection responses predicted by different types of elements under three length-to-depth ratios. The results are filtered by SAE1000 in LS-PrePost. ................................................................................................................. 102 Figure 5-15. Plate crash tests by using three types of elements..................................... 104 Figure 6-1. Detailed FE crash model without plug initiator and the drop tower test set-up. [22] .................................................................................................................. 106 Figure 6-2. Schematic of the cross- 2 tube model ................. 106 Figure 6-3. Comparison of experimentally observed and predicted damage/failure morphologies of (a) 2 and (b) 4 tubes with plug initiator. ................................................................................................................. 108 Figure 6-4. Progressive evolution of predicted damage/failure morphologies in the crash front of 4 tube with plug initiator from cross-section views. . 109 Figure 6-5. Comparison of experimentally observed and predicted damage/failure morphologies of (a) 2 tube and (b) 4 tube without plug initiator. ......................................................................................................... 111 Figure 6-6. Progressive evolution of predicted damage/failure morphologies in the crash front of (a) 2 tube and (b) 4 tube without plug initiator from cross-section views. ........................................................................ 112 Figure 6-7. Comparison of force-displacement responses for both 2-ply and 4-ply  tubes crushed with and without plug initiators. ........................................... 113 Figure 6-8. Detailed comparison between predictions and experiment results in terms of a) peak force, b) average plateau force, c) load ratio, d) crush length and e) SEA. ...................................................................................................................... 116 Figure 6-9. Comparison of damage/failure morphologies of (a) 2, (b) round 2 and (c) 2 tubes without plug initiator ................. 119 Figure 6-10. Comparison of force-displacement responses for (a) 2, (b) round 2 and (c) 2 tubes without plug initiator. ................ 120 Figure 7-1. Schematic of 5 different triggers and comparison of the predicted failure morphologies of 2×22 tube models under conditions both with and without a plug initiator. ......................................................................................... 125   xv Figure 7-2. Force-displacement responses of 2 tube models with different triggers (a) with plug initiator; (b) without plug initiator. ..................... 127 Figure 7-3. Different combinations of eratios result in different crash front failure morphologies......................................................................................................... 129 Figure 7-4. Comparison of different failure morphologies of 2 tube models with different eratio combinations. 0.5-0.75-0.5 corresponds to the case with eratio1=0.5, eratio2=0.75 and eratio12=0.5, and so on. ............................... 130 Figure 7-5. Force-displacement responses of 2 tube models with different eratio combinations. ............................................................................... 131 Figure 7-6. Comparison of different failure morphologies of simulation tubes with different tiebreak contact strength. 0.015-0.025 corresponds to the case with NFLS=0.015, SFLS=0.025, and so on. ................................................................. 133 Figure 7-7. Force displacement responses of 2 tube models with different tiebreak contact strength. ....................................................................... 133 Figure 7-8. Crash front morphologies and delamination conditions under different combinations of tiebreak strength parameters. In LS-DYNA, the contact gap fringe gives delamination fringe valued from 0-1. 0 refers to intact area while 1 refers to fully delaminated area. ............................................................................ 134 Figure 7-9. Comparison of different failure morphologies of 2 tube models with different friction coefficient. ............................................................ 136 Figure 7-10. Force displacement responses of 2 tube models with different friction coefficient. ................................................................................. 136 Figure 7-11. Transverse stress-strain relationships with (a) varied mechanical properties (MP), (b) varied residual strengths. ..................................................... 138 Figure 7-12. Comparison of different failure morphologies of 2 tube models with different mechanical properties. ....................................................... 139 Figure 7-13. Force displacement responses of 2 tube models with different mechanical properties............................................................................. 140 Figure 7-14. Comparison of different failure morphologies of 2 tube models with different residual strength................................................................. 141 Figure 7-15. Force displacement responses of 2 tube models with different residual strength. .................................................................................... 141   xvi Figure 7-16. Comparison of different failure morphologies of  tube models with different element deletion strains. .................................................... 143 Figure 7-17. Force displacement responses of 2 tube models with different element deletion strains. ......................................................................... 145 Figure 7-18. Comparison of three major damage states of 2 tube models both with and without plug initiator during crash. ................................................ 148 Figure 7-19. Comparison of the elastic moduli of 2 tube models both with and without plug initiator during crash. All the plots are normalized by their initial elastic moduli. .................................................................................... 149 Figure 7-20. Comparison of the delamination states of 2 tube models both with and without plug initiator...................................................................... 150 Figure 7-21. Energy breakdown of 2 tube models (a) with and (b) without plug initiator............................................................................................. 152 Figure 8-1. Schematic of the amount of bias tows griped in (a) a straight-sided specimen and (b) a bowtie-shaped specimen [1] .................................................. 157 Figure 8-2. Working schematic of EDGE function. ....................................................... 159     1 CHAPTER 1:  INTRODUCTION   Fiber reinforced composites are widely used in aerospace and automotive industries due to their high stiffness & strength-to-weight ratio, corrosion resistant and energy absorption ability. However, the progress on the crashworthiness prediction of fiber reinforced composite structures has been relatively slow which to some extent limited the application of composites for massive productions [1]. This research is dedicated to improve the crashworthiness prediction of fiber reinforce composite structures.  This chapter firstly briefly introduces the fiber reinforced composites, then reviews studies on crashworthiness of composites and proposes the research objectives and finally introduces the thesis structure.  1.1. Fiber Reinforced Composites  Fiber reinforced composites refer to the materials in which a polymer matrix component is reinforced by a stronger fiber component. For simplicity they will also be called composites or composite materials in this study. Fiber reinforced composites can be classified by the fiber reinforcements which are available in a wide range of size and forms. As shown in Figure 1-1, there are two categories of fibers: discontinuous and continuous fibers. The forms of discontinuous fibers have short and long fibers, aligned and random chopped strand mats. The forms of continuous fibers have continuous fiber mats, unidirectional or multi-axial laminates, woven, knitted, and braided fabrics. Among   2 them, unidirectional/multiaxial laminates, woven and braided fabrics are most widely used on crash structures.    Figure 1-1. Classification of fabric forms [1].   The mechanical properties of these composites are affected by the number of fiber bundles, or tows, the number of fibers per tow, and the weave pattern. The reinforcements in the form of biaxial, triaxial, and multiaxial fabrics exhibit up to 30 percent higher tensile strength than woven fabrics which accumulate stresses due to fiber waviness (undulation). Multiaxial reinforcements can be designed to meet specific requirements and perform multiple tasks such as providing good surface finish and structural integrity besides impact resistance [1]. For example, triaxial braided composites can offer an isotropic design by utilizing axial and angled fiber bundles in a single plane. In general, braided composites can offer better damage resistance, torsional stability, and bearing strength compared to unidirectional or weaved composites [1,2]. This research is focused on 2D triaxial braided (2D3A) composites. However, it should be noted that the material model and the modeling technique developed in this work are applicable to other composites as well.     3 1.2. Crashworthiness of Composites  Crashworthiness is the ability of a structure to protect its occupants during crash through controlled failure modes that enable the maintenance of a gradual decay in the load profile during energy absorption [3,4].  The energy absorbing ability of a material is measured by its specific energy absorption (SEA) value. As shown in Figure 1-2, composites offer higher SEA over the traditional materials. The automotive and aircraft industries have attempted to use composites in a wide range of crashworthiness applications [5,6]. Belingardi and Chiandussi [7] presented four applications of composite material to structural parts of a vehicle like the front longitudinal beams. It has been shown experimentally that when properly designed, composite materials can absorb more energy per unit mass of material specific energy than the conventional metals [8,9]. Park et al. [1] investigated the opportunities of using advanced plastics and composites in vehicle lightweighting. It was found that by using 2D3A composites in light weight trucks such as Silverado, the vehicle weight could be reduced efficiently while maintaining equivalent crash and structural performance in the frontal impact condition.    4  Figure 1-2. Typical SEA values of different materials. [9]  An axial crush simulation of composite tubes is often used as a benchmark to evaluate the predictive capability of composite crash simulation models [11]. The tubes are representative structures for energy management system in vehicles, such as the front rail, as shown in Figure 1-3 [10].    Figure 1-3.  Application of fiber-reinforced composite tubes in the front of a vehicle [10].  A tube structure under axial crushing force may fail catastrophically or progressively, as depicted in Figure 1-4 [18]. Progressive failure at the crush front is preferred since it is stable and able to absorb more energy, as indicated by the area under the load-displacement curve. Catastrophic failure is generally caused by global buckling or rupture   5 of the tube away from the crush front. This may be avoided by incorporating a trigger mechanism at the front end of the tube. The damage/failure process of a crushing tube with a trigger is as follows: damage and breakage of the trigger, delamination, separation and splaying of the inner and outer layers of composite walls and formation of debris wedges in the corners. The whole process is accompanied by contact friction.   Figure 1-4. Schematic of axial crushing with and without a trigger mechanism [18].  A large amount of work have been done on such composite tubes both experimentally and numerically in order to have a better understanding on the damage and failure mechanisms of composites. This section reviews some of the most cited experimental and numerical works in this field separately.  1.2.1. Experimental Work on Crashworthiness of Composites    6 In the early days, the experiments on composites were focused on the comparison of materials and geometries. Thornton and Edwards (1982) [12] examined the energy absorbing capacity of a wide range of structures with variation in both composite materials and geometry. It was shown that the fracture collapse depending on the tube geometry. Cylindrical wall sections were found to be more stable than planar wall sections. Below a certain thickness to section size ratio, the collapse of a structure tended to become unstable.  Glass and graphite fiber reinforced plastic (FRP) composites were more reliable than Kevlar since Kevlar or hybrids containing Kevlar FRP tubes showed unstable collapse by buckling. Farley (1983) [8] found that chamfered and notched tubes resulted in reduced initial peak loads without affecting the sustained crushing loads. The static and dynamic test produced similar failure modes and post crushing integrity.  Gradually, more researches started to concentrate on the failure mechanisms. Hull (1991) [13] demonstrated the two basic failure mechanisms were fragmentation and splaying. He firstly found a close relationship between the load-bearing capacity of a tube undergoing progressive crushing and the micromechanisms associated with crush zone morphology by examining (0/90/90/0) glass fiber-polyester resin tube and woven  glass cloth-epoxy resin tube, which typified many of the basic features observed in a wide range of materials. Hull also studied the micromechanisms of crush by changing the distribution of the fibers. It was concluded that an increase in the proportion of hoop fibers reduce the splaying of the axially oriented fibers and leading to fragmentation increases. Following this work, Starbuck et al. (2000) [14] studied the energy absorption mechanisms by isolating the damage modes associated with the frond formation (splaying mode) using a   7 new designed test fixture. Some damage mechanisms were activated by adjusting a profile constraint. This fixture was also used by Lee et al. (2002) [15] and Jacob et al. (2006) [16].  Meanwhile, more researches have been done on the effects of some detailed factors, like temperature, trigger shape, crushing speed, microstructure and etc. Mamalis et al. (1997) [5] reviewed a large amount of work on the crashworthy capability of composite material structures. It was pointed out: 1) In general, specific energy absorption of various composite materials decreases with increasing temperature above about 0°C; 2) Specimen geometry has a strong effect on the energy absorbing capability  of composite shells and corners have a negative influence on it; 3) The energy absorbed by tulip triggered specimens show significantly higher crashworthy behavior than the bevel triggered ones of same geometry and materials; 4) Crushing speed affects the energy absorption capability of axially loaded shells but the change of the  specific energy depends on the material properties, etc. Ramakrishna (1997) [9] reviewed the effect of microstructural variables, like fiber/matrix material properties, fiber architecture on the energy absorption characteristics of composite material. It was pointed out that the crushing characteristics of a material could generally be described in terms of: 1) crushing mode & crush zone morphology; 2) specific energy absorption capability. Lau et al. (2012) [17] also reviewed the effect of geometrical designs on axial crushing of composites. The structure shape, geometry and triggering design could largely affect the peak load and specific energy absorption. Various variables that may influence the energy absorption abilities of composite materials are listed in Figure 1-5.   8  Figure 1-5. Various variables that may influence the energy absorption abilities of composite materials [9].  In some later work, more details properties were characterized. Mamalis et al. (2004) [19] have done a series of experiments on woven carbon fiber reinforced composites tubes with square cross-section to investigate the influence of some geometric features, e.g. tube axial length, aspect ratio and wall thickness on the static crash response. Without trigger mechanisms, three distinct collapse modes were observed, i.e. progressive end-crushing, unstable local tube wall buckling and mid-length collapse. Only the progressive end-crushing resulted in high SEA, and the main failure mechanisms associated with   9 energy absorption were the frictional effects in the crush zone and the contact areas. Brimhall (2005) [20] studied the friction energy absorption in fiber reinforced composites by designing a roller fixture to quantify the energy change under quasi-static and dynamic conditions. It was found that the sliding friction was responsible for the decrease in overall SEA from quasi-static compression to dynamic crush. Feraboli et al. (2009) [21] designed a test method to characterize the SEA of fabric carbon/epoxy square tubes, C-channels and corner elements of varying dimensions. The fiber tensile fracture and tearing at the corners of a square tube was found to be responsible for more energy absorption than frond formation and splaying of the large flat segments.  In summary, the composite tubes deform and fail in a manner different from similar structural components made of conventional materials like steel and aluminum because of different micro-failure modes include fiber breakage, matrix cracking, fiber matrix debonding, delamination, etc. As will be shown in 1.2.2, these complex failure mechanisms make it even more difficult to model the collapse behavior of fiber reinforced composites.  1.2.2. Modeling Work on Crashworthiness of Composites  Composite material models for crash simulation have been investigated in a number of studies. The major modeling works on crashworthiness of composites started in 1990s by using finite element method. As is well known, laminated composites exhibit a highly hierarchical structure [50] and the modeling techniques can be classified by different length scales: micro-scale, meso-scale and macro-scale [51]. The micro-scale approaches   10 and multi-scale approaches are very useful in terms of understanding some underlying damage and failure mechanisms but too computationally expensive and far from complete to be used on crash simulation of structures. It is much reasonable to build a direct relationship between experiment and simulation models in macro-scale since this scale provided enough details on the damage/failure mechanisms and the possibilities to capture every necessary parameter used in simulation by experiment. Maimi et al. [52] also pointed out that if a structure exhibits stable crack propagation, macro-scale models are more appropriate to predict the structural collapse.  This section reviews the macro-scale modeling works on crashworthiness of composites in terms of two major aspects: the constitutive material model and the finite element (FE) model.  1.2.2.1. Composite constitutive material models  Botkin et al. (1998) [23] summarized the research works on dynamic crush simulation of fiber glass composite tubes. Three different material models for three different materials were evaluated. It was pointed out that all the models were based on non-physical adjustable parameters and thus more work needed to be done for truly predictive capability.  A large amount of composite material models have been implemented in commercial software programs like LS-DYNA and PAM-CRASH. Schweizerhof et al. (1998) [25] reviewed and enhanced three composite material models in LS-DYNA, i.e. MAT 54,   11 MAT 58, MAT 59. It was pointed out that any set of parameters used to fit the experimental results was only valid in combination with a characteristic element size and a particular geometry, not valid for the investigated material in general. Deleo et al. [27] also investigated the major parameters in MAT54 of LS-DYNA on the crash response of UD tape sinuisoidal specimens. This study showed that many parameters including the contact definitions could have a significant effect on the numerical results, thus the simulation could match the experiment but could not predict it. Mamalis et al. (2006) [26] simulated crash performance of CFRP square tubes by using MAT55 in LS-DYNA. Despite of a good correlation under static condition, strain rate effect was claimed to be a necessary factor to consider in order to simulate dynamic crush. Both Obradovic et al. (2012) [29] and Belingardi et al. (2013) [30] analyzed the crash behavior of impact attenuator structure for a formula SAE car. It was shown that the difference between numerical results and experimental ones was due to the complexity to simulate the initialization of the fragmentation phase and the defects in the construction of the tubes.  McCarthy et al. (2001) [24] also investigated the capability of a 'Bi-Phase' model in PAM-CRASH by simulate the crash performance of a composite helicopter subfloor structure. In order to match the experimental results, a significant amount of effort had to be spent on tuning several parameters. It was concluded that the calibration of post-failure material properties was more crucial than the behavior up to initial failure.  In conclusion, a truly predictive composite constitutive material model is still required.    12 1.2.2.2. FE model  A number of studies have attempted to model composite material tubes under axial impact [26,34,47,49,136-138]. However, very few were able to correctly predict both the failure morphologies and force-displacement responses. Crash simulations have the tendency of instability, particularly when it comes to axial crash of composite tubes without plug initiators. McGregor et al. [47] attempted to model the composite tube with a single layer of shell elements, while Xiao et al. [11] used multiple layers of shell elements. Although both works showed a good correlation with the experiment in terms of force-displacement responses, none of them was able to simulate the unique frond-like failure morphology in the crash front as shown in Figure 1-6. Furthermore, these simulations tend to have convergence and global buckling issues.   (a)                                                                           (b) Figure 1-6. Simulation of a [0/±30]2-braided composite tube without plug initiator by (a) one layer of shell elements [47] and (b) two layers of shell elements [11].  There are some other works looked into the simulation of failure morphologies. Mamalis et al. [26] studied different collapse modes with different tube geometries, like the wall thickness and length. Using a quarter-model with three layers of shell elements,   13 simulations were able to produce a stable, progressive splaying of the tube walls under static crushing as shown in Figure 1-7a. However, the predicted force-displacement response did not match the experimental results. Furthermore, a time-step failure parameter has to be adjusted in order to achieve a specific collapse mode. McGregor et al. [22] attempted to improve the stability of the simulation model by introducing fixed debris wedges even though real debris wedges were formed gradually during a crash. Artificial rigid walls were also introduced on the symmetric boundaries through the whole tube to constraint the bended inner plies as shown in Figure 1-7b. Kiani et al. [28] modeled a composite tube using three layers of shell elements. To ensure a correct crush mode, prior to loading, the trigger was modeled in such a way that the first row of the inner-layer elements bent inward while the first row of the outer-layer elements bent outward. A crash front softening parameter SOFT, defined in composite material models in LS-DYNA to reduce the strength of the elements at crash front, was also tuned to obtain a satisfactory result. However, the simulated failure morphology was still not satisfactory as shown in Figure 1-7c.       14  (a)                                          (b)                                        (c) Figure 1-7. (a) Progressive collapse of a CFRP quarter-model in [26]; (b) use of fixed debris wedge and rigid walls to assist the formation of failure morphologies in [22]; (c) final deformation of the simulation model obtained in [28].  In conclusion, a new modeling technique is also required in order to correctly predict both the failure morphologies and force-displacement responses of a composite structure.  1.2.3. Summary  There is a large amount of work done both experimentally and numerically on the crashworthiness of composites. However, due to the complex damage/failure mechanisms of composites, the development of current simulation models is not in pace with the understandings on the damage and failure mechanisms of composite structures. The current simulation models are still not able to correctly predict both the failure morphologies and force-displacement responses of a composite structure. As a result, the slow development of the modeling techniques may largely limit the application of composites materials for mass production in the automotive industry [1, 30]. A robust crash model with real predictive ability is in urgent need.    15 1.3. Research Objectives  This work is dedicated to develop a robust, accurate crash simulation model for composites which is able to both evaluate the deformation behavior and determine the energy absorbing efficiency of various composite structures. The proposed crash model is composed of two major components: a new composite material model and a new composite modeling technique.  For the new composite material model, the main objectives are: 1) Be able to fulfill energy balance; 2) Could derive most input parameters from testing results; 3) Consider both pre-failure and post-failure behaviors; 4) Consider plasticity and any other irreversible strains.  For the new composite modeling technique, the objectives are: 1) Easy to model; 2) Computationally efficient; 3) Be able to correctly capture different failure morphologies naturally.  Both one element and structural simulation tests will be performed to fully evaluate the new crash model. Structural crash tests have been performed on two different groups of triaxial braided composite tubes: one group of composite tubes is tested with plug initiator, and the other group is tested without plug initiator as shown in Figure 1-8. It should be noted that a plug initiator was introduced in some tests [22] to prevent   16 catastrophic failure. With a plug initiator, the composite tube is restricted from bending inward but splitting at the corners and bending outward. The plug provides a smooth guidance which guarantees a stable and highly controllable crash performance. A tube crushed with a plug initiator, however, generally dissipates less energy than the one without a plug initiator as will be discussed in chapter 7.     Figure 1-8.  Schematic of damage/failure mechanisms of braided composite tubes in the crash front with and without plug initiator [22].      17 1.4. Thesis Structure  This thesis is organized in 8 chapters. This section presents a brief overview of different chapters.  Chapter 2 is intended to find out the difficulties of existing composite material models. A detailed review on CDM models for pre-failure behaviors and different methods for post-failure behaviors are presented.  Chapter 3 proposes a new composite material model as Enhanced Continuum Damage Mechanics (ECDM) model. The detailed theory of this constitutive model is described in detail. The general implementation procedures of ECDM model in LS-DYNA is also introduced.   Chapter 4 introduces the 2D3A material used in this research. Validations of ECDM model are performed on the quasi-static one element models and the dynamic tube models with the plug initiator.  Chapter 5 investigates the simulation difficulties of composite tubes without the plug initiator. Three test cases are performed to compare different element formulations in LS-DYNA. A new shell-beam element method is proposed and evaluated.    18 Chapter 6 examines both ECDM model and shell-beam element method by crash simulation of thin-walled composite tubes both with and without the plug initiator. Five different composite tubes are used in this study.  Chapter 7 analyzes the sensitivities of the tube crash performance on studied parameters, including triggers, eratios, tiebreak contact strength, friction coefficients, mechanical properties and residual strength and element deletion strains. An investigation is also done on the detailed simulation results, including damage conditions, residual stiffness, delamination and energy breakdown.  Chapter 8 concludes this thesis and provides outlooks into future researches.  The appendix provides the following information: derivation procedures of yielding & loading surfaces and conditions based on strain increment, detailed ECDM material code and ECDM input cards in LS-DYNA.     19 CHAPTER 2:  COMPOSITE MATERIAL MODELS   This chapter is intended to assess the effectiveness and shortcomings of the existing composite material models and solutions in crashworthiness simulations. After a brief discussion on the requirements for crash simulation, detailed reviews on CDM models and their capability and the methods for modeling the post-failure behaviors are presented.  2.1. Considerations of a Composite Material Model for Crash Simulation  For a crash model, one of the most important requirements is to correctly evaluate and predict the energy distributions. The total energy consists of kinetic energy, internal energy and dissipated energy. The fundamental mechanical behaviors of a material determine how these energies evolve.  To support the development of physics based model for crashworthiness prediction of composite structures, DeTeresa et al [45] measured the post-peak response of a braided composite by loading a specimen in compression till substantial damage developed, then unloading and finally loading to failure in tension. Xiao et al [46] attempted to use a traditional CDM model, MAT 58 in commercial explicit FEA code LS-DYNA, to simulate such material behavior. As shown in Figure 2-1, in experiment, a significant amount of irreversible strain remains in the system after unloading of compressively damaged composite. However, MAT 58 only considers elastic softening due to damage. It does not take the irreversible deformation into account. As a result, the energy   20 dissipated by the damaged material is treated as the elastic energy, which is subsequently released back to the system upon unloading. This model underestimates not only the absorbed energy but also the residual stiffness of the composite materials which increases the tendency towards instability in simulations [11, 34].   McGregor et al. [47], using an in-house CODAM model, significantly improved the crash simulation results by describing a correct unloading path for compressively damaged composites. Xiao [48] developed a simple algorithm to extend an existing CDM model into a coupled damage plasticity model to include the irreversible strain. The model was capable of modeling substantially damaged materials and improved crash simulations. This model, however, has some shortcomings. For example, it has the same CDM framework with MAT58 where the damage is described with an exponential power law that coupled the pre-failure and post-failure responses. Such coupling effect makes it hard to precisely describe both the pre- and post-failure responses. Although the models with irreversible strains significantly improved the stability of crash simulations, the robustness of the models is still unsatisfactory [22, 48-49].     21  Figure 2-1. Comparison of the stress-strain response predicted by MAT 58 to experiment data of a [0/±30]2 triaxial braided composite. The specimen was loaded in compression, then unloaded and loaded to failure in tension [11].   In summary, to correctly predict the crash performance of a composite structure, the constitutive material model should be able to correctly describe the plastic deformation of the matrix and any other irreversible behaviors of the composite. Decoupling of pre-failure and post-failure behaviors is also required.  Among all the approaches used to describe pre-failure behaviors, CDM models are the most-widely used ones due to their capabilities of modeling progressive damage/failure of composites as pointed out by both Tay et al. (2008) [53] and Liu et al. (2010) [54]. Numerous CDM models are available in literature. This chapter reviews CDM models for pre-failure behaviors and different methods for post-failure behaviors separately.      22 2.2. Literature Review on CDM Models 2.2.1. Concept and Characters of CDM Models  CDM was initially introduced by Kachanov (1958) [55] to predict creep rupture of metals, and then further developed by Hult (1979) [56], Chaboche (1981) [57] and Krajcinovic (1984) [58]. It provides a measure of the micro-scale material degradation or damage at macro-scale, thus couples -strain behavior. It can be used to account for the continuous stiffness degradation and to predict 54-61]  The general concept of CDM can be found in [57-58, 60, 62-63]. The stress-strain relationships can be rigorously derived on the basis of the thermodynamics with internal variables. It has been largely developed in recent decades. The differences between various CDM models are: the choice and definition of damage variables, formulation of the energy potentials, damage evolution laws and the considered phenomena (creep, plasticity, fatigue, anisotropy, etc.) [50, 59, 61].  A lot of damage models are based on principles of thermodynamics, because the inelastic mechanical behavior of engineering materials is intrinsically linked to irreversible thermodynamic processes accompanying energy dissipation and physical changes of the microstructure. A CDM model typically consists of three major components: (a) definition of damage variables, (b) definition of strain energy density (as a function of damage variables and other state variables) and (c) definition of damage evolution laws. [64] The general relationships are available in [50] and [54].    23  2.2.2. Popular CDM Models for Composites  Many CDM models are available in literature. Kaddour et al. (2013) [65] compared twelve different damage/failure models for continuous fiber-reinforced polymer composites under multi-axial loading in the third world-wide failure exercise. Table 2-1 summarized some well documented CDM models. The following is a short review of these material models.  Murakami and Kamiya [66] developed a damage model based on irreversible thermodynamics by assuming that Helmholtz free energy was a quadratic function of elastic strain tensor and second rank symmetric damage tensor.  Schapery [67] proposed a thermodynamically based work potential theory to model material damage. This theory was developed for structures that exhibit limited path-independent work so that the total applied work can be partitioned into elastic work and the work of structural change, both of which were functions of internal state variables. Pineda et al. [68] applied s Theory (ST) to model progressive microdamage along with generalized method of cells (GMC) to determine the nucleation of transverse cracks and the truncation of microdamage, as well as fiber failure. One single internal state variable was assumed to account for all damage presented in the matrix up to mesoscopic matrix cracking. Damage evolution law was found by fitting polynomials with experimental data of degraded stiffness. Instantaneous transverse matrix cracking and fiber breakage was assumed. This theory was further improved by Pineda and Waas   24 [69] by introducing multiple internal state variables, i.e. transverse matrix cracking, shear matrix cracking and axial fiber fracture.  Talreja [70-71] developed a CDM model based on a set of vector fields, each representing a damage mode. This approach was called synergistic damage mechanics (SDM) since it is enriched with micromechanics [51]. The material damage had two contributions: the damage entity vector, represents the type of damage in the material; and the damage influence vector, accounts for the influence of the flaw on material behavior [50]. As a result, each damage mode was depicted by a second order tensor. The damage evolution law was formulated based on the evolution and increase of flaws. This model automatically considers crack opening and closing with the change of loading directions. However the process of determining certain properties, like crack opening displacement and crack density, were very tedious [50]. Since this method requires tracking individual damage mechanisms, it becomes formidably difficult when different damages start to interact. As a result, this model is unable to predict the lamina or laminate failure [51].  Puck et al. [72] developed a damage model by introducing empirical degradation laws to This model was enhanced by Schuecker, Pettermann and Flatscher [50, 73] by considering plastic shear deformation. Pinho et al. [74] also proposed a 3D physically based failure criterion for laminated fiber-reinforced 7275] work. The matrix compression failure was based on the MohrCoulomb criterion. This model was able to predict the fracture   25 angle. Due to further rotation during the compressive loading, an initial fiber-misalignment angle was used to consider fiber kinking and thus to trigger the failure. Pinho et al. [76] implemented this model in LS-DYNA. The damage initiated at the onset of failure and then evolved linearly until final failure. Later, the damage evolutions were modified by correlating them with the fracture angle so that both wedge effect and progressive cracking could be considered [77]. The parameters used in this model can be identified from tests. This model was evaluated in WWFE-II [77] and WWFE-III [78].  Williams et al. [79] developed a CDM model known as the COmposite DAmage Model (CODAM). This is a sublaminate-based constitutive model with damage parameters defined in the sublaminate coordinate system instead of principal lamina coordinate system. The sublaminate behavior was decoupled into fiber and matrix, both of which assumed linear response between damage initiation and damage saturation. A damage potential function was also assumed, which however, introduced many parameters that were difficult to characterize. The original CODAM model only considers elastic softening. To model the irreversible strains needed in crash simulations, McGregor et al. [80] improved CODAM by introducing an 'analog model' to describe the unique compressive damage behavior. The analog model included two sub-models: Laminate sub-model and Rubble sub-model. The Laminate sub-model depicts the behavior of the intact laminate during damage propagation and results in overall strain-softening response in both tension and compression, while the Rubble sub-model describes the compressive behavior of any material that is fully damaged but still able to transfer some loads. The damage evolution was assumed to evolve linearly between the damage   26 initiation and saturation. CODAM was initially designed for laminated composites. It is questionable to be applied on braided composites [47] since the behaviors of braided composites are too complicated to be divided into independent fiber and matrix behaviors. The corresponding physical foundation is doubtful.  52, 81] developed a CDM model to predict the onset and evolution of intralaminar failure mechanisms. LaRC04 failure criteria [75] were used as damage activation functions. Four ply fracture modes were considered: tensile fracture, fiber kinking and matrix cracking with fracture planes oriented at 0 and 53 degrees with respect to the ply thickness direction. Unilateral effect was considered by defining different damage variables corresponding to active and passive damage in both longitudinal and transverse directions. Exponential laws were used to describe the evolution of the damage variables. For longitudinal fracture, Bazant and Planas's model [144] were used to define two damage evolution laws: fiber bridging with linear softening and fiber pull-out with exponential softening. As a result, the tensile damage in longitudinal direction can decrease the compressive stiffness corresponding to the misalignment of fractured fibers, while compressive damage in transverse direction can increase tensile damage related to larger crack. This model is also able to consider thermal effects.  Ladevèze model [82], was developed by the research group of Ladevèze at LMT Cachan, France [82-91], is one of the most widely used CDM models for fiber reinforced composites based on energy potentials. Uniform state of damage was assumed within   27 each meso-constituent [84, 50]. It was proved that this mesomodel can be interpreted as the homogenized result of micromodels involving common microdamage mechanisms like mircracking, fiber-matrix debonding and delamination. Since it provides a general formalism, it can be transposed more easily into commercial codes [92]. Ladevèze model took into account stiffness recovery and inelastic strains [50]. It was sufficient to describe the nonlinear or plastic behaviour that some thermoset or thermoplastic composites exhibit, especially under transverse and shear loadings [93]. All the parameters needed in the elementary ply of a Ladevèze model can be measured by experiment as listed in [82]. 94] carried out a detailed experimental test series to determine the input parameters for Ladevèze model on both unidirectional carbon fiber reinforced plastics and glass fiber reinforced plastics. By using these measured parameters, it was proved that Ladevèze model can capture the full extent of the non-linearity and accurately predict the responses in axial, transverse and shear directions for both materials.  Ladevèze model are not limited to a specific loading and configuration. Allix et al. extended Ladevèze model to consider compression [95] and temperature-dependent scenarios [96]. Phillips et al. [97] demonstrated that Ladevèze model was able to predict damage regardless of fiber orientations. Johnson et al. [98-100] modified Ladevèze model by decoupling the responses of transverse and shear directions. It was shown that this modified model was able to predict the impact performance of fabric reinforced composites. Pickett et al. [101-102] further improved modified Ladevèze model by changing the evolution law of shear damage in terms of shear driving force   28 from linear to polynomial function of rank two so that it could be use on biaxial braided composites.   Table 2-1 compared the detailed characters of the aforementioned CDM models. In summary, the material characterization process is rather complicated for some material models like CODAM [4751, 70-71]. Most CDM models ignore the existence of irreversible strains in composites. By comparison, the Ladevèze model [82] appears to be a good candidate as a starting point of a crash model. Firstly, as a thermodynamic phenomenological model, it is able to predict consistent thermodynamic energy. Secondly, its inputs are all measureable scalar damage parameters. Last but not least, plasticity is considered, which is very important to correctly capture the energy dissipation during damage and failure scenarios.  2.3. Literature Review on Post-Failure Models  A variety of materials, such as paper, some tissues and fabrics, and fiber-reinforced composites, are able to sustain certain load after losing the initial strength, i.e. their stress-strain curves exhibit a post-failure region. This region is marked by the decline of stress with increasing strain or, generally, when the tangential elastic modulus ceases to be positive definite. This behavior has also been observed in rocks, concrete and some soils under compression [132-133,135,142]. The heterogeneity and brittleness of these materials is the cause of strain softening. The mechanisms consist of progressive distributed damage, such as dispersed microcracking, void formation or loss of interparticle contacts. For a crash model, the post-failure behavior can affect the stability of a simulation model and, more   29 importantly, the predicted SEA to a large extent [11, 24, 34, 46-49]. This section reviews some of the notable works on strain-softening.  2.3.1. Experimental work on strain-softening  The experimental works on strain-softening of composites are concentrated on compression scenario. For laminated composites, fiber kink banding has been identified as a compressive strength limiting mechanism in aligned fiber reinforced composite. [107] However, most of the experiment is only able to capture the phenomenon up to the peak load, the so-called failure. The behavior after the peak load, i.e. the post-failure region is largely unknown.   One of the few experiments designed to study post-failure behavior is the Overheight Compact Tension (OCT) test developed by Kongshavn and Poursartip [108]. This test investigated the damage evolution in front of the crack tip by sectioning the process zones into thin strips roughly 2mm wide. It was found that the failure strain of the damaged material can exceed the fiber failure strain by as much as a factor of two and a strain softening curve is necessary to capture the progressive growth of damage across the specimen width during loading.  Zobeiry et al. [109] recently proposed a new compact compression test method to get the strain softening response of a quasi-isotropic laminate composite [90/45/0/-45]4s. This test was able to capture the main compression failure modes and their evolution process. It was concluded that during compressive failure, the composite laminate went through   30 the softening and broadening steps before the ultimate failure due to the formation of slip surfaces. However, the obtained strain-softening responses were too scattered to generate a persuasive optimized response.  2.3.2. Modeling work on strain-softening  The ideal way to model strain-softening is to build a micromechanical model that is able to include all the detailed scenarios so that the strain-softening behavior can be naturally captured. However, most micromechanical models like kinking band model are designed to predict the measurable strength only and they are strongly sensitive to some adjustable micro scale geometries. Basu et al. [110] examined the compressive strength of a fiber reinforced lamina under multi-axial stress states using a new model for fiber kink band. It was shown that the in situ compressive strength of a lamina is strongly dependent on the initial fiber misalignment, load multi-axiality and the complete, nonlinear lamina shear response. Feld et al. [111] also pointed out that the nature of the laws and couplings introduced in the kink band model depend strongly on the target meso-model. Thus the micromechanical models are not mature enough to support the strain-softening modeling.  By comparison, phenomenological models are more practical and workable to model strain-softening. In Nahass review [106], twelve different methods were compared for their abilities to model the post-failure behavior of laminated fiber-reinforced composites beyond first ply failure. Hahn-Tsai [112] assumed that any lamina would still support the load which it was carrying when it failed until the total laminate failure occurred. Chiu [113] considered the case of instantaneous unloading. The same concept was used by   31 many failure criteria, like the maximum stress, maximum strain, Hashin, Tsai-Wu, Tsai-Hill, etc. Petit-Waddoups [114] assumed the failed lamina unloaded linearly. Nahas [115] used an exponential function to depict the unloading curve. It was pointed out that the failure criterion and nonlinear considerations were important factors that affect the post-failure behavior. (Exponential softening was also used in (Flatscher, Pettermann 2011) [116].) However, these models have the drawbacks of being effective only for specific types of laminate [117] or being unrealistic as it disregards the constraints that are imposed on the failed lamina by the adjacent lamina and the undamaged part on the damaged part [118]. Mazzeranghi and Vangi [117] proposed a post-failure model that considered some detailed damage phenomena. But there were some very fundamental assumptions need to be made, like interface thickness, crack thickness etc. It is worth noting that these models are all designed to evaluate the mechanical characteristics of a laminate with failed lamina.   The post-failure strain-softening behavior can be more rigorously modeled by using CDM. In CDM, the evolution of damage is determined from assumed functional relationships which include the current state of damage. Thus, softening behavior is always occurring after initial failure with increasing strain. Randles [119] assumed that the damage threshold at a finite stress value and then increased in a simple form as:                                                   (2-1)    32 where  is the initial youngs modulus,  is a damage parameter. The ultimate failure was identified as a complete loss of load-carrying capability through loss of modulus as the coefficients in front of the initial modulus approaching zero. This model was further evaluated by Nemes and Speciel [120].  A more famous strain-softening model was developed by Matzenmiller et al. [121], referred as MLT (Matzenmiller, Lubliner, Taylor) model. The damage variables in this the strain-softening behavior naturally. This model is widely used [23] and implemented in LS-DYNA as MAT 58 [44] and extended by Schweizerhof et al. [25]. Kongshavn et al. [108] demonstrated significant improvement in the predictions of damage progression in notched compact tension specimens by using MLT model. Williams et al. did a thorough investigation on this model [118, 122]. It was proved that MLT model is able to predict impact failure in FRC laminates more accurately. However, the greatest difficulty in applying this model is to obtaining the characterization data for the damage growth law.                                              (2-2)  where  is the damage parameter, m is a value obtained by fitting measured stress-strain curve,  is the strain,  is the strength.     33  Figure 2-2. Effect of variations in the exponent, m, on the longitudinal stress-strain behavior predicted by the MLT model. [118]  As shown in Figure 2-2, as m approaches infinity, the predicted response approaches that of the instantaneous failure models. As m decreases, not only does the strain-softening behavior become more gradual but also the damage growth prior to the peak load  become increasingly non-linear in the loading part of the curve [118]. Even though these behaviors cannot be independently determined or calibrated [108], the strain-softening model was continuously been used [123-129]. William et al. pointed out that the coupling of pre- and post-failure behavior of MLT model largely limited its improvement on mesh sensitivity and some other issues like strain rate dependence. Yen [123] overcame some of these weaknesses of the MLT model by: 1) separating the elastic loading and post-failure softening into two models; and 2) considering the layer strength and elastic moduli as logarithmic functions of the strain rate. This model was implemented in LS-DYNA as MAT 162 [124]. Deslauriers et al. [125], Brown et al. [126] and Xiao et al.   34 [127] demonstrated that MAT162 was useful on modeling the progressive failure of composite materials consisting of unidirectional and woven fabric layers under impact. Xiao et al. [127] even established a methodology to systematically determine the model function to quantify damage evolution and corresponding stiffness reduction were adopted in many other studies, like (Tabiei, Aminjikarai 2009) [130].  In summary, it is reasonable to use two different damage models to depict the pre-failure and post-failure behaviors since they are undergoing different damage mechanisms. Ladevèze model is ideal to be used as pre-failure material model while distribution function can be used to describe post-failure behaviors.      35  Table 2-1. Comparison of different CDM models. Model Damage form Damage type Unilateral effect* Damage evolution law Tension/ compression Axial damage Ther./ Phe.*** Irreversible strains Comments CODAM [47] scalar Damage No Linear Yes Yes Fake Phe. Yes  Xiao [48] scalar Damage No Power No (different strength) Yes Phe. Yes Coupled damage-plasticity model Ladeveze [82] scalar Damage Yes Exponential Yes Yes (compressive) Ther.+Phe. Yes  Maimi [52] scalar Damage Yes Polynomial Yes Yes Ther.+Phe. No Thermal regularization MLT [44,121] Scalar Damage No Power No (different strength) Yes Phe. No MAT 58 in LS-DYNA Murakami [66] 2nd order tensor Damage Yes - Yes Yes Ther.+Phe. No Based on assumed form of Helmholtz free energy Pinho [74,76-77] scalar Cracking No Linear Yes Yes (kink) Phe. No Smeared softening Schapery [67,68] scalar ISV** No Polynomial No No Ther.+Phe. No  Talreja [51,70-71] 2nd order tensor Cracking Yes - Yes No Ther.+Phe. No Complicated crack quantification  *Unilateral effect: the closure of transverse cracks under load reversal [50, 52]. The damage mode is either active or passive. **ISV: internal state variable. ***Ther.: thermodynamical model; Phe.: Phenomenological model.  Phenomenological models in nature require experiments to characterize the damage parameters, not derived by consideration of physical modeling.      36 CHAPTER 3:  ECDM MODEL  This chapter is dedicated to developing a robust, accurate material model for crash simulation of composite structures. A new composite material model, the so-called enhanced continuum damage mechanics (ECDM) model is proposed. The theories of the ECDM model are introduced in the following order: pre-failure material model, initial failure criteria, strain softening laws, post-failure residual model and final failure criteria. The general implementation procedures of ECDM model in LS-DYNA is described at the end of this chapter.   3.1. Proposed ECDM Material Model  The proposed ECDM model employs two sub-models, a pre-failure model and a post-failure model, to describe the pre-peak and post-peak regions of the stress-strain curve. Figure 3-1 is a schematic of uniaxial tensile stress-strain response of ECDM in one material direction. The pre-failure model uses a modified Ladevèze model to describe the nonlinear behaviors of the material before it reaches its peak strength. The post-failure model covers two stages: a strain softening and a residual state stage. Two different laws are employed to describe the property evolution in these two stages. There are two criteria in the ECDM model: an initial failure criterion and a final failure criterion. The initial failure criterion determines the peak strength, while the final failure criterion determines the element deletion conditions.   37    Figure 3-1. Schematic of the uniaxial stress-strain response of the ECDM model. The model is composed of two sub-models: a pre-failure model and post-failure model. The post-failure model considers two stages:  strain softening and residual state with two different property evolution laws.  As is mentioned in Chapter 1, post-failure behaviors like strain-softening occur not only in tension, but also in compression and shear [131]. This study considered them in all material directions under different loading conditions. A detailed description of each sub-model is presented next.  3.2. Pre-failure Material Model  The pre-failure model uses a modified Ladevèze model to describe the nonlinear behaviors in pre-failure region. The original Ladevèze model was designed for unidirectional composites [82]. To predict the impact performance of fabric-reinforced   38 composites, Johnson et al. [98-100] implemented a modified Ladevèze model where the responses of transverse and shear directions were decoupled and the plasticity is considered only in the shear direction. The same approach is adopted in the current study.   Figure 3-2. Elementary ply [82]  3.2.1. Damage kinematics of the elementary ply  ECDM model only considers damage at the elementary-ply level to describe the matrix microcracking and fiber/matrix debonding. Delamination is not considered inside the material model. As shown in Figure 3-2, a plane-stress state is assumed. The damaged material strain energy is written as:                                                                                 (3-1) with     39 where   and  are the scalar damage variables that remain constant throughout the ply thickness; ,  and  are the stresses; ,  and  are the moduli of elasticity; f can be any variable such as E; the subscript t denotes tension and c for compression.   The tensile behavior and compressive behavior are described separately since they have different damage and failure mechanisms. Under compression, as is shown in Figure 3-3, the formation of kink bands is generally unavoidable for continuous fiber-reinforced composites. The waviness of the tows in fabric or braided composites can trigger the formation of kink bands, as shown in Figure 3-3b. Meanwhile, the weak fiber matrix interface allows fiber movement inside the composites [103]. These two factors result in micro-buckling, which is not seen under tensile loading.  The damage forces, , ,  , are conjugate quantities, which are defined by:    (3-2)    40  a) b)  c) Figure 3-3. Under compression, a) kink bands formed in a unidirectional carbon reinforced composite [104]; b) kink bands formed in a 2D3A material [105]; c) micro-buckling observed in a 2D3A material [103]   The damage forces are analogous to the energy release rates in Fracture Mechanics. They govern the damage development as the energy release rates govern the crack propagation. It was shown in [82] that the square root of the damage forces can be quantified much easily from testing data, thus the following parameters are defined. They are called  damage driving forces in this study.   41                                                                                                                                    (3-3)  where represents any previous time. The damage evolution laws are written as:   if   if   if                                              (3-4)  , ,, , ,  and b are material characteristics which can be determined by experiment. As shown, the material behaviors in all directions are decoupled.  3.2.2. Plasticity modeling and damage-plasticity coupling  In order to model the inelastic or irreversible deformation of a composite ply, plasticity is considered in the damaged material. Experimental results in longitudinal and transverse directions of triaxial braided composite in pre-peak region show almost linear elastic behavior in tension and compression, with negligible plastic strains, and thus plasticity is   42 only considered in the shear direction. A classical plasticity model is used with an elastic domain function and hardening law applied to the effective stresses in the damaged material. Inelastic or strain increments are assumed to be normal to the elastic domain function [82, 99].  The elastic domain function is defined by:                                                (3-5)  where  is the initial threshold value for inelastic strain behavior, p is the accumulated effective plastic strain over a complete loading cycle,   is a material characteristic function determined from cyclic loading tests.  corresponds to a stress state inside the elastic domain where the material may be purely elastic or elastic damaging. Generally, an index function is assumed as the hardening function:                                                                (3-6)  where  and power index m are material properties which can be determined by cyclic tests. The plasticity is assumed to evolve isotropically due to the isotropic nature of the matrix. It should be noticed that this plasticity theory is only used to define pre-failure irreversible strains. Post-failure irreversible strains are calculated using a separate theory as is stated in 3.5.2.    43 3.3. Initial Failure Criteria  Failure occurs in a material when the applied load reached a threshold that is the limit of its loading carrying capacity. Since it is an expensive procedure to establish strength characteristics of materials for every complex stress state, the idea of "failure criterion" or "failure theory" has been introduced to predict the strength of materials under multiaxial loading conditions using strength data obtained from uniaxial tests. The failure criterion (being a function of the stresses, is a surface called "failure envelope" determined in the stress space) encloses the stress states that the material can sustain without failure. In biaxial stress state conditions this envelope is a curve [106]. To be consistent with the damage evolution law, the axial failure criterion is independent from the other two directions. The same as in the MLT model [44, 121], a set of simplified Hashin criteria is used:   Axial direction:                                      (3-7) Transverse and shear direction:              (3-8)  where ,  and are the axial strength, transverse strength and shear strength, respectively. The transverse direction and shear direction are coupled in this model since these two directions share some common damage modes like matrix cracking. Thus, the damage in one direction will affect the damage in the other direction.    44 3.4. Strain Softening Laws  The ability to smoothly describe the sudden drop of load-carrying abilities. It has been used to quantify the damage evolution and corresponding stiffness reduction in other studies [121, 130]. In the ECDM model, the pre-failure damage evolution law and post-failure damage evolution law are separated due to different damage mechanisms. The following damage evolution law is used for post-failure damage:                                              (3-9)  where  and  are the damage and strain at the peak load, i.e. the beginning of the post-failure region, respectively, and m is an index used to control the falling slope of a stress-strain curve. The slope drops more quickly with a higher value of m. It can be shown that before initial failure, while  after initial failure. Thus, the damage state is continuous even though the pre-failure damage evolution law and the post-failure damage evolution law are decoupled.  3.5. Post-Failure Residual Model 3.5.1. Post-failure residual state model  In a displacement controlled experiment or in a structure that permits local unloading and load redistribution before catastrophic failure such as under crash and impact loading, a   45 material can exhibit a residual load carrying capability after initial failure.  Some models, like MAT 58 in LS-DYNA, give a fixed residual stress value to model this phenomenon. A closer examination of experimental results reveals that, for some materials such as braided composites, the residual load carrying capability may still increase a little before complete failure. Figure 3-4 explains this phenomenon. Assuming 90% of the fibers in a braided composite rupture at a given point, the remaining 10% can still take load. The load carried by the unbroken fibers increases with deformation till final rupture. The overall behavior of the composite is the summation of these two groups of fibers, Figure 3-4b.     (a)                                                               (b) Figure 3-4. Schematic of the phenomenon of increasing residual load carrying capability with deformation. (a) Assume 90% of the fibers fail earlier than the remaining 10% fibers; (b) the response of the material is the summation of the two responses in (a).   In the ECDM model, two options are available for the description of the post-failure residual state. The first option uses a constant residual stress, while the second one uses a constant residual stiffness in each major direction. For the first option, the minimum stress in the   46 post-failure region of the material stress-strain curve is the residual stress. The corresponding damage state can be described by:                               (3-10)  For the second option, the maximum damage in the post-failure region is the residual damage determined by residual stiffness. The corresponding damage state can be written as:                                  (3-11)  where  is the residual stress and  is the residual damage. It should be noted that force-displacement responses were used to determine these residual values instead of stress-strain responses. This is because most strain measuring techniques, such as strain gage, laser extensometer and digital image correlation (DIC), are unable to provide accurate results after initial failure of the composites.  3.5.2. Post-failure irreversible strains  One of the most important functions of a crash model is to correctly evaluate and predict the energy distributions. Irreversible strain energy is one of the most important energies that is neglected by most crash models. The irreversible strains are classified as the pre-failure   47 irreversible strain and post-failure irreversible strain as is shown in Figure 3-5. Without considering any one of them, the simulation model may underestimate the energy absorption ability. Pre-failure irreversible strain is considered by the plastic model discussed in 3.2.2. This section focuses on modeling the post-failure irreversible strain.  For a material under tension, its residual state is generally determined by its residual stiffness and thus the incremental strain in the post-failure region is considered to be reversible. For a material under compression and shear, the microbuckling and dislocation of the fibers are permanent and thus most of the corresponding incremental strain is considered as irreversible strain. Figure 3-5b compares the simulation results of a CDM composite model, MAT58 in LS-DYNA [44] with LLNL compression tension test data [11]. Since a major portion of the strain accumulated after initial failure is irreversible, using the traditional CDM models may underestimate a large portion of the absorbed energy (irreversible strain energy) if the material is unloaded in the post-failure region.  Generally, a strain can be sectioned into two parts:                                              (3-12)  where  and  represent the reversible and irreversible strains, respectively. In the post-failure model, a set of new constants ERATIOs is introduced to divide these incremental strains. ERATIOs are defined as follows:    48                                             (3-13)  ERATIOs are user-defined inputs. They have a significant influence on the bending curvature of crash front material. It is hard to measure these parameters by tests due to severe damage of the material after initial failure. However, they can be determined by trial and error by matching the crash front shapes between simulations and experiments.   (a)                                                                    (b) Figure 3-5. Energy comparison between (a) CDM model and plastic CDM model; (b) predicted results by MAT58 in LS-DYNA and LLNL compression tension test data [11].   3.6. Final Failure Criteria  In explicit FE simulations, the final failure criteria are also the element deletion criteria. As far as a volume of material is able to carry some residual load, it should not be deleted from the model. On the other hand, element deletion is used to simulate the separation of   49 structures. Final failure criteria are purely designed for numerical models without a sound physical foundation; thus, simple equivalent strains are used to build these criteria. Two options are provided in the ECDM model: option 1 uses individual strain failure criteria, while option 2 uses an equivalent strain failure criterion, as shown below:  Option1: Axial direction:        Transverse and shear direction:          (3-14)  Option 2:             (3-15)  where , ,  and  represent the axial, transverse, shear final failure strain and equivalent final failure strain, respectively.  3.7. Implementation of ECDM Model in LS-DYNA  The ECDM model is written as a user-defined material model for LS-DYNA. In LS-DYNA, it is the strain increment that is passed into the user-defined material subroutines; thus, all the equations are rewritten and organized based on strain increment. The general working flow chart is shown in Figure 3-6.    50  Figure 3-6. The general working flow chart of the ECDM model implemented in LS-DYNA.        51 CHAPTER 4: VALIDATION OF ECDM MODEL  This chapter is dedicated to evaluate the predictive ability of the ECDM model. Both quasi-static one-element tests and dynamic tube crash tests were simulated. For the dynamic tube crash tests, composite tubes with a plug initiator were used because such tests are relatively more stable as compared with those without a plug initiator [11, 46]. A commonly used CDM model, MAT 58 in LS-DYNA, was used for comparison in both evaluations.  4.1. Material and Material Properties 4.1.1. 2D3A for quasi-static tests  The 2D3A used for the quasi-carbon fibers and and the bias tows contained 12K fibers. This composite has been evaluated for aerospace applications. The identification procedures of the damage parameters are detailed in [82, 94]. To obtain the plasticity parameter, cyclic loading-unloading tests are required. In this study, the tensile experiments were performed with specimens cut in 0°, 45° and 90° directions in respect to the axial tow. All the tests were performed with an Instron load frame (model 1321). A laser extensometer and digital image correlation (DIC) were used to measure the strains. For DIC, ISTRA software was used for image acquisition and analysis, the same as in [134].  Detailed input data for 2D3A (quasi-static) are listed in Table 4-1. As is shown, only tensile tests are done for this material.   52  Table 4-1. Input data of 2D3A for quasi-static simulation. Material Properties Units 2D3A Experiment    48.8  Longitudinal tensile test on 0 specimens  - 0.32 Y1o  0.01943 Y1c  0.7092 X11t GPa 0.632   23.86    Cyclic tensile test on 45 specimens   0.009568   0.139074   0.005   34.269  - 0.7685 X12s  0.135   46.353   Cyclic tensile test on 90 specimens   0.008748   0.16569 X22t  0.302  4.1.2. 2D3A for dynamic tests  Triaxial braided composite tubes supplied by General Motors Corporation are used in this study. They are made of Fortafil #556 80k carbon axial tows and Grafil 34-700 12k carbon biaxial tows in an Ashland Hetron 922 resin with [0/±45] braid architecture [10]. Five different tube configurations are investigated in this study. Special attention is paid to [0/±45]2 and 4 tubes since they are tested both with and without plug initiators. The other three configurations, 2×42 , 4×42 and round 2 [0/±45]2 are only tested without plug initiator. The material properties are derived from [10].  The dimensions of composite tubes tested are listed in Table 4-2.   53  Table 4-2. Dimensions of different composite tubes  LxW (mm) T (mm) 2×22 355.7 x 55.64 x 55.64 2.41 (1.205/layer) 2×24 355.5 x 55.55 x 55.55 6.18 (1.545/layer) 4×4[0/±45]2 356x102x102 3.50 (0.875/layer) 2×42 365x50.95x101.84 2.77 (0.693/layer)  LxD (mm) T (mm) round 2 360.5x63.60 2.04 (0.51/layer)  Detailed input data for 2D3A (GM) are listed in Table 4-3. For 2D3A (GM), much larger strength values are used. This is because the tested tubes are closed structures. Traditional coupon tests tend to underestimate the material properties of composites. Pickett et al. [101] studied the difference between biaxial braided composites with cut and uncut edges. The preparation processes of uncut and cut specimens are shown in Figure 4-1a. The mechanical response is shown in Figure 4-1b. The initial mechanical properties were found to be the same, but the uncut specimen had much larger strength (27.5%) and failure strain (39.1%). This is because they have different damage modes: tow breakage initiated at the stress concentration area adjacent to the tabs for uncut specimens while tow pull-out and tow failure initiate at the free edges for cut specimens. The differences between uncut and cut specimens are also related to the braiding angle which can change the proportion of load transferred by fiber mechanisms and that of shear matrix mechanisms.   54     (a)   b) Figure 4-1. a) Schematic of the preparation of uncut and cut specimens; b) comparison of results for the cut and uncut [25]s braided specimens. [101]    55 McGregor et.al [10] found that wide four-point bending specimens are less sensitive to edge effects. To estimate the tensile/compressive strengths based on four-point bending tests is able to obtain more accurate results. The same approach is used in this study. Generally, cyclic tensile tests on 45 specimens are required to determine the shear material properties and the plasticity parameters. Without such tests in this study, a cyclic tensile stress-strain curve was assumed. It will be proven in 4.3.2.1 that pre-failure irreversible strains are much less important than post-failure irreversible strains, thus such assumption wont affect the results too much.  Table 4-3. Input data of 2D3A (GM) for dynamic simulation Material Properties Units 2D3A (GM)   Material Properties Units 2D3A (GM)    62.5    59.6   0.6568    0.3456  - 0.012   - 0.007   0.1314    0.3110   0.030    0.016   0.286    0.1327  - 0.36   - 0.35     Eratio1 - 0.75   7.52    0.01   0.07    1.6355  - 0.014   - 0.4345   0.035   - 0.5   0.0115  Eratio12 - 0.75   0.05275       12.5    14.4   0.1100    0.0859  - 0.012   - 0.009   0.0550    0.04295   0.01164    0.0038   0.080    0.0604      - 0.6     Eratio2 - 0.75 * if: initial failure ** FF: final failure *** Res: residual   56 4.2. Validation of ECDM Model on Quasi-Static One Element Model  This section focuses on quasi-static test of ECDM model on a single shell element.   a) shear damge response  b) shear plasticity response Figure 4-2. Shear damage and plasticity responses and fitted curves.    57 4.2.1. Growth of damage and plasticity  Special attention is focused on the shear direction which undergoes both damage and plastic deformation. In this study, cyclic tensile tests are performed on 45 2D3A specimens to obtain the required shear parameters. Figure 4-2 is obtained by following the data acquisition procedures. As shown, linear function fits the relationship between shear damage and shear damage driving force very well, and the power law function perfectly fits the growth of plastic strain.  The plasticity parameters are derived from the fitted curve in Fiure 4-2b. In ECDM model, the plasticity is assumed to grow independently.   4.2.2. Tensile responses of 2D3A  After obtaining all the elastic, damage and plastic paramters, detailed simulation are performed in three major material directions by using ECDM and MAT58 material models.   For the axial direction, with negligible plasticity, no cyclic tensile tests are perfomed on the specimens in this direction. The corresponding material properties are directly derived from Littells [103] tensile test results. Failure and post-failure behaviors are assumed in axial direction. For transverse and shear directions, cyclic tensile tests on 45 specimens and 90 specimens are used to derive the damage and plastic material properteis. The   58 cyclic tests are highly repeatable in both directions and they merely change the overall mechanical responses.  Figure 4-3 compares the simulation results with experimental results for simple tension with 0°, 45° and 90° specimens. Generally, the simulation results match the experimental results very well for the pre-failure stress-strain responses. The major difference is observed in the post-failure region. The ECDM model is able to accurately simulate the material behaviors, while MAT 58 largely overestimates the stress responses. The difference between these two models is due to the damage evolution law. As an example, Figure 4-4 compares the predicted damage evolution by MAT 58 and the ECDM model in the longitudinal direction. MAT 58 uses a single damage evolution law for pre-failure and post-failure material behaviors, while ECDM uses three different damage evolution laws in different regions. In experiments, the composite loses a large portion of load-carrying ability after reaching the initial failure point, which corresponds to a sudden increase of damage. The ECDM model is able to simulate such a physical phenomenon. It was noted that, with the residual stiffness, the ECDM model can better fit the post-failure behavior. On the other hand, MAT 58 predicts a significant amount of extra energy that is not really absorbed by the composite materials. This is effect will be further discussed in 4.3.     59  (a) 0 results   (b) 45 results Figure 4-3. Comparison of simulation results of ECDM model and MAT58 and experimental results on  (a) 0, (b) 45 and (c) 90 specimens.     60 Figure 4-3   (c) 90 results    Figure 4-4. Comparison of predicted damage evolution in longitudinal direction by MAT58 and ECDM model.   61  4.3. Validation of ECDM Model on Dynamic Crash Model with Plug 4.3.1. Simulation model for composite tubes with plug initiator 4.3.1.1. Dynamic crash model  The majority composite material models were developed for shell element due to its efficiency and practical potential to be used on large structures. If solid elements are used, at least 4 layers of elements are needed in order to simulate transverse shear and bending deformation which is computationally too expensive. Flesher et al. [136] developed a 3D crash model for braided composites by using two layers of solid elements. Despite good correlations on force-displacement responses of plug-initiated crash, this model had difficulty predicting the crash morphology and was unable to predict no plug-initiated crash.   Shell elements are used in this study to model the composite tube. Previous studies like Mamalis et al. (2006) [26], Zarei et al. (2008) [32], Huang et al. (2009) [33] using MAT54, Xiao (2009) [34] using MAT58 and McGregor et al. [47] using CODAM have shown that model one layer of shell elements for each layer of composites is able to simulate the delamination and splaying of textile composites tubes. Thus this study will initially use the same strategy. But as will be discussed in chapter 5, shell elements are not well designed to simulate the crash behavior for composite tubes without plug initiator. A new modeling technique will be introduced.     62 Four-node fully integrated shell elements are used to avoid hourglass issues (zero energy deformation modes). Figure 4-5 shows the detailed FE crash model built according to the drop tower crash test setup at University of British Columbia (UBC) [22]. The FE model is composed of four major components: a drop mass, a composite tube, a plug initiator and a rigid wall. The drop mass represents the mounting plate and all attached weights on the top with a total mass of 535kg. Two layers of shell elements were used for two-layer composites. A 45° trigger in the crash front was modeled by staggering the length of the two layers of elements. Both the drop mass and plug initiator were modeled with a rigid material model. The thickness of the tubes was 2.41mm.    Figure 4-5. Detailed FE crash model and the drop tower test setup at UBC [22].  4.3.1.2. Delamination model    63 The interface between the two layers of composites is modeled by a tiebreak contact: *CONTACT_AUTOMATIC_ONE_WAY_SURFACE_TO_SURFACE_TIEBREAK with contact option 8. There are multiple phenomena occur in the crash front: longitudinal bending, delamination between layers, splitting in transverse direction and debonding between fiber tows and matrix within each layer. In a homogeneous shell model, it was found that if the tiebreak strength is set at a low value, the delamination can easily spread to the whole tube and limit the development of the other damage modes. In this study, the tiebreak strength is raised to a much higher value than the results measured by mode I and model II tests in order to compensate for such defect of a macro-scale homogeneous FE model. A detailed sensitivity study on the tiebreak strength is performed in 7.1.3.  4.3.1.3. Boundary conditions and contact definitions  The top of the composite tube is attached to the drop mass by using *CONSTRAINED_EXTRA_NODES_SET. The drop mass is only movable in the vertical direction while the plug initiator is constrained in all directions. The interaction between different components are modeled by *CONTACT_AUTOMATIC _SINGLE_SURFACE.   There are extensive contact between the inner surface of the composite tube and the plug initiator. Thus a special experiment setup is built to measure the friction coefficient between them. The experiment setup is presented in Figure 4-6. A rectangular steel plate   64 is used in place of the plug initiator. It was in direct contact with the inner surface of a which is directly related to the friction coefficient by the following equation:                                                            (4-1)  The interior of the tube has a much rougher surface compared to the exterior due to the special molding method. The friction coefficient is measured as 0.34 between the tube and the plug initiator.    Figure 4-6. Experiment setup designed to measure the friction coefficient between the plug initiator and the inner surface of a composite tube  Gravity is considered in this study to account for the potential energy transformed during the crash event. The gravity acceleration is considered as a constant: 9.81m/s2. Both the drop mass and the composite tube are assigned with initial impact velocity measured right before impact. The velocity varies from 2.3m/s to 2.9m/s for 2-ply tubes and from 3.2 m/s to 5.0m/s for 4-ply tubes.   65  4.3.2. Crash simulation results of thin-walled composite tubes with plug initiator  In order to examine the influence of the irreversible strains and to investigate the capability of the ECDM model, two test cases were performed: the first case was to compare different versions of ECDM models by considering all irreversible strains (ECDM-ir), only pre-failure irreversible strains (ECDM-pre-ir) and no irreversible strains (ECDM-no-ir); the second case was to compare the ECDM model with MAT 58 in LS-DYNA. The stress/strain responses and energy responses of the crushing tube models were investigated.  4.3.2.1. Compare different ECDM models  Figure 4-7 compares the crash front morphologies predicted by different ECDM models with the one obtained by experiment. As shown, ECDM-no-ir and ECDM-pre-ir have very similar results, both of which predict relatively loose crash front fronds. Only ECDM-ir is able to predict the highly warped crash front fronds. This difference is more obvious at the beginning of crushing as shown in Figure 4-12. Both ECDM-no-ir and ECDM-pre-ir predict almost horizontal movement of the crash front fronds.  The major damage modes of a crushing composite tube with a plug initiator are longitudinal compression, transverse tearing, delamination and longitudinal bending. They are accompanied with micromechanical phenomena like fiber-matrix debonding, matrix cracking, fiber breakage and internal ply separation. A crushing tube firstly undergoes collision along the crashing direction. Then it initiates contact with the plug initiator, which   66 is followed by sliding, bending and transverse tearing. After leaving the transition area, the movement direction is changed from vertical to horizontal. The material undergoes unloading.   Figure 4-7. Comparison of the crash front morphologies predicted by different ECDM models with the one obtained by experiment.  It was found that the crash front morphologies after unloading are determined mainly by the amount of irreversible strains or residual strains.  As is shown in Figure 4-8, the outer part of the composite wall undergoes compression while inner par undergoes tension during bending. After leaving the transition area, only the model simulated with ECDM-ir has a large amount of longitudinal residual strains, which maintained the originally highly curved shape. The difference in simulations with different ECDM models is obvious, as is shown by the longitudinal strain responses of two adjacent elements from the outer surface and inner surface in Figure 4-9.   67  (a)  (b) Figure 4-8. Predicted failure morphology is highly affected by residual strains. (a) Prediction with ECDM-ir; (b) Prediction with ECDM-no-ir and ECDM-pre-ir.    68  (b) Figure 4-9. Longitudinal strain responses of two adjacent elements from the outer surface and inner surface for cases predicted with ECDM-ir and ECDM-pre-ir; positive and negative strains corresponding to tension and compression load.  Figure 4-10 compares the force-displacement and velocity responses of different ECDM models. Figure 4-11 lists some detailed comparisons. As shown, ECDM-no-ir and ECDM-pre-ir have very similar response in almost every aspect. ECDM-ir, however, is quite different.  It predicts higher values for the peak force, average plateau force and SEA, and lower values for the load ratio and crash length. According to the law of conservation of energy, the crash length is inversely proportional to the average plateau force. ECDM-ir predicts 10.4% higher average plateau force, 14.9% lower crash length and thus 10.9% higher SEA than ECDM-no-ir. This can be explained by Figure 4-12 and Figure 4-13.   69  (a)  (b) Figure 4-10. a) Predicted force-displacement and b) predicted velocity responses by different ECDM models.    70  (a)                                                                  (b)  (c)                                                                  (d) Figure 4-11. Detailed comparison of a) peak force, b) average plateau force, c) load ratio, d) crash length, e) SEA and f) average plateau force vs. crash length.         71 Figure 4-11   (e)                                                                  (f)  Figure 4-12 compares the predicted internal energy densities of the crash front elements. As shown, most of the predicted internal energy concentrates in the elements in the crash front and the areas close to the corners, where large transverse tearing and damage occurs. Generally, ECDM-ir predicts a higher level of internal energy density. Similar to Figure 4-9, Figure 4-13 traces the development of the longitudinal strain and internal energy density of the same element in the inner layer of a tube model. As shown in Fig. 4-13a, all the ECDM models predict the same level of stretching (peak value). The curves separate after unloading since only ECDM-ir predicts obvious residual strain which results in 10.6% higher internal energy density as presented in Figure 4-13b. It should be noted that the predicted energy difference from the crushing tube is different from the one obtained from the one-element model. This is because in a crushing tube model, the elements are undergoing different levels of loading-unloading.    72   Figure 4-12. Comparison of the predicted internal energy densities by ECDMs    (a) Figure 4-13. Predicted longitudinal strains and internal energy densities of the same element in the inner layer of a tube model by three ECDM models.     73 Figure 4-13 (contd)  (b)  In conclusion, considering the irreversible strains can improve the crash predictive ability of a composite constitutive material model. Both the pre-failure and post-failure irreversible strains are investigated. The results show that the post-failure irreversible strain is the dominant factor that results in higher internal energy stored in damaged materials and thus affects the predicted crash performance. (In the current model, pre-failure irreversible strain is relatively small and it is only considered in shear direction.)  4.3.2.2. Compare ECDM with MAT58  Figure 4-14 compares the crash front morphologies predicted by MAT58 and ECDM models with the one obtained by experiment. MAT58 predicts very soft crash front response that the fronds are loosely extended as is also presented in Figure 4-17. The curvatures of the fronds are much smaller than the one predicted by ECDM model and   74 the experimental result. In fact, MAT58 has even poorer prediction than ECDM-no-ir and ECDM-pre-ir. Nevertheless, MAT 58 and ECDM predicted very similar force-displacement responses as shown in Figure 4-15. Both models predicted a slightly higher peak force, load ratio, and crash length but a lower average plateau force than the experiment. The predicted SEAs are almost the same. There are several potential reasons for the mismatch: 1) the quality of the tube material is fair, with a moderate number of voids, resin rich areas and variable braid pattern [10]; 2) the thin shell element is not able to consider real through thickness deformation. The similar response of MAT 58 with ECDM is due to the fact that, as shown in Figure 4-3, the stress-strain curve described by MAT 58 is not as accurate which results in extra energy stored in the material.   Figure 4-14. Comparison of the crash front morphologies predicted by MAT58 and ECDM models with the one obtained by experiment.    75   Figure 4-15. Comparison of experimental results with the predicted force-displacement responses using the ECDM model and MAT58.     (a)                                                              (b) Figure 4-16. Detailed comparison of a) peak force, b) average plateau force, c) load ratio, d) crash length and e) SEA.     76 Figure 4-16 (contd)  (c)                                                              (d)  (e)  Figure 4-17 compares the predicted internal energy densities by MAT58 and ECDM. MAT58 predicts relatively uniform distribution of internal energy density while ECDM predicts higher concentrated internal energy density close to the corners and the crash front. Two elements are selected for a detailed comparison: one in the center (E1), one at the edge (E2). Figure 4-18 records the development of longitudinal strains and the   77 variation of internal energy density in these elements. As shown, the predicted evolutions of longitudinal strains are almost the same for the two elements by both models. However, the predicted internal energy densities are quite different due to different contributions of the transverse/shear stresses and strains. MAT58 predicts higher internal energy densities in the center and lower values at the edge than ECDM.    (a) MAT58                                    (b) ECDM Figure 4-17. Comparison of the predicted internal energy densities by MAT58 and ECDM model.     78  (a)  (b) Figure 4-18. Predicted longitudinal strains and internal energy densities of the same elements in the inner layer of a tube model by MAT58 and ECDM.  In conclusion, this test case demonstrated that different material models can be calibrated to match the experimental results in terms of force-displacement responses. However, the detailed local stress/strain and energy responses may not be accurate. Correctly   79 recovering real physical phenomena is vital for the improvement of a material model. The ECDM model is able to predict not only a good crash front morphology and force-displacement response but also more accurate local stress/strain and energy responses. Recover physically existing phenomena should be the only way to improve the accuracy of any prediction.         80 CHAPTER 5:  SHELL-BEAM ELEMENT METHOD  This chapter firstly examines some of the modeling difficulties of composite tubes without plug initiator, then performs detailed investigation on different element formulations in LS-DYNA and finally proposes and evaluates a new modeling approach, the shell-beam element method.  5.1. Simulation Difficulties of Composite Tubes without Plug Initiator   As is discussed in 1.2.2.2, a number of studies have attempted to model composite material tubes under axial impact [26,34,47,49,136-138]. The major composite models were developed for shell element due to its efficiency and practical potential to be used on large structures.   Park et al. [1] summarized three most popular approaches to model composites by using shell elements in LS-DYNA in the macro-scale level. 1) Single-layer approach  It was used by Bisagni et al. (2005) [31] and Deleo et al. (2010) [27]. In this approach, the composite is represented by one layer of shell elements. It is computationally efficient. However, this approach is not very stable, unable to simulate delamination and splaying of textile composites and to capture the local damage deformation.  2) Multi-layer approach   81 It was used by Mamalis et al. (2006) [26], Zarei et al. (2008) [32], Huang et al. (2009) [33] with MAT54, and Xiao (2009) [34] with MAT58 and CODAM. In this approach, the composite is represented by multi-layer of shell elements. It is able to simulate the delamination and splaying of textile composites by using certain contact algorithms. However, the simulated failure morphology was still not satisfactory for composite tubes without plug initiator. 3) Unit cell approach It was used by Tabiei and Chen (2001) [42], Cheng et al. (2008) [35], Li et al. (2009) [36] with MAT54, and Littell et al. (2008) [37], Littell et al. (2009) [38], Blinzler et al. (2010) [39], Xiao et al. (2011) [40], Blinzler et al. (2012) [41] with MAT58. The unite cell is a unit of repeated fiber architecture. It can be divided into several sub-cells, with each sub-cell consisting of fiber tows with varying size and orientation based on actual geometrical shape and location. Each sub-cell can be represented by an integration point through the thickness of the shell element. This approach is able to simulate more details of composites. However, it has the same instability issue as all the other models using shell elements.  In summary, very few models were able to capture both the failure morphologies and force-displacement responses of the crushing composite tubes without plug initiator. Crash simulations using shell elements have the tendency of instability, particularly when it comes to axial crash of composite tubes without plug initiators.  5.2. Investigation of Different Element Formulations in LS-DYNA    82 In order to find out the problem of existing modeling approaches, similar models with different types of elements are investigated. It should be noticed that even though the symmetric boundary conditions applied in non-full size models (half model or quarter model) can improve the stability of the simulation model, it is impossible to use a non-full size model on a complex structure, like a commercial vehicle. Thus only full size models are investigated in order to find a practical solution.   5.2.1. Different element formulations in LS-DYNA  As discussed in Chapter 5.1, traditional models using single-layer and multi-layer shell elements are not able to have a stable and consistent crash performance. This is mainly because the shell element is not well defined for in-plane compression. A little perturbation of any node could result in quite different bending behavior. As a result, different failure morphologies (bending inward, bending outward or buckling and catastrophic failure of the whole model) may be formed as shown in Figure 5-1. The problem of shell element is that it is a two dimensional element without real through thickness definition. Shell element is free to bend to either side without constraint. Conquering this deficiency can potentially solve the instability problem.     83    Figure 5-1. Unstable shell element model could predict different failure morphologies under the same loading conditions.  The existing element types with real through thickness definitions are solid element and thick shell element. In this research, only thick shell elements are investigated due to two reasons: firstly, if solid elements were used, at least 4 layers of elements are needed in order to correctly simulate the transverse shear and bending deformation. This is computationally expensive considering the fact that the in-plane dimension of a single element is 5mm×5mm while the thickness is only 1.2mm/layer; secondly, there is one type of thick shell element (thick-thick shell) which behaves just like solid element.    84  (a)                                               (b)                                          (c) Figure 5-2. Element discretization and nodal degree of freedom in local coordinate system for (a) fully integrated shell element; (b) thin-thick shell element; (c) thick-thick shell element. [139]  As shown in Figure 5-2, unlike the shell element, the thick shell element has a solid discretization. It can be classified into two major categories: the thin-thick shell and the thick-thick shell. The thin-thick shell element used a plane stress based constitutive law. -thick shell element used a 3D based constitutive law. Its thickness change comes naturally from the corresponding degree of freedom. So the thin-thick shell element behaves as a thickness enhanced shell element while the thick-thick shell element is like a solid element. In order to capture bending, the element stacking is not needed for thin-thick shell element while a minimum of 2 layers elements are needed for thick-thick shell element [139]. A comparison between these two categories of elements is provided in Table 5-1.    85 Table 5-1. Comparison between thin-thick shell and thick-thick shell elements. Element Type Constitutive Law Thickness Change Algorithm Element Layers to Capture Bending Element Formulations Thin-thick shell 2D plane stress based ratio 1 layer ELFORM=1,2 Thick-thick shell 3D based Corresponding degree of freedom  ELFORM=3,5  It should be noticed that since the thick-thick shell element requires 3D based constitutive law, a new material model is needed. For simplicity, a linear elastic behavior is assumed in the added out-of-plane directions without damage and failure.  5.2.2. Three test cases  In order to examine the capabilities of different types of thick shell elements, this section compares them with shell elements by three groups of test cases. All the models are built for a small size composite tube crushing into a rigid wall. In test case I, the tubes are assumed to be made of only 1 layer of triaxial braided composites without trigger in the crash front. In test case II, the tubes are the same as the ones in test case I but with trigger. In test case III, the tubes are assumed to be composed of 2 layers of triaxial braided composites with trigger. The element types selected for these test cases are fully integrated shell element, thin-thick shell element with ELFORM 2 and thick-thick shell element with ELFORM 5 in LS-DYNA.  Since the tests are designed to examine the capabilities of the element types, the classic material models instead of ECDM are used to avoid unnecessary instabilities. The   86 material models are based on the widely used Matzenmiller-Lubliner-Taylor (MLT) model [121]. Both 2D-MLT material model and 3D-MLT material model are implemented in LS-DYNA for thin-thick shell element and thick-thick shell element separately.   5.2.2.1. Test I: One layer of triaxial braided composites without trigger  In this test case, one layer of triaxial braided composites without trigger in the crash front is simulated by the selected three different types of elements. In order to correctly capture bending of the tube walls, 2 layers of continuous thick shell elements are used for thick-thick shell model instead of 1 layer of elements like the other two models. As shown in Figure 5-3, both two thick shell models predict severe buckling in the middle of the composite tubes. This is because thick shell elements have real thickness, which results in face contact in the crash front. A larger amount of instability is generated in the middle rather than at the crash front, thus severe buckling and failure occur in this region. The shell element model as shown in Figure 5-4 also predicts buckling in the middle of a tube at the beginning. However, since it is a line contact in the crash front, the shell elements can easily slide to the sides (not necessarily inside). The weakest part of the whole tube quickly concentrates in the crash front which results in stable crash afterwards.     87    (a)                                          (b)                                         (c) Figure 5-3. Comparison of the simulation results of 1 layer of triaxial braided composite tubes without trigger by (a) fully integrated shell elements; (b) thin-thick shell element with ELFORM 2; (c) thick-thick shell element with ELFORM 5.     Figure 5-4. Progressive crash of a l-layer composite tube made by fully integrated shell elements without trigger in the crash front.  From this test case, it can be concluded that: 1) The thick shell element is more stable than shell element due to real through thickness definition; Thick shell element models are able to predict the severe buckling responses of crushing tubes without trigger.   88 2) The defect of shell element can transfer some numerical instability from other regions of a model to the crash front region; Shell element model is not able to predict stable buckling response in this case. 3) The models with thin-thick shell elements and thick-thick shell elements predict similar results.  5.2.2.2. Test II: One layer of triaxial braided composites with trigger  In real tests, triggering is a process to form stress concentration on the edges of profile geometry to originate localized failure [17]. Triggers on composites tubes can be used to prevent catastrophic crushing, like severe buckling in test case I. In this test case, the same simulation models are used with triggers in the crash front. Detailed schematic of the triggers are shown in Figure 5-5. As presented, the thick shell element models are able to directly model the 45° triggers while the shell element model can only define a thinner shell element in the crash front to approximate the trigger.       89  (a)                                          (b)                                         (c) Figure 5-5. Schematic of the triggers used in the composite tubes in test case II modeled by (a) fully integrated shell elements; (b) thin-thick shell element with ELFORM 2; (c) thick-thick shell element with ELFORM 5.      (a)                                          (b)                                         (c) Figure 5-6. Comparison of the simulation results of 1 layer of triaxial braided composite tubes with trigger by (a) fully integrated shell elements; (b) thin-thick shell element with ELFORM 2; (c) thick-thick shell element with ELFORM 5.    90                  (a)                    (b) Figure 5-7. Progressive crash of a l-layer composite tube made by (a) thin-thick shell elements and (b) shell elements with triggers in the crash front.  Figure 5-6 shows the simulation results of 1 layer of triaxial braided composite tubes with trigger by the selected three different types of elements. As presented, the two thick shell element models predict large buckling. The locations with the largest out-of-plane deformation are closer to the crash front comparing to the results in test case I. This can be explained by Figure 5-7. As shown in Figure 5-7a, for the thick shell element model, the trigger is firstly forced to bend inward and then be deleted due to large deformation before any mature damage mode is formed. For a shell element model as presented in   91 Figure 5-7b, the trigger also tends to bend inward initially and then is forced to fold against the rigid wall. The rotational degrees of freedom of the shell elements save the trigger from being deleted. Nevertheless, the trigger is not functional after folding. For all models tested in this case, the loads are initially concentrated at the crash front region but then transferred to the whole structure after losing the trigger. As a result, the triggers in this test case are able to delay the instability but not sufficient to lead to desired crashing morphology, like splaying. In conclusion, no significant improvement is made by introducing triggers on 1-layer composite tube models.  5.2.2.3. Test III: Two layers of triaxial braided composites with trigger  This test case is focused on performance of selected elements on multiple layers of composites. Ideally, separated frond-like failure morphologies are expected in the crash front. Trigger is also used in all the models. Detailed schematic of the triggers are shown in Figure 5-8. The trigger for shell element model is composed of 2 layers of staggered elements while the triggers for thick shell element models are directly modeled with 45°. CONTACT_AUTOMATIC_ONE_WAY_SURFACE_TO_SURFACE_TIEBREAK is used to simulate delamination between layers. It should be noticed that the thick-thick shell element model is composed of two continuous inner layers of elements and two continuous outer layers of elements. They are connected together by tiebreak contact on the middle surface.    92  (a)                                          (b)                                         (c) Figure 5-8. Schematic of the triggers used in the composite tubes in test case III modeled by (a) fully integrated shell elements; (b) thin-thick shell element with ELFORM 2; (c) thick-thick shell element with ELFORM 5.    (a)                                          (b)                                         (c) Figure 5-9. Comparison of the simulation results of 2 layers of triaxial braided composite tubes with trigger by (a) fully integrated shell elements; (b) thin-thick shell element with ELFORM 2; (c) thick-thick shell element with ELFORM 5.  Figure 5-9 shows the simulation results of 2 layers of triaxial braided composite tubes with trigger by the selected three different types of elements. As presented, both thick   93 shell element models have large buckling and humped crash front. For all three models, even with delamination, the inner layers of elements and the outer layers of elements behave in a synchronous fashion, that is, they bend toward inside and outside together. It should be noticed that the simulations of thick shell models tend to automatically terminate very early due to convergence issues like negative element volume and infinite small timesteps. No mature failure morphology can be formed before termination. The shell element model, in the other hand, does not have such issue.   5.2.2.4. Test case summary  After all, none of the three types of elements are able to predict the separated frond-like failure morphologies of a crushing tube with multiple layers of triaxial braided composites. The following conclusions can be drawn from these three test cases: 1) Both shell element model and thick shell element models are not able to correctly predict the crush performance of composite tubes. 2) Modeling real thickness (by using thick shell elements in these test cases) in a model can largely improve its stability. 3) Large deformation on elements with solid discretization may encounter convergence issues like negative element volume and infinite timesteps. This is especially true for thin-thick shell element models. 4) Shell element is unstable but very flexible. It does not have any convergence issues in these test cases. The shell element models are always able to conquer   94 instabilities at the beginning of a crash event and reach a stable state afterwards. They may lead to good prediction with proper guidance.  5.3. Proposed New Modeling Method  Based on the conclusion in 5.2.2.4, this section is dedicated to develop a new modeling technique to improve the stability of crush simulation of composite tubes.   5.3.1. Shell-Beam element method  In this study, a new type of element, called shell-beam (SB) element, is designed to have the flexibility and efficiency of 2D-discretized elements and the stability of 3D-discretized elements. As shown in Figure 5-10, a SB element is composed of 2 shell elements and 4 beam elements. They are connected at the nodes.  Modeling structures using shell elements reinforced with beam elements have been used in a number of applications, like to simulate a hang-glider or all kinds of buildings, in which cases the beam elements are used to model beam-type structures. In the proposed SB element, however, the beams serve as the through-thickness reinforcement for shell elements. One layer of SB elements corresponds to one ply of composite. In Figure 5-10, the top and bottom shell elements have the same properties. Each shell element has one half the thickness of one composite layer. For a regular SB element, the beam length is equal to the thickness of each shell element. In a single SB element, each beam element has a cross section area equal to a quarter of the area of the shell. In a structure, each   95 beam element is shared by four pairs of shell elements.  Its cross section area is equal to the area of a shell element.   (a)  (b) Figure 5-10. Schematic of the proposed shell-beam element  Since both the shell element and beam element have six degrees of freedom at each node, the SB element is very flexible to adapt to different deformations. As a result, the shell elements are able to function just like the traditional shell elements, while the beam elements take the through-thickness tension, compression and out-of-plane shear. The SB element is  Shell thickness Shell-beam thickness Shell element Beam element Shell element Shared node  2 shell elements  4 beam elements   96 expected to improve the stability of a simulation model subjected to in-plane compression. In essence, a SB element is classified as a 2½D element. This methodology is called the SB element method.  5.3.2. Evaluation of shell-beam element method  The material properties of shell elements in Shell-Beam elements are determined by the same way as for single layer shell element models. The material properties of the beam elements, ideally, should be measured in the through thickness direction. However, these tests are very difficult to perform. In this study, the beam elements are assumed to have the same elastic modulus as the shell elements in the transverse direction. An elastic-plastic piecewise linear material model is used with the beam elements in SB elements without failure, considering that the beam represents the bond within a composite ply which does not delaminate even after severe damage. The beam elements are defined by *SECTION_BEAM in LS-DYNA [140]. *MAT_SPOTWELD is used as the material model for beam elements.   As shown in Figure 5-10, the beam elements overlap with the shell elements. To maintain the mass of the structure, the material density for the SB element is adjusted:                                                  (5-1)    97 where  is the measured material density while  is the effective density of the shell-beam element; ,  and  represent the total volume, beam element volume and shell element volume. As shown in Figure 5-10, , . Substituting them into Eq. 5-1:                                                               (5-2)  Detailed material properties for beam elements are listed in Table 5-2. Figure 5-11 shows the stress-strain response of the beam elements.  Table 5-2. Material properties of beam elements   7.52 GPa   0.3 Initial yield stress  125 MPa Hardening modulus  1.0 GPa Density  0.88 kg/mm2 * Assumed to be equal to the transverse modulus of the composite.    98  Figure 5-11. Stress-strain response of beam elements  5.3.2.1. One element tests  In order to fully evaluate the SB elements, one element tests are conducted in which a single element is subjected to load along major material directions. As shown in Figure 5-12, there are twelve test cases, which include tension/compression in every major material direction and shear in every major plane. The red arrows indicate the loading direction, and the blue triangles represent the constraints.  Table 5-3 compares the input material properties and the corresponding simulation results with one SB element. It should be noted that the inputs are for either shell elements or beam elements while the simulation is for the whole SB element. Generally, the simulation results match the inputs very well. The SB element offers relatively low resistance to the out-of-plane shear deformation (in 13 and 23 planes). The two out-of-plane shear moduli are close and strongly related to the material properties of the beam   99 elements. They are merely affected by the in-plane material properties of shell elements. A weak out-of-plane shear response is a limitation of the current SB element method.                                          (a)                                   (b)                                 (c)                                        (d)                                   (e)                                 (f) Figure 5-12. One element tests on single shell-beam element. T, C, S refer to tension, compression and shear separately.        100 Table 5-3. Comparison of input moduli and simulated results by one shell-beam element  Input Simulation   Input Simulation E11t 64.7 65.2  G12 7.52 7.45 E11c* 64.7 64.6  G21 7.52 7.43 E22t 12.5 12.5  G13* - 1.73 E22c* 12.5 12.6  G31* - 1.75 E33t 12.5 12.5  G23* - 1.72 E33c* 12.5 12.6  G32* - 1.75 * Input compressive moduli are assumed to be same with input tensile moduli. **No experiment is done for out-of-plane shear. The out-of-plane shear properties are generated automatically by shell-beam elements without corresponding inputs.  5.3.2.2. Cantilever Beam Tests  The cantilever beam test is chosen to further evaluate the SB element since it includes major loading conditions like in-plane compression, tension and out-of-plane shear seen in axial crush of a tube structure. The SB element is compared with shell element and solid element. Three series of virtual cantilever beam tests are done with different length-to-depth ratios. Figure 5-13 shows the cantilever beam modeled by the three different types of elements with length-to-depth ratio as 8:1 as an example. Table 5-4 summarized the model geometries used in the three tests.    101   Figure 5-13. Cantilever beam finite element models with a length to depth ratio 8:1.  Table 5-4. Summary of specimen configurations. Test Length L (mm) Width W (mm) Thickness T (mm) Support Span-to-Depth Ratio 1 19.2 9.6 2.4 8:1 2 38.4 9.6 2.4 16:1 3 76.8 9.6 2.4 32:1  In order to compare the performance of the elements alone, an elastic model MAT_01 in LS-DYNA [44] is used since it is well-defined for both the shell element and the solid element. The Youngs modulus is assumed to be 65GPa for all the elements tested in this case. The virtual test is designed to evaluate the element response under complex loads, especially out-of-plane shear load. Thus, low length-to-depth ratios, 8:1 and 16:1 are used to amplify the corresponding effects on purpose.    102 For a cantilever beam with end load, the end deflection (D) and end load (P) have the following relationship [141]:                                                            (5-3)  where I is the moment of inertia of area. Equation 5-3 can be rewritten as:                                                (5-4)  where r is the length-to-depth ratio, r=L/h. Equation 5-4 provides the analytical solution for this test case.   Figure 5-14. Force-deflection responses predicted by different types of elements under three length-to-depth ratios. The results are filtered by SAE1000 in LS-PrePost.    103 Figure 5-14 compares the predicted force-deflection responses with analytical solutions. As shown, the predictions of the shell element model and solid element model match the analytical solution very well under all conditions. The prediction of the SB element model is very close with the analytical solution at high r ratios (16:1 and 32:1). Under the low r ratio condition (8:1), nevertheless, the SB element model predicts a little lower stiffness as a result of lower out-of-plane shear properties.   5.3.2.3. Plate Crash Tests  This virtual test case is designed to test the stability of three types of elements under dynamic in-plane compression. Two parameters, a boundary condition and a material property are slightly varied in three tests. These are the initial crash velocity and the friction coefficient between the plate and the right wall. As shown in Figure 5-15, each plate model is attached to a mass block and crushed into a rigid wall. Test1 and Test2 differ slightly in velocity (5m/s and 4.8m/s), while Test2 and Test3 differ slightly in the friction coefficient between the plate and the right wall. The results show that the plate modeled with shell elements buckles easily under in-plane compression load, and a small variation of any one of the two tested parameters can change the bending direction of the plate. The plate with SB elements, on the other hand, behaves like the plate modeled with solid elements that is crushed by in-plane compression in all three tests. It can be concluded that, under in-plane compression, the SB element is as stable as the solid element and much more stable than the shell element.    104  Figure 5-15. Plate crash tests by using three types of elements  In brief, the SB element is much stable than the shell element and as accurate as the other elements besides out-of-plane shear directions. This limitation does not affect the performance of SB element in simulations of thin-walled structures which are generally built with high length-to-thickness ratios. SB element method is valid to be used in this study.   105 CHAPTER 6: CRASH SIMULATION OF COMPOSITE TUBES WITH ECDM AND SHELL-BEAM ELEMENT METHOD  This chapter is dedicated to evaluate the ECDM model and shell-beam element method in axial crash simulations of composite tubes. In general, composite tubes with plug initiator are used to evaluate the ECDM model first due to their ability to isolate some instable factors. Then, the same calibrated material properties and modeling parameters are used in the simulations of composite tubes without plug initiator.  6.1. Crash Simulation Models  As mentioned in Chapter 4, different tube configurations are investigated in this study. Special attention is paid to 2 and 4 tubes first since they are tested both with and without plug initiators. Both the ECDM material model and shell-beam finite element modeling method will be used to predict their crushing performance.  The tube model with plug initiator is generally the same with the one depicted in 4.3.1.1 except that the shell elements are replaced by shell-beam elements. Figure 6-1 shows the detailed finite element crash model without plug initiator. It is composed of three major components: a drop mass, a composite tube made by shell-beam elements and a rigid wall. The drop mass is modeled with a rigid material model (*MAT_RIGID).    106   Figure 6-1. Detailed FE crash model without plug initiator and the drop tower test set-up. [22]    Figure 6-2. Schematic of the cross-section of a  2 tube model   Each ply of composites 2 tube has two composite plies and thus it is modeled with two layers of SB elements which contains 4 layers of shell elements and 2 layers of beam elements as shown in Figure 6-2. Similarly, 4 tube is modeled with 4 layers of SB elements. The interface between any two composite plies is modeled by a tiebreak contact: *CONTACT_AUTOMATIC _ONE_WAY_SURFACE_TO_SURFACE_TIEBREAK   107 with contact option 8. The interaction between different components is modeled by *CONTACT_AUTOMATIC _SINGLE_SURFACE.  Both the drop mass and composite tube are assigned with an initial impact velocity measured right before impact. The velocity varies from 2.3m/s to 2.9m/s for 2 tubes and from 3.2 m/s to 5.0m/s for 4 tubes. In simulations, initial velocities of 3.0m/s and 4.19m/s are assigned to 2-ply tube and 4-ply tube respectively.  In this study, all the simulations were done using single precision LS-DYNA ls971s R6.1.0 on WINDOWS 64bit operating system. The processor is an Intel Core 2 Quad CPU Q8200 @2.33GHz. Table 6-1 compares the computational time of 2 tube models with plug initiators using the shell element and the shell-beam element with just one CPU. Since the shell-beam element model introduces more shell elements and beam elements, the computational time is about 3 times as much as the shell element model.    Table 6-1. Comparison of computational times of 2 tube models with plug initiators using the shell element and the shell-beam element. Model CPU Number of Elements Computational Time Shell element model 1 5076 Shell elements 15h36m Shell-beam element model 1 10261 Shell elements 5184 Beam elements 48h31m     108 6.2. Crash Simulation of 2×22 Tubes and 2×24 Tubes 6.2.1. Damage/failure morphology  In this study, simulations on composite tubes with plug initiators are performed first due to their ability to isolate unstable factors such as oscillation. These models are ideal for testing the ECDM material model and calibration of some FE modeling parameters like the effective strength of the tiebreak contact. Composite tubes without plug initiators share the same material properties and modeling parameters with the models as the plug initiator.  As shown in Figure 6-3, the predicted failure morphologies match the experimental results very well for both tubes. Figure 6-4 shows that the progressive crash in the front of the tubes is very smooth with the help of plug initiators. All the tube walls are forced to bend outward.   (a) 2 Figure 6-3. Comparison of experimentally observed and predicted damage/failure morphologies of (a) 2 and (b) 4 tubes with plug initiator.    109 Figure 6-3 (contd)  (b) 4    Figure 6-4. Progressive evolution of predicted damage/failure morphologies in the crash front of 4 tube with plug initiator from cross-section views.    110 Figure 6-5 compares the failure morphologies predicted by the simulation with the tested composite tubes. As shown, the simulations are able to capture the experimentally observed ply-splitting failure mechanism and the frond-[0/±45]2 without plug initiator, the outer one layer bends outward while the inner one 4 without plug initiator, the outer two layers bend outward while the inner two layers bend inward.  The progressive evolution of the failure morphologies in the crash front is shown through cross-section views in Figure 6-6. The damage/failure process of a crushing tube with trigger is as follows: damage and breakage of the trigger, delamination, separation and splaying of inner and outer layers of composites, formation of debris wedge in the corners. The whole process is accompanied with contact friction. For 4 simulation model, delamination is also observed between inner two and outer two layers due to severe bending. For the case without plug initiator, no debris wedge is formed in the splitting initiation area. Without the help of a debris wedge, both models run smoothly without any instability problems.    111  (a) ]2  (b) 4 Figure 6-5. Comparison of experimentally observed and predicted damage/failure morphologies of (a) 2 tube and (b) 4 tube without plug initiator.    112  (a) 2  (b) 4 Figure 6-6. Progressive evolution of predicted damage/failure morphologies in the crash front of (a) 2 tube and (b) 4 tube without plug initiator from cross-section views.    113  6.2.2. Force-displacement and SEA responses  The predicted force-displacement responses for both 2 and 4 tubes crushed with and without plug initiators are presented in Figure 6-7 along with experimental results. As shown, most of the simulation results match the experimental results very well. The 4 models are more stable than the 2 models since a thicker wall experiences less radial oscillation during crushing than a thinner wall.   (a) 2 with plug Figure 6-7. Comparison of force-displacement responses for both 2-ply and 4-ply  tubes crushed with and without plug initiators.        114 Figure 6-7 (contd)  (b) 2 without plug  (c) 4 with plug       115 Figure 6-7 (contd)  (d) 4 without plug  Comparisons between simulation and experimental results for the peak force, average crush force and SEA are shown in Figure 6-8. The predicted initial peak forces match the testing results for the cases without a plug initiator, but are much higher for the cases with plug initiator. The mismatch may have been resulted from poor definition of contact and the rough modeling of the plug initiator. In experiments, the plug initiator was attached to the tube using hot glue. Such details are neglected in the current FE model. After the formation of stable crash modes, the predicted average plateau forces are very close to experimentally measured results. As a result, predicted load ratios for cases with plug initiators are much larger than the testing results. However, the predicted crush lengths and SEA values match the experimental results very well, because the energy absorption ability is mainly determined by the average plateau force.     116  (a) Peak force  (b) Average plateau force Figure 6-8. Detailed comparison between predictions and experiment results in terms of a) peak force, b) average plateau force, c) load ratio, d) crush length and e) SEA.     117 Figure 6-8 (contd)  (c) Load ratio  (d) Crush length     118 Figure 6-8 (contd)  (e) SEA  6.3. Crash Simulation of Some Other Triaxial Braided Composite Tubes without Plug Initiator  To examine the robustness of the SB element method, three other different tubes are simulated using the same material properties and parameters as in the previous investigations. These are a symmetric square tube with larger edge length 4×42, a round tube 2 and an asymmetric tube 2×42. All tests are performed without a plug initiator (thus less constraint).   The predicted damage/failure morphologies of these three tubes compare well with the tested results as presented in Figure 6-9. For 2 tube with round cross section, the axial, biaxial tows and matrix separated from each other in crush test and thus the tested   119 tube has a net-liked morphology. It is impossible to predict such morphology using FE models with homogenized properties. Nevertheless, the SB element method captured the general deformation trend. For the asymmetric 2×42 tube, the number of plies bending inward initially was equal to the number of plies bending outward in simulation. However, after a while, some parts of the inner layers were forced to bend outward by accumulated composite materials inside the tube. The final morphology obtained by simulation matches closely with the experimental result.                       (a)                                 (b)                                 (c) Figure 6-9. Comparison of damage/failure morphologies of (a) 4×42, (b) round [0/±45]2 and (c) 2×42 tubes without plug initiator   120  Detailed force-displacement responses are plotted in Figure 6-10. All simulations match the experimental results very well in terms of the average plateau force. The higher predicted peak force for the round tube may be due to different crush initiation mechanisms and the variation in tube quality and hence material properties. The round tubes had a lower fiber volume fraction than other shapes [40, 136]. As will be shown in Chapter 7, different modeling approach of trigger may have a significant effect on the predicted peak force. The square and rectangular tubes usually have four pairs of fronds whereas the round tubes form more fronds.  The SB element model predicted this trend without artificial help such as seeding initial cracks around the tube edge in some earlier studies.   (a) 4×42 Figure 6-10. Comparison of force-displacement responses for (a) 4×42, (b) round 2 and (c) 2×42 tubes without plug initiator.   121 Figure 6-10 (contd)  (b) round 2  (c) 2×42     122 In summary, the FE simulations using SB element method and ECDM model are able to correctly capture most of the damage modes and force-displacement responses of a composite tube with and without a plug initiator in a crash test. The SB element is effective in terms of conquering the difficulty of instability under in-plane compression observed in the models with shell elements and the convergence issues in models with the solid elements, like negative element volume and reduced timesteps. This method is also sufficiently robust to be used in different tubes with round or asymmetric cross sections.      123 CHAPTER 7: INVESTIGATION ON THE SIMULATION MODELS AND DETAILED RESULTS  The ECDM model and shell-beam element method have been demonstrated to be relatively accurate and robust in the simulations of axial crash of triaxial braided composite tubes. As shown in previous chapters, there are many input parameters in the simulations. To have a better understanding of ECDM model and shell-beam method, a sensitivity study is carried out for key model parameters. This includes the parameters which have been observed to have a significant impact on the crash prediction (trigger, tiebreak contact definition), the ones which are not measurable by experiment (eratios, element deletion strains), and the ones which may vary with time/speed or space (friction coefficient, mechanical properties).   The sensitivity of axial crash simulation of a composite tube to each one of the studied parameters is analyzed separately. In each study, only one parameter is varied at a time. Both damage/failure morphologies and force-displacement responses are used in the evaluation. All the virtual tests were done on 2×22 tube models. Both conditions with and without plug initiators were considered.  7.1. Sensitivity Studies  7.1.1. Trigger  As mentioned before, a trigger mechanism is introduced in the front of a crushing tube to avoid catastrophic failure and thus encourage smooth crash. Triggers are generally in the   124 form of chamfer edge. It is modeled by staggering the length of the front layers of elements. Even though only 45º chamfer was used in experiment, this virtual test studied three different chamfer angles (30º, 45º, 60º) and three different modeling strategies.   As shown in Figure 7-1, five different triggers are investigated. A two-ply 2×22 tube is modeled by 2 layers of shell-beam elements which are composed of 4 layers of shell elements, marked by different colors. It should be noted that since the shell elements used in this study are Belytschko-Tsay shell elements, tapered (nonuniform) thickness is not supported. It is impossible to build a perfect, sharp trigger.  Figure 7-1 also compares the predicted failure morphologies with different triggers under the conditions both with and without a plug initiator. It is shown that the simulation results of models with trigger 1, trigger2 and trigger3 are similar and all agree with the experiment observations very well. The trigger angle has little effect on the failure morphology as long as the trigger elements are staggered uniformly. On the other hand, trigger 4 and trigger 5 caused catastrophic failures in the models without the plug initiator. Without a staggered edge, the inner layer and outer layer did not separate. For all models with plug initiator, the trigger does not affect the formation of failure morphologies since the inward bending of the composite layer is restricted by the plug initiator.      125 Trigger Type With Plug Initiator Without Plug Initiator  Trigger1      Trigger2      Trigger3      Trigger4      Figure 7-1. Schematic of 5 different triggers and comparison of the predicted failure morphologies of 2×22 tube models under conditions both with and without a plug initiator.   126 Figure 7-1 (contd)  Trigger5       Figure 7-2a compares the force-displacement responses of a 2×22 tube model with a plug initiator for five different triggers. The models with trigger 1~3 have very similar force responses, independent of the trigger angle. The model with trigger 4 has a peak force 13% smaller than the others. This is most likely a contact issue as such peak force drop is not observed in the curves without plug initiator. The force-displacement responses of the models without the plug initiator is shown in Figure 7-2b. The simulations with trigger 4 and trigger 5 have a much lower load carrying ability after reaching the peak values. This is because these two triggers introduced local and global buckling into the structure instead of stable desired damage/failure modes.    127  (a)  (b) Figure 7-2. Force-displacement responses of 2×22 tube models with different triggers (a) with plug initiator; (b) without plug initiator.    128 In summary, trigger modeling strategy can affect the crash performance of a tube structure to a large extent in terms of both failure morphologies and force responses. Modeling trigger elements with equal stagger is the best strategy to reproduce experimental observations. The trigger angle, however, has little effect on the overall performance. Trigger 2 is generally used in this research.   7.1.2. Eratios  Eratio defines the amount of strain that is reversible. In ECDM model, there are three eratios in three material directions. In this study, different combinations of eratios were tested. It should be noted that since the shear strain is relatively small, only the longitudinal and transverse eratios were investigated. Eratio12 was fixed as 0.5. Firstly, a study was done by varying eratio1 and eratio2 separately on a crash tube with plug initiator. Figure 7-3 shows eratio2 merely affect the results, while eratio1 is the dominant parameter to determine the crash front morphology since the majority of damage and bending occur in this direction. A smaller eratio1 results in a larger irreversible strain and a higher residual stiffness and thus a larger bending curvature.    129  (a)                               (b)                            (c)                              (d) Figure 7-3. Different combinations of eratios result in different crash front failure morphologies.   Secondly, the dominant eratio, eratio1 is studied. In this test case, eratio12 was fixed as 0.5 and eratio2 was fixed as 0.75. Crash models both with and without the plug initiator were used. As shown in Figure 7-4, for the crushing tubes with plug initiators, eratio1 only changes the bending curvature of the crash front material. For the crushing tubes without plug initiators, however, the damage modes are also affected by eratio1. As in Case I, a large bending curvature forces both the inner and outer wall bend outward.          130  Case I: 0.5-0.75-0.5 Case II: 0.75-0.75-0.5 Case III: 0.95-0.75-0.5    With plug        Without plug     Figure 7-4. Comparison of different failure morphologies of 2×22 tube models with different eratio combinations. 0.5-0.75-0.5 corresponds to the case with eratio1=0.5, eratio2=0.75 and eratio12=0.5, and so on.  Figure 7-5 compares the force-displacement responses of these test cases. For the simulation with plug initiator, the average plateau force slightly increases with decreasing eratio1 due to a larger amount of residual energy, while peak force is merely affected. For the simulation without plug initiator, the peak force increases with decreasing eratio1. The average plateau forces vary without clear trend due to varied damage modes.   131  Figure 7-5. Force-displacement responses of 2×22 tube models with different eratio combinations.  In summary, eratios determines the amount of strain that is reversible under unloading conditions. They largely affect the residual status of the crash front materials and thus the bending curvatures. Eratio1 is the dominant parameter to determine the crash front morphology since the majority of damage and bending occur in this direction.   7.1.3. Tiebreak Contact Strength  *CONTACT_AUTOMATIC_ONE_WAY_SURFACE_TO_SURFACE_TIEBREAK with contact option 8 was used to model the delamination between composite layers in this study. The failure criterion of this contact is:                                                   (7-1)   132  where NFLS refers to the normal failure stress while SFLS refers to the shear failure stress. Theoretically, the normal failure stress and shear failure stress can be obtained by correlating simulations with mode I and mode II experimental results. However, the values determined in such correlations tend to predict a softer interlaminar connection and premature delamination. In this study, much larger failure stress values were used. They were calibrated using models with plug initiator and then applied to all models.   This section focuses on evaluating different combinations of NFLS and SFLS. Both NFLS and SFLS are increased from Case I to Case II to Case III. As shown in Figure 7-6, larger NFLS and SFLS make it harder to separate the composite layers. Case III shows catastrophic failure on the model without plug. Figure 7-7 presents detailed force responses. For all the cases with progressive failure, both the peak force and average plateau force increase with tiebreak contact strength. Compare Case III to Case I on the models with plug initiator, the peak force and average plateau force increase 54% and 60% respectively. Thus, the calibration of tiebreak contact strength is vital to have a good prediction.         133  Case I: 0.015-0.025 Case II: 0.035-0.050 Case III: 0.055-0.075    With plug       Without plug    Figure 7-6. Comparison of different failure morphologies of simulation tubes with different tiebreak contact strength. 0.015-0.025 corresponds to the case with NFLS=0.015, SFLS=0.025, and so on.   Figure 7-7. Force displacement responses of 2×22 tube models with different tiebreak contact strength.   134  Meanwhile, a test case was done only on SFLS since it was observed to have the most significant effect on the tube crash performance. As shown in Figure 7-8, little variation of SFLS (from 0.06 to 0.055 to 0.05) changes the crash front shape and delamination condition dramatically. A higher SFLS results in less delamination and thus smaller curvature of the crash front fronds.   Figure 7-8. Crash front morphologies and delamination conditions under different combinations of tiebreak strength parameters. In LS-DYNA, the contact gap fringe gives delamination fringe valued from 0-1. 0 refers to intact area while 1 refers to fully delaminated area.  In summary, the tiebreak contact strength has a significant effect on the crash performance. It changes the delamination conditions and thus the failure morphologies and force responses. The calibration of tiebreak contact strength is vital to have a good prediction. Among the two failure stresses, SFLS is the dominant factor.   135  7.1.4 Friction coefficient  The friction coefficient is a measure of the resistance of one object to the other moving object at their contact surfaces. In the test case studied in this section, friction occurs between 1) the composite tube and the plug initiator & rigid wall; 2) the bended composite material and the tube wall. For simplicity, this study assumed all the friction coefficients are the same in all the contact definitions in one model, varying from 0.2 to 0.4.  Figure 7-9 compares the predicted results of 2×22 tube models with different friction coefficients. As shown, the general damage modes and predicted failure morphologies are almost the same for both simulations with and without plug initiator. Only minor variation on the bending curvature of crash front materials is observed. However, as is plotted in Figure 7-10, the force responses are quite different. Both the peak forces and average plateau forces increase with the friction coefficient. Compare the case with f=0.4 to the case with f=0.2, with and without plug, the peak force and average plateau force both increase about 6KN, independent of the loading conditions.       136  f=0.2 f=0.3 f=0.4     With plug        Without plug    Figure 7-9. Comparison of different failure morphologies of 2×22 tube models with different friction coefficient.   Figure 7-10. Force displacement responses of 2×22 tube models with different friction coefficient.   137  In summary, if the friction coefficients stay in a reasonable range, like 0.2~0.4, the general damage modes and predicted failure morphologies of 2×22 tube models are the same. However, larger friction coefficients between different contact components result in larger force responses and thus more energy is expected to be dissipated per unit length of the tube.   7.1.5 Mechanical Properties and Residual Strength  In reality, the measured mechanical properties may have a certain variation. This is especially true for composite materials, like triaxial braided composite studied in this research. The fiber tows in the braid may be affected by the dynamic resin flow during moulding. Thus, voids, resin rich area and varied braid pattern can be found in all the tested panels and tubes. It is necessary to investigate these variations on the overall crash performance.  Two separate investigations were performed on mechanical properties. The first virtual test focuses on the properties have apparent impact on the pre-failure behaviors. All the moduli and initial failure strengths are changed by ±10%. The second virtual test focuses on the properties have more effect on the post-failure behaviors. All the residual strengths are varied by ±10%. Figure 7-11shows an example of the influence of these variations on the transverse stress-strain relationships.    138   (a) varied mechanical properties  (b) varied residual strengths Figure 7-11. Transverse stress-strain relationships with (a) varied mechanical properties (MP), (b) varied residual strengths.    139 The failure morphologies and force responses for the first test case are shown in Figure 7-12 and 7-13. Varying the mechanical properties by 10% has little effect on the crash performance of the models. The only exception is the case with -10% MP on a model without plug initiator. Softer material properties change the damage modes of one side wall and result in a little lower force response.    +10% MP 0% MP -10% MP    With Plug       Without Plug     Figure 7-12. Comparison of different failure morphologies of 2×22 tube models with different mechanical properties.    140  Figure 7-13. Force displacement responses of 2×22 tube models with different mechanical properties.  The failure morphologies and force responses for the second test case are shown in Figure 7-14 and 7-15. As shown, a higher residual strength changes the damage modes of one model without plug initiator. The other failure morphologies are not affected. Significant variation on the force responses is observed as all the forces increase with residual strength. For the cases without plug initiator, 10% increase of residual strength results in 6.2% increase in the overall force response while 10% decrease of residual strength results in 5.0% decrease in the overall force response.        141  +10% RS 0% RS -10% RS    With Plug       Without Plug     Figure 7-14. Comparison of different failure morphologies of 2×22 tube models with different residual strength.   Figure 7-15. Force displacement responses of 2×22 tube models with different residual strength.   142  In summary, the current crash models are able to sustain minor changes (at least ±10%) of mechanical properties. But the variation of residual strength may bring in significant change on the overall force response. This can be explained by the fact that post-failure region is generally much larger than the pre-failure region.  7.1.6 Element deletion strains  Element deletion strains are used to define when an element should be deleted. Element deletion is important in terms of simulating material fracture and crack propagation. It should be noted that since the elements can still take some load before final deletion, element deletion strains may largely affect the predicted load carrying ability of a crushing structure. There are two element deletion criteria in ECDM model, one is based on an effective strain, and the other is based on the strains in each major material direction. Only the second element deletion criterion is studied in this chapter. Since shear strain is found to be small for the test cases, only longitudinal and transverse deletion strains are studied with five different combinations.  Figure 7-16 compares the failure morphologies of 2×22 tube models with different element deletion strains. Case I, III and IV have the same transverse element deletion strain. Smaller longitudinal element deletion strain (Del1) in these test cases results in more deleted elements. When Del1 is lower enough like in Case IV, the elements are deleted quickly without the formation of stable fond-like failure   143 morphology. Case I, II and V have the same longitudinal element deletion strain. Varied transverse element deletion strain has little effect on the deformation shape. These observations can be explained by the fact that: longitudinal crushing occurs on every element while transverse tearing only has major effect on the corner elements which will be deleted regardless of the element deletion strains due to large deformation.    With Plug Without Plug    Case I: Del1: 0.35 Del2: 0.5      Case II: Del1: 0.35 Del2: 0.35    Figure 7-16. Comparison of different failure morphologies of 2×22 tube models with different element deletion strains.          144 Figure 7-16 (contd)     Case III: Del1: 0.50 Del2: 0.50      Case IV: Del1: 0.20 Del2: 0.50      Case V: Del1: 0.35 Del2: 0.75    Figure 7-17 presents the corresponding force responses. By comparing Case I, III and IV, it is found that the longitudinal element deletion strains have little effect on the force responses besides the extreme case IV. A comparison of Case I, II and V shows that the forces increase with the transverse element deletion strains for both cases with and   145 without plug initiator. This is because most of the elements are deleted at the corners governed by the transverse element deletion strain.  (a) with plug initiator  (b) without plug initiator Figure 7-17. Force displacement responses of 2×22 tube models with different element deletion strains.    146 In summary, longitudinal crushing is the dominant damage/failure mode. A proper selection of the longitudinal element deletion strain is vital to the formation of fond-like failure morphology. However, if the value is large enough, it has little effect on the crushing force. On the other hand, the transverse element deletion strain does not influence the failure morphology, but has a notable effect on the force responses.  7.2. Investigation on detailed simulation results  In order to have a better understanding of the crash process inherently, this section examines some detailed simulation results based on all optimized results in 7.1. ECDM model traces many variables through history functions in LS-DYNA. A detailed history list can be found in the ECDM material code. This section only examines the following four states on 2×22 tube models both with and without a plug initiator: damage, residual stiffness, delamination and energy.   7.2.1. Damage conditions  The damage states are quantified by damage variables, range from 0 to1. 0 means undamaged while 1 means completely damaged. Figure 7-18 compares three major damage variables on 2×22 tube models both with and without the plug initiator. In order to have a better view, inner layer and outer layer are displayed separately. The following observations and conclusions can be drawn from this plot:   147 1) By comparing different damage variables, it is found that the dominant damage mode is in the longitudinal direction for most materials. The transverse damage and shear damage are concentrated around the corners. 2) By comparing inner layer with outer layer, the outer layer is more severely damaged than the inner layer for the model with plug initiator, while the opposite is true for the model without plug initiator. In the latter case, the inner layer is under severe bending and squeezing inside a highly constrained space. 3) By comparing the same layer of composites on two models, almost all the materials are damaged more severely on the model without plug initiator, especially for the materials in the middle of each wall. As a result, the materials are better utilized to contribute to the crash resistance. This explains why the same type of tube exhibits a significantly higher crushing force when crushed without a plug initiator as compared to when crushed with a plug initiator.   It should be noted that away from the crash front, the material may also be damaged to some extent as shown by the transverse damage state on both outer layers. However, these damage states are not large enough to break the material.         148  With plug  Without plug   Inner layer Outer layer  Inner layer Outer layer                         Figure 7-18. Comparison of three major damage states of 2×22 tube models both with and without plug initiator during crash.   7.2.2. Residual stiffness  The residual stiffness refers to the transient stiffness during crash. It is inversely proportional to the corresponding damage variable which is also governed by eratios. Figure 7-19 compares the residual stiffness of 2×22 tube models both with and without plug initiator during crash. All plots are normalized by their initial elastic moduli. Generally, more severe damage results in smaller elastic modulus. The residual stiffness can drop to 10% of the initial stiffness or even lower. Conclusions similar to that of 7.2.1 can be drawn based on these plots.    149  With plug Without plug   Inner layer Outer layer Inner layer Outer layer                         Figure 7-19. Comparison of the elastic moduli of 2×22 tube models both with and without plug initiator during crash. All the plots are normalized by their initial elastic moduli.  7.2.3. Delamination  Delamination is modeled by tiebreak contact in this study. Interface variables are used to trace the development of delamination: the resultant interface force and the contact gap. The resultant interface force is used to trace the interaction of the inner and outer layers. The contact gap is directly used to measure the delamination state, range from 0 to1, where 0 means the intact state while 1 means the fully delaminated state.  Figure 7-20 compares the delamination states of the 2×22 tube models both with and without plug initiator during crash. It should be noticed that the interface is attached   150 to the outer layer, thus it has a similar shape with the outer layer. Meanwhile, the interface related to deleted elements is still displayed in the interface plots (gray corner areas on the plots). It does not affect the results. From the resultant interface force response, it is found that interaction occurs between the layers all over the tube. However, it is not severe enough to cause delamination until the material reaching the crushing end of the tube. The contact gap plot indicates almost no delamination transition zone exists in both models. The delamination state switch from intact to fully delaminated right at the crushing end of the tube.     Resultant Interface Force Contact Gap     With plug       Without plug    Figure 7-20. Comparison of the delamination states of 2×22 tube models both with and without plug initiator.   151  7.2.4. Energy breakdown  Figure 7-21 displays the growth of different energy components. The total energy is a constant which is initially composed of the kinetic energy and potential energy. They are gradually transformed to the internal energy and dissipated by the sliding energy. The internal energy is mainly accumulated through material deformation, damage and fracture while the sliding energy is accumulated through friction and delamination since the tiebreak contact is defined to model delamination in this study. It should be noted that the two models were crashed with different initial velocities. Thus the initial kinetic energy and the total energy are slightly different.  For both models, the majority of the energy is absorbed through damage growth, which converted to the internal energy: 56% and 66% for the model with and without a plug initiator respectively. The remaining energy is dissipated by sliding energy. With a plug initiator, more energy is dissipated through friction. Without a plug initiator, more energy is absorbed by severe damage and deformation of the composite materials.   It should be noted that the delamination energy only composes 30% of the sliding energy. However, delamination is vital to determine the bending behavior of the crash front fronds and thus the force-displacement responses.    152  (a) 2×22 tube model with plug initiator  (b) 2×22 tube model without plug initiator Figure 7-21. Energy breakdown of 2×22 tube models (a) with and (b) without plug initiator.     153 CHAPTER 8:  CONCLUSION AND FUTURE WORK  8.1. Conclusion  In order to improve the predictive capability of crash simulation models for thin-walled composite structures subjected to axial crush loading, this work developed a new crash model. It is composed of two major components: an enhanced continuum damage mechanics (ECDM) material model and a new Shell-Beam (SB) Element modeling technique.   The ECDM model is composed of two sub-models, a pre-failure model and a post-failure model, to describe the pre-peak and post-peak behaviors separately. The pre-failure model uses a modified Ladevèze model while the post-failure model is composed of a strain-softening model and residual state model. An algorithm of calculating post-failure irreversible strains and two failure criteria (an initial failure criterion and a final failure criterion) are developed. The initial failure criterion determines the peak strength, while the final failure criterion determines the element deletion conditions.  The ECDM model was implemented as a user material model in LS-DYNA and evaluated for triaxially braided composites. The separated damage evolution laws in the pre- and post-peak regions allow the ECDM model to describe the stress-strain response of braided composites from a stress-free state to final failure in all three major material directions more accurately than the existing CDM models. This work shows that the post-  154 failure irreversible strain is a dominant factor affecting the predicted crash performance. In simulations of axial crash of tubes with a plug initiator, only the model with post-failure irreversible strain was able to simulate a crash frond morphology resembling that of experiment.  A SB element has been developed to improve the stability of FE models for crashworthiness simulations of thin-walled composite structures. By introducing beam elements in the through-thickness direction between two layers of shell elements, this method enriches the overall through thickness deformation capability of the element. The SB element was evaluated with one element, cantilever beam, plate in-plane compression, and tube crash tests. The results show that, along with the ECDM material model, the SB element method is able to correctly capture most of the damage modes and force-displacement responses of triaxial braided composite tubes with and without a plug initiator in a crash event. This method is effective in terms of conquering the difficulties of some traditional elements such as instability and non-convergence issue.   The sensitivities of composite crushing tube models on different parameters used by ECDM model and shell-beam method are investigated at the end of this study. It is found that: 1) Modeling the trigger region by staggering the shell layers with equal length difference is the best strategy. Varying the trigger angle from 30 to 60 degree has little effect on the overall performance. 2) Eratios, a set of parameters define the amount of strains that are reversible, determine the residual status of the crash front materials and thus the bending curvatures. Eratio1,   155 defines the reversible strains in the longitudinal direction, is the dominant parameter among them.  3) The tiebreak contact strength changes the delamination conditions and thus the failure morphologies and force responses. SFLS, the shear failure stress, is the dominant factor. 4) The friction coefficients, when their values are assigned within a reasonable range of 0.2~0.4, change the force responses, not the failure morphologies. 5) The current crash models are able to sustain a variation of ±10% in mechanical properties with little changes in crash performance. A ±10% variation in residual strengths, however, may bring in significant change on the overall force response.  6) The longitudinal element deletion strain is vital to form good, fond-like failure morphologies while the transverse element deletion strain has notable effect on the force responses.  8.2. Future Work  The works to be considered in future research are summarized as follows: 1) Investigate the mesh-size dependent issue; 2) Investigate the high initial peak force problem in the simulation models with plug initiators; 3) Introduce a better delamination model with parameters that can be quantified by experiments; 4) Add strain-rate sensitive parameters in the material model for strain-rate dependent composite materials;   156 5) Write SB element as a user defined element for more convenient applications; 6) Evaluate the crash model developed in this work under different loading conditions such as off-axis and lateral impact, and verify this approach in other structures; 7) Investigate the edge effect encountered in measuring the material properties with standard and modified test methods. A corresponding EDGE function can be introduced in ECDM model. A brief explanation of this work is provided below.  As was observed, the crash front material is weaker than the other part of a tube. The traditional straight-sided specimen, as shown in Figure 8-1a, can be used to measure the material properties of crash front materials since it has discontinuous fibers from one end to the other. On the other hand, the bowtie-shaped specimen shown in Figure 8-1b has continuous fibers connecting the two ends. It can be used to measure the material properties of the intact tube walls away from the crash front. Both Bowman et al. [143] and Park et al. [1] reported the use of bowtie-shaped specimen for 2D3A. It should be noted that in order to grip 100% of the axial and bias fibers in the gage section, the bowtie specimen has much shorter length than the straight-sided specimen.     157  (a)                                (b) Figure 8-1. Schematic of the amount of bias tows griped in (a) a straight-sided specimen and (b) a bowtie-shaped specimen [1]  The separation of material properties used for edge and intact areas could be very important for closed structures like composite tubes since they experience evolution of edge size and edge location during crash. Two different sets of tests should be done on straight-sided specimens and bowtie-shaped specimens separately.   In terms of modeling, a function SOFT is used in many composite material models in LS-DYNA [44]. It is used to reduce the strength applied on the crush front elements and it has been treated as one of the most influential parameters to calibrate in order to maintain a stable simulation [27]. However, SOFT is only designed to be adjusted artificially without any solid physical foundation. Considering the only way to calibrate this parameter is to compare the crash force-displacement responses obtained experimentally and numerically, it is very likely to be overused so that some of the other   158 important phenomena occurred at the same time are largely affected, like contact, delamination and any other parameter that is hard to be quantified by experiment.   In order to solve this problem, an EDGE function can be introduced in ECDM model to detect the position of an element is at the edge of a structure or not. The EDGE function can guide the edge elements which group of material properties (edge or intact material properties) to use instead of just reducing the material strength like in SOFT function.  The EDGE function works on a numerical model with element failure. As shown in Figure 8-2, once an element is completely failed, the corresponding failure flag (represented by the red numbers) will be set as 1. This failure flag is defined on the nodes which are shared with adjacent elements. The adjacent elements will be recognized as edge elements if any of their failure flags are detected as 1. For example, if element 2 is failed, element 4, 5, 6 will become edge element (element 1, 3 are already edge element). If element 5 failed, element 4, 6, 7, 8, 9 will become edge element.  It should be noticed that the initial edge elements (like element 1, 2, 3 in Fig. 13) can only be defined before simulation. They should be classified into one edge section in the finite element model with the edge properties. The rest of the elements can be automatically updated as edge element by EDGE function if any one of their adjacent element is failed.     159  Figure 8-2. Working schematic of EDGE function.    E1 E2 E3 E4 E5 E6 E7 E8 E9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 E1 E2 E3 E4 E5 E6 E7 E8 E9 0 1 1 0 0 1 1 0 0 0 0 0 0 0 0 0 After E2 failed   160           APPENDICES    161 APPENDIX A: ECDM HISTORY VARIABLES  The ply level ECDM model is written as user defined material models in a commercial explicit code LS-DYNA. In LS-DYNA, it is strain increment that passed into the user-defined material subroutines; thus all the equations used in ECDM model are written based on strain increment. ECDM model is written in Fortran 77 format in order to work properly in LS-DYNA. LSTC provides dyn21.f for the insertion of user defined material models. Each UMAT is written within a specified subroutine. In this work, subroutine umat47v is used. History variables are used to trace the development of interested parameters. 64 history variables are used in ECDM model. The Fortran code related to the definition of these variables is listed as follows:         subroutine umat47v(cm,d1,d2,d3,d4,d5,d6,sig1,sig2,      . sig3,sig4,sig5,sig6,eps,hist,lft,llt,dt1siz,capa,etype,tt) c****************************************************************** c****************************************************************** c|  livermore software technology corporation  (lstc)             | c|  ------------------------------------------------------------  | c|  copyright 1987,1988,1989 john o. hallquist, lstc              | c|  all rights reserved                                           | c****************************************************************** c c--03/03/2015--New Input Orders c       include "nlqparm" c-----SOFT------------------       include 'iounits.inc' C_TASKCOMMON (aux33loc) C_TASKCOMMON (maxsvloc) c      common/maxsvloc/mxsave,ipt_type(300) c      common/aux33loc/ c     & ix1(nlq),ix2(nlq),ix3(nlq),ix4(nlq),ixs(nlq,4),mxt(nlq)       common/aux33loc/ix1(nlq),ix2(nlq),ix3(nlq),ix4(nlq),ix5(nlq),      * ix6(nlq),ix7(nlq),ix8(nlq),mxt(nlq)       logical icrash,icrash_loc c      common/crash/icrash,iflagd,iflagm,iflagf,iflagc,jf c      common/bk00/numnp,numpc,numlp,neq,ndof,nlcur,numcl,numvc, c     1 ndtpts,nelmd,nmmat,numelh,numelb,numels,numelt,numdp, c     2 grvity,idirgv,nodspc,nspcor,numelh10,numels8 c       dimension ncrshf(1000000)       data ncrshf/1000000*0.0/   162       real qq1(nlq),soft c-----SOFT------------------ c c material property array and strain increments       dimension cm(*)       dimension d1(nlq),d2(nlq),d3(nlq),d4(nlq),d5(nlq),d6(nlq) c stresses and effective plastic strain (not used) returned by routine       dimension sig1(nlq),sig2(nlq),sig3(nlq),sig4(nlq),sig5(nlq),      & sig6(nlq),eps(nlq) c history variables       dimension hist(nlq,*) c step size       dimension dt1siz(nlq) c damage parameters       real w1(nlq),w2(nlq),ws(nlq) c damage thresholds       real r1(nlq),r2(nlq) c total strains         real e1tot(nlq),e2tot(nlq),e4tot(nlq),e5tot(nlq),e6tot(nlq) c original total strains       real e1toto(nlq),e2toto(nlq),e4toto(nlq) c constitutive tensor       real c(nlq),c11(nlq),c12(nlq),c22(nlq),c44(nlq),c55(nlq),c66(nlq) c failure flags c added in line 'common/failuloc/sieu(nlq),fail(nlq)'       logical failur,failgl       common/failuloc/sieu(nlq),fail(nlq)       common/failcm/failur,failgl,mtfail(1000)       common/failcmloc/ifail(nlq)       dimension failu(nlq),efail(nlq) c failure criteria and threshold variables       real const1(nlq),const2(nlq), g1(nlq),g2(nlq),gh1(nlq),gh2(nlq),      &  ge11(nlq),ge12(nlq),ge21(nlq),ge22(nlq),ge24(nlq) c element type flag       character*(*) etype c material constants       real bmod,y1t,y1com,y2t,y2com,n12,n21,n23,n13,g12,g23,g13       real y1(nlq),y2(nlq)       real s1t,s1c,s2t,s2c,ssc c damage variable function constants c added in dummys for ws equation       real factor, m1c,m1t,m2c,m2t,mss,      1     im1c,im1t,im2t,im2c,imss,      2     y1s1t,y1s1c,y2s2t,y2s2c,g12ssc       real dummy, dummys   163 c erosion variables       real erode,epsmax,wicrit1,wicrit2,wdcrit1,wdcrit2 c load and display on first call flag c added new variables for residual strength,  c           peak strain and peak damage values c added toteps variable for strain based erosion       logical once       data once /.false./       real res1t,res1c,res2t,res2c,ress      1     e1tpk,e1cpk,e2tpk,e2cpk,espk      2     w1tpk,w1cpk,w2tpk,w2cpk,wspk       real toteps(nlq) c--------------------------------------------------------- c decoupled model       real pm1t,pm1c,pm2t,pm2c,pms,check1(nlq),check2(nlq) c Ladeveze damage--------------       real y1o,y1c,y2c,y12c,y2o,y12o,y2oc,y2cc,b,gama1,y1oc,y1cc       real yd1(nlq),yd2(nlq),yd12(nlq) c new inputs       real input(16),exinp       data input/16*0.0/,exinp/0.0/ c Ladeveze plasticity in shear direction-------       real r0,beta,alfa       real f(nlq),rp(nlq),p(nlq),eratio(nlq),e4toteo(nlq),fps(nlq)       real de4e(nlq),de4p(nlq),e4tote(nlq),e4totp(nlq),dsig4(nlq),      &     fsig4(nlq),sig4o(nlq),drp(nlq),dp(nlq) c      real eratio(1000000) c--- initial failure point ----       real x11s(nlq),x22s(nlq),fc1(nlq),fc2(nlq),fc12(nlq),      &    w1if(nlq),w2if(nlq),wsif(nlq),e1if(nlq),e2if(nlq),e4if(nlq)       real fail1(nlq),fail2(nlq),fail12(nlq) c  Post-failure irreversible strain       real erat1(nlq),pirre1(nlq),e1tote(nlq),e1totp(nlq),de1e(nlq),      &     de1p(nlq),e1toteo(nlq),e1max(nlq)       real erat2(nlq),pirre2(nlq),e2tote(nlq),e2totp(nlq),de2e(nlq),      &     de2p(nlq),e2toteo(nlq),e2max(nlq)       real pirre4(nlq),e4max(nlq)       real ratio1,ratio2,ratio4 c -- erosion criteria -----       real e1dele,e2dele,e4dele,dele1(nlq),dele2(nlq) c . . . .   164 . . c c load history variables c       do 70 i=lft,llt         hist(i, 1)=   w1(i)         hist(i, 2)=   w2(i)         hist(i, 3)=   ws(i)         hist(i, 5)=   r1(i)         hist(i, 6)=   r2(i)         hist(i, 8)=e1tot(i)         hist(i, 9)=e2tot(i)         hist(i,11)=e4tot(i)         hist(i,12)=efail(i)         hist(i,13)=sig1(i)         hist(i,14)=sig2(i)         hist(i,15)=sig4(i)         hist(i,16)=eratio(i)         hist(i,17)=e4tote(i)         hist(i,18)=e4totp(i)         hist(i,19)=p(i)         hist(i,20)=f(i)         hist(i,21)=rp(i)         hist(i,22)=fps(i)         hist(i,23)=dp(i)         hist(i,24)=de4e(i)         hist(i,25)=de4p(i)         hist(i,26)=erat1(i)         hist(i,27)=e1tote(i)         hist(i,28)=e1totp(i)         hist(i,29)=pirre1(i)         hist(i,30)=g1(i)         hist(i,31)=gh1(i)         hist(i,32)=e1max(i)         hist(i,33)=erat2(i)         hist(i,34)=e2tote(i)         hist(i,35)=e2totp(i)         hist(i,36)=pirre2(i)         hist(i,37)=e2max(i)         hist(i,38)=g2(i)         hist(i,39)=gh2(i)         hist(i,40)=pirre4(i)         hist(i,41)=e4max(i)         hist(i,42)=y1(i)*(1-w1(i))   165         hist(i,43)=y2(i)*(1-w2(i))         hist(i,44)=g12*(1-ws(i))         hist(i,45)=toteps(i)         hist(i,46)=dele1(i)         hist(i,47)=dele2(i)         hist(i,48)=fail1(i)         hist(i,49)=fail2(i)         hist(i,50)=fail12(i)         hist(i,51)=w1if(i)         hist(i,52)=w2if(i)         hist(i,53)=wsif(i)         hist(i,54)=e1if(i)         hist(i,55)=e2if(i)         hist(i,56)=e4if(i)         hist(i,57)=e5tot(i)         hist(i,58)=e6tot(i)         hist(i,59)=sig5(i)         hist(i,60)=sig6(i)         hist(i,61)=ncrshf(ix1(i))         hist(i,62)=ncrshf(ix2(i))         hist(i,63)=ncrshf(ix3(i))         hist(i,64)=ncrshf(ix4(i))    70 continue       return   900 format (A30, e10.3, A5, e10.3, A5, e10.3, A1)       end c      166 APPENDIX B: MATERIAL CARDS FOR ECDM MODEL IN LS-DYNA   *MAT_USER_DEFINED_MATERIAL_MODELS is the keyword in LS-DYNA used to define all the input material properties for ECDM model. The detailed input cards are listed as follows:    The extra inputs are listed in order in a txt file. These inputs are defined in the code with the meanings as shown at the left hand side:   167       168                     BIBLIOGRAPHY     169 BIBLIOGRAPHY  1. Park, C-K., Kan, C-D., Hollowell, W., & Hill, S.I. (2012, December). Investigation of opportunities for lightweight vehicles using advanced plastics and composites. (Report No. DOT HS 811 692). Washington, DC: National Highway Traffic Safety Administration. 2. ASM Handbook Volume 21 Composites, ASM International, Material Park, OH, 2001. 3. Nageswara R. Janapala, Fu-Kuo Chang, Robert K. Goldberg, Gary D. Roberts and Karen E. Jackson, Crashworthiness of Composite Structures with Various Fiber Arichitectures, 11th international LS-DYNA Users Conference, 2010. 4. Libo Yan, Nawawi Chouw, Crashworthiness characteristics of flax fibre reinforced epoxy tubes for energy absorption application. 5. Mamalis, A.G; Robinson, M.; Manolakos, D.E.; Demosthenous, G.A.; Ioannidis, M.B.; Carruthers, J. (1997): Crashworthy capability of composite material structures. In Composite Structures 37 (2), pp. 109134. 6. Waimer, M.; Kohlgrüber, D.; Hachenberg, D.; Voggenreiter, H. (2013): Experimental study of CFRP components subjected to dynamic crash loads. In Composite Structures 105, pp. 288299. 7. Belingardi G, Chiandussi G. Vehicle crashworthiness design  general principles and potentialities of composite material structures. In: Abrate S, editor. Impact engineering of composite structures. Wien: Springer; 2011. p.193264. 8. Farley, G. L. (1983): Energy Absorption of Composite Materials. In Journal of Composite Materials 17 (3), pp. 267279. 9. Ramakrishna, S. (1997): Microstructural design of composite materials for crashworthy structural applications. In Materials & Design 18 (3), pp. 167173. 10.  11. carbon tubes using MAT58 in LS-hin-walled structures, Volume 47, Issues 6-7:740-749. 12. Thornton, P.H.; Edwards, P.J. (1982): Energy Absorption in Composite Tubes. In Journal of Composite Materials 16 (6), pp. 521545.   170 13. Hull, D. (1991): A unified approach to progressive crushing of fibre-reinforced composite tubes. In Composites Science and Technology 40 (4), pp. 377421. 14. J.M. Starbuck, S. Simunovic, G.C. Jacob, Test methodologies for determining energy absorbing mechanisms of automotive composite material systems, SAE Tech Pap Ser, 2000-01-1526 (2000), pp. 16. 15. Lee, H. K.; Simunovic, S. (2002): A damage mechanics model of crack-weakened, chopped fiber composites under impact loading. In Composites Part B: Engineering 33 (1), pp. 2534. 16. Jacob, George C.; Starbuck, J. Michael; Fellers, John F.; Simunovic, Srdan; Boeman, Raymond G. (2006): Crashworthiness of various random chopped carbon fiber reinforced epoxy composite materials and their strain rate dependence. In J. Appl. Polym. Sci. 101 (3), pp. 14771486. 17. Lau, Saijod T.W.; Said, M.R.; Yaakob, Mohd Yuhazri (2012): On the effect of geometrical designs and failure modes in composite axial crushing: A literature review. In Composite Structures 94 (3), pp. 803812. 18. gVolume 48, Issues 1-3:107-111. 19. Mamalis, A.G; Manolakos, D.E; Ioannidis, M.B; Papapostolou, D.P (2004): Crashworthy characteristics of axially statically compressed thin-walled square CFRP composite tubes: experimental. In Composite Structures 63 (3-4), pp. 347360. 20. Brimhall, Thomas Jay (2005): Friction energy absorption in fiber reinforced composites, PhD dissertation, Michigan State University, Publication Number: AAI3189616. 21. Feraboli, Paolo; Wade, Bonnie; Deleo, Francesco; Rassaian, Mostafa (2009): Crush energy absorption of composite channel section specimens. In Composites Part A: Applied Science and Manufacturing 40 (8), pp. 12481256. 22. inite element modeling of the Journal of Impact Engineering 37:662-672. 23. Botkin, Mark E.; Johnson, Nancy L.; Zywicz, Ed; Simunovic, Srdan (Eds.) (1998): Crashworthiness Simulation of Composite Automotive Structures. 13th Annual Engineering Society of Detroit Advanced Composites Technology Conference and Exposition. Detroit, MI, USA, Sep. 28-29. 24. McCarthy, M.A.; Wiggenraad, J.F.M. (2001): Numerical investigation of a crash test of a composite helicopter subfloor structure. In Composite Structures 51 (4), pp.   171 345359. 25. Schweizerhof, K.; Weimar, K.; Munz, Th.; Rottner, Th. (Eds.) (1998): Crashworthiness analysis with enhanced composite material models in LS-DYNA merits and limits. LS-DYNA World Conference. Detroit, MI, USA. 26. Mamalis, A.G.; Manolakos, D.E.; Ioannidis, M.B.; Papapostolou, D.P. (2006): The static and dynamic axial collapse of CFRP square tubes: Finite element modelling. In Composite Structures 74 (2), pp. 213225. 27. Deleo F, Wade B, Feraboli P, Rassaian M, Byar A, Higgins M (Ed.) (2010): Crashworthiness of Composite Structures Modeling of the Crushing of UD Tape Sinuisoidal Specimens using a Progressive failure model. AMTAS Fall Meeting. Edmonds, WA, Oct. 2010. 28. Kiani, Morteza; Shiozaki, Hirotaka; Motoyama, Keiichi (2013): Using Experimental Data to Improve Crash Modeling for Composite Materials. Composite Materials and Joining Technologies for Composites, Volume 7. New York, NY: Springer New York, pp. 215226. 29. Obradovic, Jovan; Boria, Simonetta; Belingardi, Giovanni (2012): Lightweight design and crash analysis of composite frontal impact energy absorbing structures. In Composite Structures 94 (2), pp. 423430. 30. Abrate, Serge; Castanié, Bruno; Rajapakse, Yapa D. S. (2013): Dynamic Failure of Composite and Sandwich Structures. Dordrecht: Springer Netherlands (192). 31. Bisagni, C., Pietro, G. D., Fraschni, L., & Terletti, D. (2005). Progressive crushing of fiber-reinforced composite structural components of a formula one racing car, Composite Structures, 68, 491-503.  32. Zarei, Hamidreza; Kröger, Matthias; Albertsen, Henrik (2008): An experimental and numerical crashworthiness investigation of thermoplastic composite crash boxes. In Composite Structures 85 (3), pp. 245257. 33. Huang, Jiancheng; Wang, Xinwei (2009): Numerical and experimental investigations on the axial crushing response of composite tubes. In Composite Structures 91 (2), pp. 222228.  34. Xinran Xiao (2009): Modeling Energy Absorption with a Damage Mechanics Based Composite Material Model. In Journal of Composite Materials 43 (5), pp. 427444.  35. Cheng, J., Binienda, W. K. (2008). Simplified braiding through integration points model for triaxially braided composites. Journal of Aerospace Engineering, 21(3):152-161. 36. Li, X., Binienda, W. K., & Littell,  J. D. (2009). Methodology for impact modeling of   172 triaxial braided composites using shell elements. Journal of Aerospace Engineering, 22(3):310-317.  37. Littell, Justin D.; Binienda, Wieslaw K.; Goldberg, Robert K.; Roberts, Gary D. (2008): A Modeling Technique and Representation of Failure in the Analysis of Triaxial Braided Carbon Fiber Composites. In NASA/TM (215245). 38. Littell, Justin D.; Binienda, Wieslaw K.; Arnold, William A.; Roberts, Gary D.; Goldberg, Robert K. (2009): Effect of microscopic damage events on static and ballistic impact strength of triaxial braid composites. In Composites Part A: Applied Science and Manufacturing 40 (12), pp. 18461862. 39. Blinzler, Brina J.; Goldberg, Robert K.; Binienda, Wieslaw K. (Eds.) (2010): Investigation of MAT_58 for Modeling Braided Composites. 11th Int'l LS-DYNA® Users Conf. Dearborn, MI, USA, June 6-8. 40. Xiao, Xinran; Kia, Hamid G.; Gong, Xiao-Jing (2011): Strength prediction of a triaxially braided composite. In Composites Part A: Applied Science and Manufacturing 42 (8), pp. 10001006. 41. Blinzler, Brina J.; Goldberg, Robert K.; Binienda, Wieslaw K. (2012): Macroscale Independently Homogenized Subcells for Modeling Braided Composites. In AIAA Journal 50 (9), pp. 18731884. DOI: 10.2514/1.J051342. 42. Tabiei, A. L.A.; Chen, Quing (2001): Micromechanics Based Composite Material Model for Crashworthiness Explicit Finite Element Simulation. In Journal of Thermoplastic Composite Materials 14 (4), pp. 264289. 43. Morgan Osborne, Single-Element Characterization of the LS-DYNA MAT54 Material Model, M.S. Thesis, University of Washington 2012. 44. LS-2007, Livermore Software Technology Corporation (LSTC). 45. DeTeresa, S.J., Allison, L.M., Cunningham, B.J., Freeman, D.C., Saculla, M.D., Sanchez, R.J. and Winchester, S.W. (2001): Experimental Results in Support of Simulating Progressive Crush in Carbon-Fiber Textile Composites, Lawrence Livermore National Laboratory, UCRL-ID-143287. 46. Xinran Xiao, Nancy Johnson, Mark Botkin, Challenges In Composite Tube Crash Simulation, Proceedings of the American Society for Composites 18th Technical Conference, paper 154, October 2003. 47. McGregor, Carla J.; Vaziri, Reza; Poursartip, Anoush; Xiao, Xinran (2007): Simulation of progressive damage development in braided composite tubes under axial compression. In Composites Part A: Applied Science and Manufacturing 38 (11), pp. 22472259.   173 48. Xiao X. (2010): A Coupled Damage-Plasticity Model for Energy Absorption in Composite. International Journal of Damage Mechanics, Vol. 19, pp. 727-751. 49. Xinran Xiao, Carla Mcgregor, Reza Vaziri, Anoush Poursartip, Progress in Braided Composite Tube Crush Simulation, Int J Impact Engineering, 36, 2009, 711-719. 50. Schuecker, C.; Pettermann, H. E. (2008): Fiber Reinforced Laminates: Progressive Damage Modeling Based on Failure Mechanisms. In Arch Computat Methods Eng 15 (2), pp. 163184. 51. Singh, C. V.; Talreja, R. (2013): A synergistic damage mechanics approach to mechanical response of composite laminates with ply cracks. In Journal of Composite Materials 47 (20-21), pp. 24752501. 52. Maimí, P.; Camanho, P.P.; Mayugo, J.A.; Dávila, C.G. (2007): A continuum damage model for composite laminates: Part I  Constitutive model. In Mechanics of Materials 39 (10), pp. 897908. 53. Tay, T. E.; Liu, G.; Tan, V.B.C.; Sun, X. S.; Pham, D. C. (2008): Progressive Failure Analysis of Composites. In Journal of Composite Materials 42 (18), pp. 19211966. 54. Liu, P.F.; Zheng, J.Y. (2010): Recent developments on damage modeling and finite element analysis for composite laminates: A review. In Materials & Design 31 (8), pp. 38253834. 55. Kachanov L. M. (1958): Time of the rupture process under creep conditions. Izv. Akad. Nauk. SSR, Otd Tekh. Nauk No. 8, pp. 26-31. 56. Hult, J. (1979): Continuum damage mechanics  Capabilities limitations and promises. Mechanisms of Deformation and Fracture, Pergamon, Oxford, pp. 233-347. 57. Chaboche, J. L. (1981): Continuum damage mechanics: A tool to describe phenomena before crack initiation. Nuclear Engineering and Design, Vol. 64, pp. 233-247. 58. Krajcinovic, D. (1984): Continuum damage mechanics. Applied Mechanics Review, Vol. 37, No.1. pp. 1-6. 59. Zhang, Wohua; Cai, Yuanqiang (2010): Continuum damage mechanics and numerical applications. Berlin, New York, Hangzhou: Springer; Zhejiang University Press (Advanced topics in science and technology in China). 60. Charboche, J. L. (1988): Continuum damage mechanics: Part I - General concepts. In Journal of Applied Mechanics 55 (1), pp. 5964. 61. Chaboche, J. L. (1988): Continuum damage mechanics: Part II - Damage growth, crack initiation, and crack growth. In Journal of Applied Mechanics 55 (1), pp. 65  174 72. 62. Krajcinovic D., Lemaitre J., Hult J., Murakami S., Sumarac D., Leckie FA, Najar J. (eds.) (1986): Continuum Damage Mechanics--Theory and Applications. Springer Verlag, Berlin. 63. Krajcinovic, Dusan (2000): Damage mechanics: accomplishments, trends and needs. In International Journal of Solids and Structures 37 (1-2), pp. 267277. 64. Park, Sun Woo; Richard Kim, Y.; Schapery, Richard A. (1996): A viscoelastic continuum damage model and its application to uniaxial behavior of asphalt concrete. In Mechanics of Materials 24 (4), pp. 241255. 65. Kaddour, A.; Hinton, M.; Smith, P.; Li, S. (2013): A comparison between the predictive capability of matrix cracking, damage and failure criteria for fibre reinforced composite laminates: Part A of the third world-wide failure exercise. In Journal of Composite Materials 47 (20-21), pp. 27492779. 66. Murakami, S.; Kamiya, K. (1997): Constitutive and damage evolution equations of elastic-brittle materials based on irreversible thermodynamics. In International Journal of Mechanical Sciences 39 (4), pp. 473486.  67. Schapery, R. A. (1990): A theory of mechanical behavior of elastic media with growing damage and other changes in structure. In J. Mech. Phys. Solids 38 (2), pp. 215253.  68. Pineda, Evan J.; Waas, Anthony M.; Bednarcyk, Brett A.; Collier, Craig S.; Yarrington, Phillip W. (2009): Progressive damage and failure modeling in notched laminated fiber reinforced composites. In Int J Fract 158 (2), pp. 125143. 69. Pineda, Evan J.; Waas, Anthony M. (2013): Numerical predictions of damage and failure in carbon fiber reinforced laminates using a thermodynamically-based work potential theory. In Proceeding of ASC 28th technical conference, Sep. 9-11, 2013, University Park, Pennsylvania, USA. 70. Talreja, R. (1985): A Continuum Mechanics Characterization of Damage in Composite Materials. In Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 399 (1817), pp. 195216. 71. Talreja, R. (1994): Damage characterization by internal variables. In: Talreja R (ed) Damage mechanics of composite materials. Composite materials series, vol 9. Elsevier Science, Oxford, Chap. 2. 72. Puck, A.; Schürmann, H. (2002): Failure analysis of FRP laminates by means of physically based phenomenological models. In Composites Science and Technology 62 (12-13), pp. 16331662.   175 73. Flatscher, Th.; Schuecker, C.; Pettermann, H. E. (2009): Prediction of plastic strain accumulation in continuous fiber reinforced laminates by a constitutive ply model. In Int J Fract 158 (2), pp. 145156. 74. Pinho, S.T.; Iannucci, L.; Robinson, P. (2006): Physically-based failure models and criteria for laminated fibre-reinforced composites with emphasis on fibre kinking: Part I: Development. In Composites Part A: Applied Science and Manufacturing 37 (1), pp. 6373. 75. Davila, C. G. (2005): Failure Criteria for FRP Laminates. In Journal of Composite Materials 39 (4), pp. 323345. 76. Pinho, S.T.; Iannucci, L.; Robinson, P. (2006): Physically based failure models and criteria for laminated fibre-reinforced composites with emphasis on fibre kinking. Part II: FE implementation. In Composites Part A: Applied Science and Manufacturing 37 (5), pp. 766777. 77. Pinho, S.; Darvizeh, R.; Robinson, P.; Schuecker, C.; Camanho, P. (2012): Material and structural response of polymer-matrix fibre-reinforced composites. In Journal of Composite Materials 46 (19-20), pp. 23132341. 78. Pinho, S.; Vyas, G.; Robinson, P. (2013): Response and damage propagation of polymer-matrix fibre-reinforced composites: Predictions for WWFE-III Part A. In Journal of Composite Materials 47 (20-21), pp. 25952612. 79. Williams, Kevin V.; Vaziri, Reza; Poursartip, Anoush (2003): A physically based continuum damage mechanics model for thin laminated composite structures. In International Journal of Solids and Structures 40 (9), pp. 22672300. 80. Mcgregor, C.; Zobeiry, N.; Vaziri, R.; Poursartip, A. (2008): A Constitutive Model for Progressive Compressive Failure of Composites. In Journal of Composite Materials 42 (25), pp. 26872716. 81. Maimí, P.; Camanho, P.P.; Mayugo, J.A.; Dávila, C.G. (2007): A continuum damage model for composite laminates: Part II  Computational implementation and validation. In Mechanics of Materials 39 (10), pp. 909919. 82. Ladeveze P and Le Dantec E, Damage modeling of the elementary ply for laminated composites, Composites Science and Technology 43 (1992) 257-267. 83. Allix O, Gilletta D,  Ladeveze P, Non linear mechanical behavior of laminates, Proceedings of the EUROMECH Colloquium on Structure and Crack Propagation in Brittle Matrix Composite Materials, 1986, P227-240 84. Ladeveze P, Allix O, Douchin B and Lévêque D, A computational method for damage intensity prediction in a laminated composite structure, Computation Mechanics, New Trends and Applications, Barcelona, Spain 1998   176 85. Allix O, A composite damage meso-model for impact problems, Composite Science and Technology 61 (2001) 2193-2205. 86. Ladeveze O and Lubineau G, On a damage mesomodel for laminates: micro-meso relationships, possibilities and limits, Composites Science and Technology 61(2001) 2149-2158. 87. Ladeveze O and Lubineau G, On a damage mesomodel for laminates: micromechanics basis and improvement, Mechanics of Materials 35 (2003) 763-775. 88. Ladeveze P, Multiscale modeling and computational strategies for composites, International Journal for Numerical Methods in Engineering, 2004, 60:233-253. 89. Ladeveze P, Lubineau G and Violeau D, A computational damage micromodel of laminated composites, International Journal of Fracture (2006) 137:139-50 90. Lubineau G and Ladeveze P, Construction of a micromechanics-based intralaminar mesomodel, and illustrations in ABAQUS/Standard, Computational Materials Science 43 (2008) 137-145. 91. Allix O, Feissel P and Thevenet P, A delay damage mesomodel of laminates under 81 (2003) 11771191. 92. G. Lubineau and P. Ladeveze, Construction of a micromechanics-based intralaminar mesomodel, and illustrations in ABAQUS/Standard, Computational Materials Science 43 (2008) 137-145. 93. Chen JF, Morozov EV and Shankar K, A combined elastoplastic damage model for progressive failure analysis of composite materials and structures, Composite Structures 94 (2012) 34783489. 94. Damage and Plasticity Parameters for Continuum Damage Mechanics Modelling of Carbon and Glass Fibre-Reinforced Composite Materials. In Strain 47 (1), pp. 105115. 95. Allix, O.; Ladevèze, P.; Vittecoq, E. (1994): Modelling and identification of the mechanical behaviour of composite laminates in compression. In Composites Science and Technology 51 (1), pp. 3542. 96. Allix, O.; Bahlouli, N.; Cluzel, C.; Perret, L. (1996): Modelling and identification of temperature-dependent mechanical behaviour of the elementary ply in carbon/epoxy laminates. In Composites Science and Technology 56 (7), pp. 883888. 97. Phillips EA, Herakovich CT and Graham LL, Damage development in composites with large stress gradients, Composites Science and Technology 61 (2001) 2169  177 2182. 98. Johnson, A.F. (2001): Modelling fabric reinforced composites under impact loads. In Composites Part A: Applied Science and Manufacturing 32 (9), pp. 11971206. 99. Johnson, A.F; Pickett, A.K; Rozycki, P. (2001): Computational methods for predicting impact damage in composite structures. In Composites Science and Technology 61 (15), pp. 21832192. 100. Johnson, Alastair F.; David, Matthew (2010): Failure mechanisms in energy-absorbing composite structures. In Philosophical Magazine 90 (31-32), pp. 42454261. 101. Pickett, A.K.; Fouinneteau, M.R.C. (2006): Material characterization and calibration of a meso-mechanical damage model for braid reinforced composites. In Composites Part A: Applied Science and Manufacturing 37 (2), pp. 368377. 102. Fouinneteau, M.R.C.; Pickett, A.K. (2007): Shear mechanism modelling of heavy tow braided composites using a meso-mechanical damage model. In Composites Part A: Applied Science and Manufacturing 38 (11), pp. 22942306. 103. Justin Littell (2008): The experimental and analytical characterization of the macromechanical response for triaxial braided composite materials. Ph.D. Dissertation. University of Akron, Akron, Ohio. 2008. 104. Yerramalli C. S., Waas, A. M. (2003): A failure criterion for fiber reinforced polymer composites under combined compression-torsion loading. In International Journal of Solids and Structures 40, pp. 1139-1164. 105. Song, Shunjun; Waas, Anthony M.; Shahwan, Khaled W.; Xiao, Xinran; Faruque, Omar (2007): Braided textile composites under compressive loads: Modeling the response, strength and degradation. In Composites Science and Technology 67 (15-16), pp. 30593070. 106. Nahas, M. N. (1986): Survey of Failure and Post-Failure Theories of Laminated Fiber-Renforced Composites. In J. Compos. Technol. Res. 8 (4), pp. 138153. 107. Prabhakar, Pavana; Waas, Anthony M. (2013): Upscaling from a micro-mechanics model to capture laminate compressive strength due to kink banding instability. In Computational Materials Science 67, pp. 4047. 108. Kongshavn, Ingrid; Poursartip, Anoush (1999): Experimental investigation of a strain-softening approach to predicting failure in notched fibre-reinforced composite laminates. In Composites Science and Technology 59 (1), pp. 2940. 109. Zobeiry, N.; Vaziri, R.; Poursartip, A. (Eds.) (2013): Strain-Softening Response of Laminated Composites under Compression. The 19th International Conference on   178 Composite Materials. Montreal, Canada. 110. Basu, S.; Am Waas; Ambur (2006): Compressive failure of fiber composites under multi-axial loading. In JOURNAL OF THE MECHANICS AND PHYSICS OF SOLIDS 54 (3), pp. 611634. 111. Feld, N.; Allix, O.; Baranger, E.; Guimard, J.-M. (2012): A micromechanics-based mesomodel for unidirectional laminates in compression up to failure. In JOURNAL OF COMPOSITE MATERIALS 46 (23), pp. 28932909. 112. Hahn, H. T. and Tsai, S. W., On the Behaviour of Composite Laminates After Initial Failure, Journal of Composite Materials, Vol. 8, July 1974, pp. 280305. 113. Chiu, K. D., Ultimate Strengths of Laminated Composites, Journal of Composite Materials, Vol. 3, July 1969, pp. 578582. 114. Petit, P. H. and Waddoups, M. E., A Method of Predicting the Non-Linear Behaviour of Laminated Composites, Journal of Composite Materials, Vol. 3, Jan. 1969, pp. 219. 115. Nahas, M. N. (1986): Yield and ultimate strengths of fibre composite laminates. In Composite Structures 6 (4), pp. 283294.  116. Flatscher, Th.; Pettermann, H.E. (2011): A constitutive model for fiber-reinforced polymer plies accounting for plasticity and brittle damage including softening  Implementation for implicit FEM. In Composite Structures 93 (9), pp. 22412249. 117. Mazzeranghi, A.; Vangi, D. (1994): A post-failure model for composite laminates based on phenomenologic aspects of damage. In Composite Structures 28 (3), pp. 323332. 118. Williams, Kevin V.; Vaziri, Reza (2001): Application of a damage mechanics model for predicting the impact response of composite materials. In Computers & Structures 79 (10), pp. 9971011. 119. Randles, P. W.; Nemes, J. A. (1992): A continuum damage model for thick composite materials subjected to high-rate dynamic loading. In Mechanics of Materials 13 (1), pp. 113. 120. Nemes, J.A.; Spéciel, E. (1996): Use of a rate-dependent continuum damage model to describe strain-softening in laminated composites. In Computers & Structures 58 (6), pp. 10831092.  121. Matzenmiller A, Lubliner J, Taylor R L (1995), A constitutive model for anisotropic damage in fiber composites. Mechanics of Materials 20, pp. 125152. 122. Williams K, Vaziri R (1995). Finite element analysis of the impact response of   179 CFRP composite plates. Proceedings of the 10th International Conference on Composite Materials, Whistler, Canada, pp. 647654. 123. Yen, C.F. (Ed.) (2002): Ballistic Impact Modeling of Composite Materials. Proceedings of 7th International LS-DYNA Users Conference. Dearborn, Michigan. pp. 6. 15-6.26. 124. -DYNA MAT162 Unidirectional and Plain Weave Composite Progressive Failure Models, December 2011 125. Deslauriers, Paul F.; Cronin, Duane S.; Duquette, Alex (Eds.) (2004): Numerical Modeling of Woven Carbon Composite Failure. Proceedings of 8th international LS-DYNA Users Conference. Detroit, Michigan. 126. Brown, Kevin; Brooks, Richard; Warrior, Nicholas (Eds.) (2005): Numerical Simulation of Damage in Thermoplastic Composite Materials. Proceedings of 5th European LS-DYNA Users Conference. Birmingham, England. 127. Xiao, J. R.; Gama, B. A.; Gillespie, J. W. (2005): Progressive Damage and Delamination in Plain Weave S-2 Glass/SC-15 Composites Under Quasi-Static Punch Shear Loading. In Composite Structures, pp. 163171. 128. Hoof, Jack van; Worswick, Michael J.; Straznicky, Paul V.; Bolduc, Manon (Eds.) (1999): Effects of post failure modelling on the response of ballistically impacted composites. Proceedings of the 12th international conference on composite materials. Paris, France. 129. Batra, R.C.; Gopinath, G.; Zheng, J.Q. (2012): Damage and failure in low energy impact of fiber-reinforced polymeric composite laminates. In Composite Structures 94 (2), pp. 540547. 130. Tabiei, Ala; Aminjikarai, S. Babu (2009): A strain-rate dependent micro-mechanical model with progressive post-failure behavior for predicting impact response of unidirectional composite laminates. In Composite Structures 88 (1), pp. 6582. 131. Peng (1984): Continuum Theory for StrainSoftening. In J. Eng. Mech. 110 (12), pp. 16661692. 132. Read, H.E.; Hegemier, G.A. (1984): Strain softening of rock, soil and concrete  a review article. In Mechanics of Materials 3 (4), pp. 271294. 133. Costin, L.S. (1985): Damage mechanics in the post-failure regime. In Mechanics of Materials 4 (2), pp. 149160. 134. Yalla Abushawashi, Xinran Xiao, Viktor Astakhov, A Novel Approach For   180 Determining Material Constitutive Parameters For A Wide Range Of Triaxiality Under Plane Strain Loading Conditions, International Journal of Mechanical Sciences, 74(2013)133142. 135. Alejano, L.R.; Rodriguez-Dono, A.; Alonso, E.; Fdez.-Manín, G. (2009): Ground reaction curves for tunnels excavated in different quality rock masses showing several types of post-failure behaviour. In Tunnelling and Underground Space Technology 24 (6), pp. 689705. 136. thesis, Stanford University. 137. hybrid square sandwich composite vehicle hollow body shells with reinforced core 86. 138. -Japan conference on composite materials, Dearborn, Michigan. 139. A. Haufe, K. Schweizerhof, P. Dubois, Properties & Limites: Review of shell element formulations. Sep. 2013 140. LS- 141.  142. Hautefeuille, M.; Colliat, J.-B.; Ibrahimbegovic, A.; Matthies, H.G.; Villon, P. (2012): A multi-scale approach to model localized failure with softening. In Computers & Structures 94-95, pp. 8395. 143. Bowman, C. L., Roberts, G. D., Braley, M. S., Xie, M., and Booker, M. J. (2002): Mechanical properties of triaxial braided carbon/epoxy composites. NASA TM2002-211702. 144. Bazant, Z. P.; Planas, Jaime. (1998): Fracture and size effect in concrete and other quasibrittle materials. Boca Raton: CRC Press.