MODE I AND MODE II INTERLAMINAR FRACTURE TOUGHNESS SIMULATION OF U NIDIRECTIONAL AND QUASI - THREE - DIMENSIONAL COMPOSITE S By Xinyu Mao A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Mechanical Engineering - M aster of S cience 2019 ABSTRACT MODE I AND MODE II INTERLAMINAR FRACTURE TOUGHNESS SIMULATION OF UNIDIRECTIONAL AND QUASI - THREE - DIMENSIONAL COMPOSITE S By Xinyu Mao Fiber reinforced polymer (FRP) composite materials are widely used in automotive and aerospace industries. The traditional FRP composites have a laminated construction in which layers are reinforced with fibers in unidirectional (UD) or fabric forms. Such materials possess good mechanical properties in - plane but relatively weak in the through - thickness direction. Due to their low interlaminar strength, composites are prone to delamination. The quasi - three - dimensional (Q3D) composite are designed to improve the interlaminar strength. In this work, the Q3D composite is made by a special braiding process in which the bias tows are braided into adjacent layers. Contrary to the conventional composites in which plies are bonded by the polymer matrix only, the pli es of Q3D composites are bridged by fiber tows. As a higher load is required to break the bridging tows for crack growth, the delamination resistance will increase. This has been proven by interlaminar fracture experiments under Mode I and Mode II conditio ns. This thesis is focused on numerical simulations of UD and Q3D composites under Mode I and Mode II loadings. This investigation is needed for the development of simulation methods for crash safety simulations of Q3D composite structures in automotive applications. In this work, the interlaminar delamination was modeled with cohesive elements. Both the Bilinear and Trilinear Cohesive Zone Models (CZM) were investigated. The CZM parameters were determined from the Mode I and Mode II fracture toughness va lues measured in experiments. The simulations were performed using explicit finite element (FE) code LS - DYNA. It was observed that in Mode I simulations, the Bilinear CZM predicted a stable crack growth, which agreed with the experimental observation. On the other hand, the Trilinear CZM predicted a relatively unstable crack growth. In Mode II experiment with an end - notched flexural (ENF) configuration, the delamination tends to be unstable . This behavior was better predicted with the Trilinear CZM. In FE models for Q3D composite, the cohesive elements were assigned with two sets of CZM parameters. The first set of parameters was the same as that for the UD model. The second set with higher facture toughness values was assigned to the cohesive elements cor responding to the bridging tows. This method captured the stick - slip behavior of the Q3D composites observed in Mode I experiments. In general, the prediction was better for the UD composite than for the Q3D composite. The predicted load at the delaminatio n initiation was within ±10% of the experimental values for the UD and within ±15% for the Q3D. iv ACK NOWLEDGE MENTS I would like to express my sincere gratitude to m y advisor Dr. Xinran Xiao for her help , advice and guidance throughout my graduate research work. She not only taught me how to do research a nd review literature papers, but also shared her life experience with me. I learned how to become a better employee, a better master student and a better person from her. Without her guidance and help , I would be impossible for me to complete this thesis . I also wat to thank Dr. and Dr. fo r participating in my thesis committee. I am thankful to people who helped me throughout my research. I want to give special thanks to Shutian Yan who shared her knowledge and experience on numerical simulation with me . When I have problems on simulation, s he is always willing to answer my questions and provide me a lot of valuable resources to read . I also want to thank Zhou Wu and Tony Wente who taught me a lot of experiment methods . Tony provide d me with experimental data and specimen dimensions which helped me tremendously throughout my research . Lastly, I wan t to express my appreciation to Dr. Danielle Zeng of Ford Motor Company for her support, feedback s and suggestions during my graduate research. I am so grateful to my parents for providing me with the opportunity to receive great ed ucation in United States. They sacrificed a lot for me to have this opportunity studying abroad. I appreciate my friends: Chenyu Mao, Ran Bi, Haozhen Lyu for appearing in my life. They always stand by my sid e when I am going through tough time . v TABLE OF CONTENT S LIST OF TABLES ................................ ................................ ................................ ........................ vii LIST OF FIGURES ................................ ................................ ................................ ..................... viii Chapter 1. Introduction ................................ ................................ ................................ .................... 1 1.1 Fiber - Reinforced Composite Materials ................................ ................................ ............... 2 1.1.1 Characteristics of Composite ................................ ................................ ...................... 2 1.1.2 Failure of Composite ................................ ................................ ................................ ... 3 1.2 Quasi - 3D Composite Material ................................ ................................ ............................ 4 1.2.1 Characteristics of Q3D Composite ................................ ................................ ............. 4 1.2.2 Manufacturing of Q3D Composite ................................ ................................ ............. 5 1.3 Scopes of Work ................................ ................................ ................................ ................... 7 1.4 Scopes of Thesis ................................ ................................ ................................ ................. 7 Chapter 2. Literature Review ................................ ................................ ................................ ........... 9 2.1 Interlaminar Fracture Toughness Simulation Methods ................................ ....................... 9 2.2 Cohesive Damage Models ................................ ................................ ................................ .. 9 2.3 Simulation of Z - pinned and Stitched Composites ................................ ............................ 12 Chapter 3. LS - DYNA Simulation ................................ ................................ ................................ .. 13 3.1 Numerical Model Setup ................................ ................................ ................................ .... 14 3.1.1 Model Dimensions and Structures for Mode I ................................ .......................... 14 3.1.2 Model Dimensions and Structures for Mode II ................................ ........................ 17 3.2 Cohesive Zone Law ................................ ................................ ................................ .......... 19 3.2.1 MAT_138 Bilinear Cohesive Zone Law ................................ ................................ ... 19 3.2.2 MAT_185 Trilinear Cohesive Zone Law ................................ ................................ . 25 3.3 Material Properties ................................ ................................ ................................ ............ 28 Chapter 4. Simulation Results and Discussion ................................ ................................ .............. 31 4.1 Mode I Results ................................ ................................ ................................ .................. 31 4.1.1 UD Composite with Bilinear CZM ................................ ................................ ........... 31 4.1.2 UD Composite with Trilinear CZM ................................ ................................ .......... 37 4.1.3 Q3D Composite ................................ ................................ ................................ ........ 41 4.2 Mode II Results ................................ ................................ ................................ ................. 43 4.2.1 UD Composite with Bilinear CZM ................................ ................................ ........... 43 4.2.2 UD Composite with Trilinear CZM ................................ ................................ .......... 45 4.2.3 Q3D Composite ................................ ................................ ................................ ........ 48 vi Chapter 5. CONCLUSION ................................ ................................ ................................ ............ 50 5.1 Summary and Conclusion ................................ ................................ ................................ . 50 5.2 Future Work ................................ ................................ ................................ ...................... 54 APPENDIX ................................ ................................ ................................ ................................ .... 55 BIBLIOGRAPHY ................................ ................................ ................................ .......................... 57 vii L IST OF TABLE S Table 3.1 Dimensions of Mode I DCB model .. 16 Table 3.2 Dimensions of Mode II ENF model .. 18 Table 3.3 Mode I fracture toughness values measured in DCB experiments [ 6 .. . . 22 Table 3.4 Mode II fracture toughness values measured in ENF experiments [ 6 ] 22 Table 3.5 Fracture toughness values of UD composite in Mode I and II .. 23 Table 3.6 Fracture toughness values of Q3D composite for Mode I .. 23 Table . 24 Table 3.8 Scaled Length parameters for UD composites in DCB and ENF simulations 26 Table 3.9 Scaled Length parameters for Q3D composites in ENF simulation 26 Table 3.10 Carbon fiber properties 29 Table 3.11 SC - 15 Epoxy properties 29 Table 3.12 Composite mechanical properties 29 Table 5.1 Comparison of the predicted and measured peak lo ad with Bilinear CZM 5 2 Table 5.2 Comparison of the predicted and measured peak load with Trilinear CZM 5 2 Table 5.3 Comparison of the predicted and measured displacement at the peak load with Bilinear CZM 5 3 Table 5.4 Comparison of the predicted and measured displacement at the peak load with Trilinear CZM 5 3 Table A.1 Different Cohesive Zone l aws 5 6 viii LIST OF FIGURE S Figure 1 .1 Fiber - Reinforced Composite Materials [ 8 ] . 2 Figure 1. 2 (a) Front - top view of Q3D composite (b) Cross - section view of Q3D composite . . ... 5 Figure 1.3 (a) Completed Q3D carbon fiber preform (b) Schematic of VARTM process [ 13 ]. 6 Figure 2.1 Schematic of DCZM [ 29 ] 11 Figure 2.2 load - displacement curve of stitch elements [ 33 ] 1 Figure 3.1 Numerical Model setup of UD Composite .. 14 Figure 3.2 (a) Interface of t ested Q3D c omposite l aminate (b) Cohesive l ayer interface of Q3D c omposite l aminate .. ... . ... 15 Figure 3.3 Front view of Mode II ENF model setup .. 17 Figure 3.4 Cohesive layer setup of Q3D composite laminate for Mode II 18 Figure 3.5 (a) MAT_138 Bilinear Traction - separation curve (b) MAT_185 Trilinear Traction - separation curve . 19 Figure 3.6 Displacement in cohesive element . 20 Figure 3.7 Integration points on a cohesive element 0 Figure 3.8 Architecture of the Q3D composite specimen for DCB and ENF tests [ 6 ] .. ... 28 Figure 4.1 Comparison of the experimental load - extension responses of UD DCB specimen under Mode I with the prediction with thick cohesive element with a Bilinear CZM law . 3 1 Figure 4.2 Comparison of the experimental load - extension responses of UD DCB specimen under Mode I with the prediction with thin cohesive element with a Bilinear CZM law . 32 Figure 4.3 Comparison of Bilinear CZM law with the same but different traction values . 33 Figure 4.4 Comparison of the load - extension curves by simulations with the two CZMs . .. 33 Figure 4.5 The experimental R curves for UD composites 35 Figure 4.6 Mode I Experiment picture of UD composite laminate 35 Figure 4.7 Comparison of Bilinear CZM with initial and Bilinear CZM with average ... . . 36 ix Figure 4.8 Load - extension curve of Trilinear CZM with thick cohesive layer 37 Figure 4.9 Load - extension curve of Trilinear CZM with thin cohesive layer . 38 Figure 4.10 Comparison of traction - separation curves for Bilinear CZM and Trilinear CZM .. .. . 39 Figure 4.11 Load - extension curve comparison of Bilinear CZM and Trilinear CZM . 39 Figure 4.12 Cohesive element deletion of Bilinear CZM . 40 Figure 4.13 Cohesive element deletion of Trilinear CZM 0 Figure 4.14 Load - Extension curve of Q3D composite laminate with thick cohesive layer . 4 1 Figure 4.15 Load - Extension curve of Q3D composite laminate with thin cohesive layer 42 Figure 4.16 Comparison of the ENF experimental Load - deflection curve of UD com posite with simulation using the Bilinear CZM .. 43 Figure 4.17 CZM comparison with different stiffness . 44 Figure 4.18 Load - deflection curve comparison with two CZMs 44 Figure 4.19 Comparison of ENF experimental Load - deflection curve with simulation with a Trilinear CZM 45 Figure 4.20 Traction - separation curves compar ison for Bilinear and Trilinear CZMs . 46 Figure 4.21 Load - extension curves comparison for Bilinear and Trilinear CZMs . .. 46 Figure 4.22 Load - extension curve of Bilinear CZM for Mode II 47 Figure 4.23 Load - extension curve of Trilinear CZM for Mode II 48 1 Chapter 1. Introduction Fiber reinforced polymer (FRP) c omposite materials have been widely used in automotive and aerospace industr ies since 19 6 0s. C ompared with conventional metallic materials, c omposite s have a relatively high specific strength, specific stiffness and fatigue resistance [ 1 ]. Implement ing composite materials in vehicle structures can lead to significant improvement in fuel efficiency and better acceleration. F RP composites have been used widely in the chassis and body of supercars and racing cars. FRP s have relatively strong in - plane strength . However, the through - thickness direction strength of these materials is relatively low . The most common failure mode in FRPs is delamination. Out - of - plane impact could cause invisible damage to com posit e s which result s in delamination failure in the through - thickness direction. To improve the interlaminar fracture toughness of a composite, methods of z - pinning and stitching [ 2 - 5 ] in the through - thickness direction are developed . This improvement is usually measured by interlaminar fracture toughness tests under the M ode I and M ode II conditions . In Mode I and II fracture toughness tests, z - pin n ing and stitching methods show better interlaminar fracture toughness strength . However, the degradation of composite in - plane strength is also observed in experiment s . To improve the interlaminar delamination resistance and remedy the in - plane property degradation , Quasi three - dimensional (Q3D) composite s have been developed which us es a special triaxial braiding technique to insert bias tows into adjacent laye r s . Each ply is connected by bias tows 2 throughout the thickness of composite . Wente et al. [ 6 ] and Zhou [ 7 ] found that Q3D composites significantly improve interlaminar fracture toughness in the through - thickness direction of composite s with small in - plane property degradation . The objective of this work is to develop a n umerical simulation method to model the delamination behavior of Q 3 D composite s . The use of c ohesive element to model the delamination behavior is investigated. Two cohesive laws, i.e. the Bilinear Cohesive Zone Law and Trilinear Cohesive Zone Law , are compared . T he results of this work will be useful in the development of simulation method for crashworthiness simulations of Q3D automotive structures. 1.1 Fiber - Reinforced Composite Materials 1.1.1 Characteristics of Composite Composite materials consist of fiber and matrix . Figure 1 .1 demonstrates the schematic s of fiber - reinforced composite material. Different types of fiber s , such as glass fiber and carbon fiber, ha ve different material properties . Matrix can be pick ed from various types of materials such as polymer, metal and ceramics [ 8 ] . Diverse combination s of fiber and matrix result in unique composite material properties . L amina that is made of fibers and matrix can become composite laminate by Figure 1 .1 Fiber - Reinforced Composite Materials [8] 3 overlapping each other. The fiber preforms may be manufactured by various ways such as u nidirectional, stitching, w oven , braid ing , etc. Each pattern gives speci al composite characteristics. Various orientation s of laminas when stacking to compose laminate, also develop unique material properties after the composite is fabricated. 1.1.2 Failure of Composite Composite components of airplane s and automobile s frequently undergo cyclical load that can cause degradation of the composite structure. T o prevent catastrophic failure , it is critical to study the failure mode of composite s [ 9 ] . Failure modes of composite laminates can be categorized by fiber failure, matrix failure , delamination failure and other types of failure. Composite f iber failure occur s when the composite is subjected to tension and compression in the direction of the fiber . Tensile load in the fiber direction could cause fiber breakage inside the composite . A c rack can propagate in the matrix which can lead to matrix failure. Under compressive load, the common failure modes are fiber kinking and buckling. Another common matrix failure is due to shear . Composite s ha ve relatively weak out of plane strength compar ed with conventional isotropic materials. C omposite material has a low ability to absor b kinetic energy a pplied in an ou t - of - plane direction [ 10 ] . Impact on a composite could create matrix damage inside the composite laminate . The crack occurs and propagates along the interface of the composite laminate eventually result ing in delamination of the laminate. 4 1.2 Quasi - 3D Composite Material Conventional composite lamina is made by embedding either unidirectional , 2D bi - axial , or even 2D triaxial fiber preform in epoxy . Although these fiber patterns provide relatively strong in - plane strength in the fiber direction, the delamination resistance in its out - of - plane direction become the biggest weakness. T he objective of developing Q3D composite is to improve the interlamina r strength by insert ing bias tow s into adjacent layer to create fiber bridging between adjacent layers of laminate s . Instead of having multi ple layers of unidirectional or 2D woven lamina stacking on each other , Q3D laminate integrate s each lamina to a uniform three - dimensional woven preform. Unlike conventional weaving technique, Q3D preform s has fiber laying down in three axes directions [ 11 ] . 1.2.1 Characteristics of Q3D Composite The failure mode s of composite material s are highly depend ing on interlaminar integr ation of laminate s . Q3D composite s improve the interlaminar integration by attaching adjacent plies with fiber tows. Although the in - plane properties of Q3D are slightly lower than unidirectional (UD) and 2D woven composite s , the delamination resistance improves significantly compar ed with UD and 2D woven configurations. The reinforcement in the through - thickness direction help s composite laminate absorb ing more specific energy due to impact [ 12 ] . The residual strength is high er than UD and 2D woven composite s after impact - induced damage . When Q3D composite s are unde r compressive load, buckling failure appears rather than shear failure [ 1 ]. Wente et al [ 6 ] indicates that Q3D composite s ha ve higher fracture toughness than UD and 2D woven c omposite laminate s after initiation of the crack propagation . According to Mode I DCB test, b ecause of fiber 5 bridging tows in Q3D composite laminate s , it requires extra work and energy to break the fiber bridging tows in composite laminates . Mode II ENF test that done by Wente et al. [ 6 ] also shows that Q3D composites has high fracture toughness . 1.2.2 Manufacturing of Q3D Composite The Q3D composite lamina i n this thesis consists of fiber tow s in three direction s . These are 0°, 60° and - 60°. The a xial carbon fiber tow lays down in 0° direction. Based on 2D woven preform, the 60° bias tow of 2D woven preform in current layer is braided into the next adjacent layer for the Q3D composite preform . In this case, two layers of 2D woven fabric are connect ed by 60° bias tow. Figure 1.2 shows microstructure of the Q3D composite that is created using TEXGEN. The green bias tow demonstrates inserted fiber tow which makes fiber bridging occur between adjacent lamina. The Q3D carbon fiber preform is shown in Figure 1.3 (a). Strips highlighted in red represent 0° axial tow. Yellow and green strips are 60° bridging bias tow and - 60° bias tow. The Q3D composite is manufactur ed by vacuum assisted resin transfer molding (VARTM) method . VARTM is commonly used by in dustry to manufacture large size and low volume part s . VARTM process utilizes vacuum to guide resin flow through the fiber preform. To ensure the good quality Figure 1.2 (a) Front - top view of Q3D composite (b) Cross - section view of Q3D composite 6 of composite plate s and avoid occurrence of voids , perfect vacuum condition inside the vacuum bag is crucial [ 13 ]. Figure 1.3 (b) demonstrates the schematic of VARTM process. Q3D preform that is made up of 12K , A - 42 carbon fiber from DowAksa is placed on an a luminum base plate. The matrix , that is used to manufacture Q3D composite , is composed of SC - 15A resin and SC - 15B hardener with mix ratio of 10:3. Transfer media (distribution media), which helps resin fully impregnate with carbon fiber preform, is placed on top of the peel ply. The p eel ply , which is inserted between the transfer media and composite preform , provide easy demolding after curing. Vacuum bag seals on top of the distribution media. Pressure bucket with vacuum assist pulls the resin from left to right of the Omega tube . After curing for 8 hours in auto - clave, the Q3D composite plate is manufactured. Figure 1.3 (a) Completed Q3D carbon fiber preform (b) Schematic of VARTM process [13] 7 1.3 Scopes of Work The importan ce of studying the interlaminar strength of composite material s is well recognized. To prevent catastrophic structural failure of composites in the through - thickness direction , study ing interlaminar strength , especially delamination behavior , of composite s is significant . The improvement of the delamination resistance in the out - of - plane direction of Q3D composite s is proved by experiment . T e st ing the performance of composite material s is relatively challenging . To manufacture perfect composite specimen s requires a lot of practice . For the sake of saving time and cost on conducting experiments , it is more efficient to numerically predict system of failure beyond elastic region of laminated composite material s . Commercial finite element software s such as LS - DYNA, ABAQUS, ANASYS are develop ed to predict failure mode, post - failure behavior of composite materials . The numerical simulation inc orporate s composite material properties , which can be obtained experimentally , into composite models . Th e objective of this research is using finite element method to simulat e interlaminar fracture toughness tests and investigate the delamination resistance of Q3D composite s . LS - DYNA software is used in this research to build a FE model to simulate M ode I DCB and M ode II ENF tests of Q 3D composite material s using bilinear and trilinear cohesive z one model s . The load - extension curve is compared with experimental results to check prediction accuracy. 1.4 Scopes of Thesis In this thesis, numerical models are developed using bilinear and trilinear cohesive zone model s to represent delamination resistance o f Q3D composite laminates under Mode I and Mode II 8 condition s . The input parameters of numerical models are obtained from experiments [15] and literature papers [ 1 4 ] . Th is thesis is organized as follow s: Chapter 1 introduces the general characteristics and failure modes of fiber - reinforced composite material s. The Quasi - 3D composite is introduced in this chapter. Chapter 2 provides literature review on numerical method of simulating interlaminar fracture toughness tests . The cohesive zone models in literature paper are studied and summarized in this chapter. Chapter 3 describes LS - DYNA simulation work . The numerical model s etup s and dimen sions are provided for both Mode I and Mode II simulations . The governing equations and cohesive laws for cohesive zone model s are presented . The material and mechanical properties of the composite that is used in the simulation are shown in t his chapter. Chapter 4 The results of Mode I and Mode II simulations using bilinear and trilinear cohesive zone model are presented . The load - extension curves of simulation are compared and discussed with experiment al results . The initial and average interlaminar fracture toughness values are used to investigate their influences on the prediction results. The effects on prediction results of using two different definitions of shell thickness reference plane are also investigate d. Chapter 5 concludes and summarizes the findings in this research . The future work and improvements are suggested. 9 Chapter 2. Literature Review 2.1 Interlaminar Fracture Toughness Simulation Method s The fiber - reinforced composites have low interlaminar strength in the through - thickness direction. One of the common failure modes of fiber - reinforced composite laminates is delamination. There are two common approaches to numerically predict the delamination in composite laminates. Linear Elastic Fracture Mechanics (LEFM) is one of the methods to simulat e the delamination behavior . Liu and Zou et al. [ 15 - 16 ] demonstrated that LEFM , particularly the Virtual Crack Closure Technique (VCCT) [ 16 - 19 ], can be used to model the delamination crack growth with a n initial crack . The downside of this approach is that it is time - consuming and difficult to complete the calculation of fracture parameters such as energy release rate for progressive crack growth . Another approach is to use damage mechanic s . C ohesive Z one M ethod (CZM) belongs to this approach. It is simpler to implement CZM in finite element models than VCCT. 2.2 Cohesive Damage Models In C ZM method, different traction - separation laws have been developed to study the fracture process of various materials [ 20 - 22 ]. These CZM laws have been employed in delamination modeling of composite materials. Jung et al. [ 2 3 ] buil t an intralaminar damage model using continuum damage me chanics [ 2 4 ] and interlaminar damage model using cohesive zone method for glass fiber - reinforced polypropylene (GFPP) composite s. The mix mode bending (MMB) test was simulated using LS - DYNA. The b ilinear traction separation law [ 25 ] was used in CZM . Lil jedahl et al. [2 6 ] also use d a bilinear 10 traction - separation law in CZM to model adhesively bonded joints . The CZM model was inserted between two adjacent shell layers. The results show that the CZM damage model using a bilinear traction - separation law successfully predict ed the experiment results of MMB test . Turon et al. [ 27 ] also use d cohesive elements to simulate the Mode I DCB test of the unidirectional carbon fiber - reinforced composite. Turon et al. [ 27 ] compared CZM simulation results of using different interface strength values with analytical solution LEFM. The simulation results show good agreement with the LEFM solution . It was also mentioned t hat the penalty stiffness , the fracture toughness for Mode I and Mode II, the mixed mode interaction parameter and the interface strength are important parameters for cohesive damage model to accurately predict delamination . Turon et al. [ 27 ] simulate d the Mode II ENF test using cohesive elements and plane stress elements in the through - thickness direction to predict the delamination behavior of the UD carbon fiber - reinforced composite laminate. The pre - damage cohesive elements were inserted in the pre - crack area to prevent the penetration of crack interfaces. The results were wel l - agreed with the analytical solution (LEFM) despite the use of different interface strength values . In additional to CZM , decohesive elements and discrete cohesive zone model (DCZM) have also been developed to simulate interlaminar fracture toughness tests of composites. A nu merical model has been developed with decohesive elements that has zero - thickness . Camanho et al. [ 28 ] develop ed the decohesive model which could predic t the crack initiation and propagation of composite laminates under mixed mode loading. The Mode I and Mode II interlaminar fracture 11 toughness tests were si mulated using decohesive cohesive model by Camanho el at [ 28 ] . The results show good consi stency with experiment al results. The d iscrete cohesive zone model (DCZM) was developed by treating cohesive elements as spring [ 29 ] . The schematic of DCZM is shown in Figure 2. 1 of adjacent elements in DCZM . In paper [ 30 ] , a n improved DCZM was discussed . The strain field at the crack tip was integrated in the DCZM . The DCZM can be scaled as a function of element size . The geometr ic nonlinearity and crack orientation were incorporated in mode mixity of the DCZM . The Mode I simulation in [ 29 ] shows accurate prediction of the maximum load and the slope compared to test results. Figure 2. 2 load - displacement curve of stitch elements [40] Figure 2. 1 Schematic of DCZM [42] 12 2.3 Simulation of Z - pinned and Stitched Composites To improve the interlaminar strength of composite laminates, the z - pinning method is used to increase delamination resistance . The review of z - pinned composite laminates , which was done by Mouritz [ 3 1 ] , shows that the z - pinned composite laminates have significant improvement on de lamination toughness which c ould resist crack propagation. Grassi and Zhang et al. [ 4 ] ha d used a n existing micro - mechanical material model [3 2 ] to create numerical model of z - pinned composite laminate to simulate the Mode I test . A non - linear FEA was use d with Newton - Raphson method to predict crack propagation. Results show that the z - pinned composite laminates have good energy absorption ability. The load drop phenomenon, which is caused by existence of fibers in the z direction , was observed in load - extension curve . S titching is another approach to increase interla minar strength in composite laminates. Tan el at. [ 5 ] simulates Mode I DCB test with 2D FE model of Vectran - stitched composite laminate . In Chen et al. [ 33 ], 2D FE model ha d been proved to have sufficient prediction accuracy when simulating the Mode I DCB test of stitched composite laminates . The stitches in 2D FE model were created by 3 - nodes rod elements. The behavior of the rod elements is governed by load - displacement curve in Figure 2. 2 . The stitch fracture process has four procedures: interfacial debonding, slack absorption, fiber breakage and pullout friction. Integrating all of procedures into FE model results in accurate numerical prediction . [ 33 ] 13 Chapter 3. L S - D YNA Simulation In this research, the finite element model is created using LS - DYNA software. LS - DYNA is a commercial finite element program . Its explicit solver is widely used in automotive and aerospace industr ies in crash safety simulation s . The numerical model of UD and Q3D composite s are buil t for simulating M ode I double cantilever beam (DCB) test and M ode II end - notched flexure (ENF) test. The fracture toughness parameter s for UD and Q3D model s are obtained from Wente et al. [ 6 ] . Other input parameters for material model are estimated based on the Robert et al [ 14 ] . The load - extension curve of simulation is compared with experimental results for validation of the numerical model s . The FE model for the DCB and ENF specimens were built for the unidirectional (UD) composite laminat e first . The UD composite laminate has 12 layers with 0 ° fiber direction. Shell elements are used to model composite plies . To simulate delamination resistance of UD composite laminates , cohesive elements are utilized with bilinear and trilinear cohesive zone laws. Results are compared with experimental results to pick optimal cohesive zone law for Mode I and Mode II . The FE model for Q3D composite material is based on the model for the unidirectional (UD) composite laminate. To simulate Q3D composite laminates , the interface of tested Q3D specimen s are examined to identify the area of fiber bridging tow s. The fiber bridging tows increase the delamination resistance in Q3D composites . The region enclosed by red lines in Figure 3 .2 ( a) is where fiber bridging appears . For fiber bridging area, higher fracture toughness value s are used to simulate higher delamination resistance due to fiber bridging . 14 3.1 Numerical Model Setup 3.1.1 Model Dimensions and Structures for Mode I To compare the simulation results with the experimental results of Mode I interlaminar fracture toughness test, th e dimensions of the specimen and its boundary conditions in numerical model must be consistent with the specimen in experiment. The dimensions of the specimen were selected according to ASTM D 5528 [ 34 ] , as listed in Table 3.1 . The pre - crack wa s manually created by inserting a Teflon film of 25 - micron thickness in the mid - plane of the composite laminate when the composite plate wa s manufactured . The UD composite laminate has 12 layers with fiber direction s in 0°, 60° and - 60° . Figure 3.1 shows the numerical model setup for the Mode I simulation. The DCB model consists of the top and bottom composite layers and t he interface is represented by a layer of cohesive element s . The top and bottom composite layers are modeled using shell elements. The reasons of using shell elements are: 1. The components of automobiles and airplanes are often made of thin sheet s of metals or composit e materials. It is quite common to use shell element in crash safety simulation , particularly in automo tive application . 2. Crash simulation are solved by explicit method. In this method, the solution accuracy is controlled by the timestep, which is in tur n determined by the Figure 3.1 Numerical Model setup of UD Composite 15 smallest element size . If solid elements are used to model composite layers, the resulted timestep will be small. So that the run time to solve the problem will be long. Using shell elements will save computational time. The shell element for composite layer is 1.8 mm thick, with 6 integration points through the thickness . Each integration point represents one ply of c omposite laminates. The interface where the interlaminar delamination to occur is modeled using solid cohesive elements ( Figure 3.1) . The cohesive layer thickness is depending on the definition of shell reference plane. In the first case, when shell thickness reference is defined at the mid - plane, the cohesive layer thickness equals to the shell thickness which is 1.8 mm. In the second case, the shell thickness refe rence is defined at the top surface for bottom composite layer and at the bottom for top composite layer. Then the cohesive layer thickness can be reduced to 0.336 mm. The region without cohesive elements is the pre - crack area. The loading block s , which are modeled using solid elements with rigid material properties, are placed at the tip of the laminate. A prescribed displacement is assigned at the center Figure 3.2 (a) Interface of tested Q3D composite laminate (b) Cohesive layer interface of Q3D composite laminate 16 of the top loading block with strain rate of 0.04 mm/s. The top loading block is constrained t o allow only translational degree of freedom in the z direction and rotational DOF in the x direction. The bottom loading block is fixed except the rotation a bout the x - axis. For the Q3D composite laminate s , the structure of numerical model is identical to the UD composite laminate s . The dimensions of the Q3D model is listed in Table 3.1. A unique feature of the Q3D composite laminate is the bias bridging tow which connec ts adjacent layers of the composite laminate . Figure 3.2 (a) shows a post - modern DCB specimen where the bias tow bridging area s are enclosed by red lines. Two sets of Cohesive Zone Model (CZM) parameters are used in the model . The elements in red had the same CZM parameters as in the UD model. The brown cohesive elements represent ing the 60° bias bridging in Figure 3.2(b) were assigned with a second set of CZM parameters with a higher fracture toughness value . In simulation s , the brown cohesive elements produce d a higher delamination resistance for laminate interface. I t require d larger traction force to separate the top and bottom composite layers . Table 3.1 Dimensions of Mode I DCB model Length (mm) Width (mm) Thickness (mm) Pre - crack (mm) UD 184.375 25.4 3.6 52.125 Q3D 182.924 25.4 3.6 60.0 13 17 3.1.2 Model Dimensions and Structures for Mode II The ENF specimen for Mode II experiment follows ASTM D 7905 [ 35 ] . The composite laminate is placed on two bottom supporters with 4 - inch span. The distance from bottom right supporter to the end of the pre - crack ( ) is 30 mm. The upper roller in Figure 3.3 is plac ed above the composite laminate . The UD composite laminat e dimensions are shown in Table 3. 2 . It consist s of top and bottom composite layers. The c ohesive layer is located between the top and bottom composite layers. Six pairs of automatic surface to surface contacts are defined. The contact surface of the shell element is placed at a distance equal to half of the contact thickness. The default contact thickness equals the shell thickness which is 1.8 mm [ 36 ]. The shell thickness reference plane is defined at the middle in Mode II numerical model. So that the contact surface is at the mid - plane between top and bottom composite layers. If the shell thickness reference is defined at the bottom for top composite layer and at the top f or bottom composite layer, the contact surface will have initial penetration. In LS - DYNA contact card s , the master and slave selections are not necessarily following the contact sequence [ 37 ]. The optional thickness for the master and slave surface (MST, SST) are the true element thickness of the contact components. In this model, the true thickness of the shell composite layer is 1.8 mm. The cohesive layer true thickness is zero. The Rayleigh Figure 3.3 Front view of Mode II ENF model setup 18 damping coefficient for composite layer is set to be 0.05 . The strain rate of Mode II ENF simulation is 0.005 mm/s. Table 3. 2 Dimensions of Mode II ENF model The dimensions of the Q3D composite laminate are listed in Table 3.2. T wo sets of CZM parameters are used . The cohesive layer setup is shown in Figure 3.4 . A set of CZM parameters with higher fracture toughness values is assigned to brown cohesive elements to simulate higher delamination resistance of bias bridging tows. Length (mm) Width (mm) Thickness (mm) Pre - crack (mm) UD 156.25 20 3.6 53.125 Q3D 143.75 20 3.6 53.125 Figure 3.4 Cohesive layer setup of Q3D composite laminate for Mode II 19 3.2 Cohesive Zone Law There are four typical traction - separation law s to govern the behavior of cohesive elements. Table A.1 de monstrate s the traction - separation curves for different material cards in LS - DYNA . In this research, two CZM laws, i.e . MAT_138 and MAT_185 are used to simulate Mode I and Mode II interlaminar fracture toughness tests . 3.2.1 MAT_138 Bilinear Cohesive Zone Law MAT_138 in LS - DYNA is governed by a bilinear traction - separation law with quadratic mix ed mode delamination criterion and damage formulation to simulate delamination behavior of cohesive elements. In order to use MAT_138, cohesive element formulation 20 is defined in SECTION_SOLID of cohesiv e elements [ 25 ] . The 8 - nodes cohesive element is used to generate cohesive layer. The number of integration points for a cohesive element to be deleted is set to 4. The traction s are defined a s mid dle point between each pair of nodes (1 - 5, 2 - 6, 3 - 7, 4 - 8) on a n 8 - Figure 3.5 (a) MAT_138 Bilinear Traction - separation curve (b) MAT_185 Trilinear Traction - separation curve [25] 20 nodes cohesive element in Figure 3.7 . [ 25 ] A cohesive element will be deleted if 4 integration points of the cohesive element reach their peak traction . The Bilinear Cohesive Zone Law complies mixed mode criterion [ 12 ]. The total mixed mode displacement is governed by Eq 3. 1 : ( 3 . 1 ) Figure 3. 6 illustrates the displacement inside the single cohesive element . The displacement in 3 direction represents the displacement in normal direction for mode I. The displacement in 1 and 2 directions are the displacement in tangential direction for mode II. The damage initia tion displacement for mode I is calculated by dividing peak traction in normal direction Figure 3.6 Displacement in cohesive element 1 2 3 4 5 6 7 8 Integration point Figure 3.7 Integration points on a cohesive element 21 (T) by normal to cohesive element stiffness (EN). For mode II damage initia tion displacement is obtained by dividing peak traction in tangential direction ( S) by in - plane stiffness of cohesive element (ET) . The damage initiation displacement is shown in Eq 3. 2 : ( 3. 2 ) Where . The total failure mixed - mode displacement is provided in Eq 3. 3 , where the mixed mode criterion exponential (XMU) is set to be 1. ( 3. 3 ) Figure 3. 5 (a) provides the traction - separation curve for the bilinear cohesive zone law. The area under the triangle represents the energy release rate which is also known as the fracture toughness for mode I , and for mode II. The traction within the cohesive elements increases as the separation of cohesive elements increases until reaching peak traction. The softening process starts after the peak load is reached. Eventually the traction drops to zero when failure of cohesive elements occurs. T he expe rimentally measured mode I and mode II fracture toughness value s for UD and Q3D composite s by Wente et al. [ 6 ] are listed in Table 3. 3 and Table 3.4 . The parameters in the CZM laws were determined using the average fracture toughness for Mode I and non - pre - crack (NPC) fracture toughness values for Mode II . These values are listed in Table 3.5. 22 For Q3D, a second set of CZM parameters with higher Mode I and II fracture toughness values were introduced to model the bridging bias tows. Since the DCB and ENF specimens had dif ferent width, the area ratio of the bridging tows on the fractured specimens were different. Therefore, the second set CZM parameters were slightly different for the DCB and ENF specimens. For the DCB s pecimen , t he volume fraction of bias tow s is 19.55% which calculated using Eq 3. 4 . The volume UD 2DW Q3D NPC 644.5 502.2 613.8 PC 559.5 696.1 824.3 fraction of matrix ( ) is 8 0 . 4 5 %. The average Mode I fracture toughness value of Q3D composite laminates that measured in experiment ( ) is . The Mode I fracture toughness for the bias tow elements ( ) can be calculated using Eq 3. 5 . Table 3.6 presents the fracture toughness values for DCB simulations . UD 2DW Q3D 563.6 564.1 466.3 732.8 730.4 884.7 Energy 3.345 3.756 4.994 Table 3.3 Mode I fracture toughness values measured in DCB experiments [6] Table 3.4 Mode II fracture toughness values measured in ENF experiments [6] 23 In ENF simulation of Q3D composites, the width of the numerical model is shorter than the width of the Q3D model in Mode I. The volume fraction of bias tow cohesive elements is 20.7% . The Mode I fracture toughness value is then re - calculated to be . Table 3.7 present the fracture toughness values for ENF simulations Table 3. 5 Fracture toughness value s of UD composite in Mode I and II Composite Mode I Fracture Toughness Value Mode II Fracture Toughness Value Unidirectional (UD) Table 3. 6 Fracture toughness values of Q3D composite for Mode I Q3D composite Mode I Fracture Toughness Value Mode II Fracture Toughness Value Matrix cohesive elements Bias tow cohesive elements 24 Table 3.7 Fracture toughness values of Q3D composite for Mode II ( 3. 4 ) ( 3. 5 ) ( 3. 6 ) ( 3. 7) The relationship between the ultimate displacement at failure (u) and displacement at the peak load (L) for the bilinear traction - separation law is shown in Eq 3. 6 and Eq 3. 7 . The ultimate displacement at failure (u) must be larger than the displacement at peak load (L). The and Q3D composite Mode I Fracture Toughness Value Mode II Fracture Toughness Value Matrix cohesive elements Bias tow cohesive elements 25 values in Eq 3. 6 and 3. 7 are give n , the peak traction s in the normal (T) and tangential (S) directions are estimated. The normal to cohesive element stiffness (EN) and the in - plane stiffness of cohesive element (ET) can be calculated using Eq 3.6 and 3.7 . The relationship between the fracture toughness values and the peak tractions is indicated in Eq 3. 8 and 3. 9 . The ultimate displacement in the normal direction (UND) and the ultimate displacement in the tangential direction (UTD) are calculated by using Eq 3. 8 and 3. 9 . ( 3. 8) ( 3. 9 ) 3.2.2 MAT_18 5 Trilinear Cohesive Zone Law MAT_185 in LS_DYNA follows a trilinear cohesive law . It was develo ped by Tvergaard and Hutchinson in 1992 [ 20 ] to use with cohesive element formulations [ 25 ]. The area under the trapezoidal traction - separation curve in Figure 3.5 ( b ) represents the fracture toughness value . The peak traction s in the trilinear cohesive law and the bilinear cohesive law ha ve equivalent value s . The scaled distance for failure ( ) is 1. The scaled distance to the peak traction ( ) is set to 0.2 and the scaled distance to beginning of softening ( ) is set to 0.5. The values of and can be adjusted. In DCB and ENF simulation of UD composites using Trilinear CZM , the fracture toughness parameters in the Bilinear CZM are the one used here to calculate the scaled length in the normal (NLS) and tangential directions (TLS) using Eq 3.10 and 3.11. The peak traction is 0.009 26 GPa. The NLS and TLS values for UD composites in DCB and ENF simulation s are l isted in Table 3.8 and Table 3.9 . ( 3. 10) ( 3. 1 1 ) Table 3.8 Scaled Length parameters for UD composites in DCB and ENF simulation s Composite Scaled Length in Normal Direction (NLS) Scaled Length in Tangential Direction (TLS) Unidirectional (UD) Table 3.9 Scaled Length parameters for Q3D composites in ENF simulation Q3D composite Scaled Length in Normal Direction (NLS) Scaled Length in Tangential Direction (TLS) Matrix cohesive elements Bias tow cohesive elements 27 In ENF simulation of Q3D composites, two sets of fracture toughness parameters in the Bilinear CZM of Q3D composites are used to calculate NLS and TLS using Eq 3.10 and 3.11 for bias tow cohesive elements and matrix cohesive elements . Look ing at the cohesive element model i n Figure 3. 6 , the dimensionless separation of the element ( ) combines the interaction between the displacement in Mode I ( ) and Mode II ( ) in Eq 3.12 . The trilinear traction - separation law used for stress calculation is show in Eq 3.13 . The normal and tangential component s of traction vector are given by Eq 3.14 [ 25 ]. (3.12) (3.13) (3.14) 28 3.3 Material P roperties Th e A - 42, 12k carbon fiber , which is manufactured by D ow A ksa , is used to make our UD and Q3D preform s . The properties of the fiber are listed in Table 3. 10 . SC - 15 Epoxy was used as matrix resin. The properties of SC - 15 Epoxy [ 38 ] are given in Table 3. 11 . The composite panels for interlaminar fracture testing were made with 10 layers of commercial triaxial braided composite as the outer layers and two layers of handmade composite preform at the interface , as shown in Figure 3.8 . The composites were manufactured by the VARTM process. The fiber volume fraction of the composite was estimated by acid digestion [ 6 ]. The Rule of Mixture ( Eq 3.15) is then utilized to calculate composite density . The composite density is . As the mechanical properties of the UD and Q3D composites are unknown , they were assumed to be close to the commercial triaxial braided composite. L iterature data on T700/PR520 UD composite [ 14 ] were used to compute the properties of the triaxial braided composite using Micromechanics approach . Table 3. 12 presents the estimated m echanical properties for the composite layers . Figure 3.8 Architecture of the Q3D composite specimen for DCB and ENF tests [6] 29 Table 3. 10 Carbon fiber properties Carbon Fiber Tensile Strength ( ) Tensile Modulus ( ) Strain (%) Density ( ) A - 42, 12k 4.2 240 1.8 1.78 Table 3. 11 SC - 15 Epoxy properties Matrix Density ( ) Compressive Strength (G ) Viscosity (centipoises) SC - 15 1.09 0.088 300 (3.15) Table 3. 12 Composite mechanical properties Mechanical Properties Unit Carbon Fiber Composite Density Kg/ m 3 1 460 Axial modulus GPa 38 Transverse modulus GPa 32.7 In - plane shear modulus GPa 8.96 In - plane P o _ 0. 3 Axial tensile failure strain _ 0.0216 Axial c ompressive failure strain _ 0.018 30 Table 3. 12 Transverse tensile failure strain _ 0.0168 Transverse compressive failure strain _ 0.011 In - plane shear failure strain _ 0.024 Axial tensile stress at failure GPa 1.04 4 Axial compressive stress at failure GPa 0.377 Transverse tensile stress at failure GPa 0.362 Transverse compressive stress at failure GPa 0.345 In - plane shear stress at failure GPa 0.307 31 Chapter 4. Simulation Results and Discussion 4.1 Mode I Results 4.1.1 UD C omposite with Bilinear CZM Two different thickness es of cohesive layers are investigated here. For the first case, t he shell thickness refer ence plane is defined the middle. The cohesive layer thickness is 1.8 mm which is relatively thick compared to the second case , where the shell thickness reference plane is defined at the top for bottom shell and at the bottom for top shell and the cohesive layer is relatively thin . The load - extension curve for the unidirectional composite laminate with thick bilinear CZM is compared with the experimental load - extension curve in Figure 4.1. The green curve indicates the simulation result and the blue and red curve s illustrate the experimental result. Before load reaches its peak value , the load - extension curve generally displays a linear elastic behavior . However, the Figure 4.1 Comparison of the experimental load - extension responses of UD DCB specimen under Mode I with the prediction with thick cohesive element with a Bilinear CZM law. 0 10 20 30 40 50 60 0 10 20 30 40 50 60 70 80 90 Load (N) Extension (mm) Load - Extension Curve Experiment 1 Experiment 2 Simulation 32 non - linearity is also observed at approximately 1 2 mm . A s imilar trend is also observed in experiment [ 39 ] . T he stiffness and crack growth results of the simulation agree well with the experimental load - extension curve. The numerical pre diction of the peak load is higher than the experiment. After the peak load, the softening part of the load - extension curve has the same trend compar ed with the e xperimental result . The numerical simulation is unable to precisely predict the instability o f crack propagation in the experiment due to delamination resistance of matrix [ 40 ]. The oscillation of the DCB also observed in simulation after the crack starts to propagate. The result of thin Bilinear CZM is presented in Figure 4.2 . The load - extension curve generally has the same elastic load increase to initiate the crack propagation compared with experimental results . The load decrease after the peak has the same trend as the experiment. The peak load is over - 0 10 20 30 40 50 60 0 10 20 30 40 50 60 70 80 90 Load (N) Extenison (mm) Load - Extension Curve Experiment 1 Experiment 2 Simulation Figure 4.2 Comparison of the experimental load - extension responses of UD DCB specimen under Mode I with the prediction with thin cohesive element with a Bilinear CZM law. 33 predicted . In the simulation, the FE model with a thin cohesive layer in the DCB tends to be more stable compared with the one with a thick cohesive layer . Figure 4.3 Comparison of Bilinear CZM law with the same but different traction values. parameters Figure 4.4 Comparison of the load - extension curves by simulations with the two CZMs 0 10 20 30 40 50 60 0 10 20 30 40 50 60 70 80 90 Load (N) Extension (mm) Load - Extension Curve CZM 1 CZM 2 0 0.002 0.004 0.006 0.008 0.01 0 0.1 0.2 Traction (GPa) Displacement (mm) CZM 1 CZM 2 34 Two different bilinear traction - displacement curves are compared in Figure 4. 3 . The fracture toughness values that are used for Mode I simulation remain unchanged for Cohesive Zone Model 1 and 2. CZM 1 and 2 ha ve t he same linear elastic slope of the traction - displacement curve . The peak load of CZM 2 is higher than CZM 1 . To keep the M ode I fracture toughness value unchanged for CZM 2, the ultimate displacement at failure (u) needs to be reduced for CZM 2 . Figure 4. 4 indicates that load - extension curves for Bilinear CZM 1 and CZM 2 are identical in stiffness, maximum load and crack growth . It is interesting to not e that , if the Mode I fracture toughness value stays the s ame for two different bilinear traction - separation law s , the numerical prediction s with two sets of C ZM parameters have no major differences. It is important to use the correct Mode I fracture toughness value in the C ZM in simulati ng interlaminar fracture toughness test s . Eq 4.1 is used to calculate the Mode I fracture toughness value from DCB experiment, where P is the measured load in Mode I experiment, is the load point displacement , b is the specimen width, a is the delamination length and is the compliance . From the M ode I experiment, two fracture toughness values are obtained: the initial value and the averaged value. In Table 3.3, the initial fracture toughness value is calculated using the measured load to initiate the crack propagation. The Figure 4.5 shows th e values that - [ 41 - 42 ] . The average fracture toughness value is higher than the initial fracture toughness value because the fiber bridging has occurred in the UD specimens in Mode I expe riment. In Figure 4.6, fiber bridging is clearly observed in the red circle. 35 Figure 4.6 Mode I Experiment picture of UD composite laminate 0 100 200 300 400 500 600 700 800 900 1000 0 50 100 150 200 G I [J/m 2 ] Crack Length [mm] UD R Curves UD1 UD3 Initial G IC Figure 4.5 The experimental R curves for UD composites 36 ( 4.1 ) The initial fracture toughness value is used in Bilinear CZM parameters, the numerical prediction result is compared with the load - extension curve result obtained by using average fracture toughness value. Figure 4.7 shows the load - extension curve results. The green curve indicates the numer ical prediction result of Bilinear CZM with the average fracture toughness value . The orange curve represents the prediction result of Bilinear CZM with the initial fracture toughness value . Th e elastic and softening behavior of the load - extens ion curves are the same for Bilinear CZM with initial and average . However, the numerical result of using a higher fracture toughness value gives a higher peak load prediction . The same phenomenon was observed in Turon et. al [ 27 ]. The CZM using the initial fracture toughness value predicts better peak load than CZM using average fracture toughness value . The CZM with average fracture 0 10 20 30 40 50 60 0 10 20 30 40 50 60 70 80 90 Load (N) Extension (mm) Load - Extension Curve Experiment 1 Experiment 2 GIC=733 GIC=564 Figure 4.7 Comparison of Bilinear CZM with initial and Bilinear CZM with average 37 toughness value overshoots the peak load but has better prediction of softening behav ior. To precisely predict the crack onset and propagation in UD composites under Mode I condition , incorporating both the initial and average fracture toughness values in CZMs is recommended. 4.1.2 UD C omposite with Trilinear CZM The UD DCB experiment s were also simulated with a Trilinear C ZM law . Thick and thin cohesive layers are also investigated with Trilinear CZM. Figure 4.8 demonstrates the load - extension curve of Trilinear CZM with thick cohesive layer. The elastic load increase to initiate the crack propagation is consistent with the experiment results. The peak load is over - predicted . The DCB tend to be unstable as the crack starts to propagate. Figure 4.8 Load - extension cur ve of Trilinear CZM with thick cohesive layer 0 10 20 30 40 50 60 0 10 20 30 40 50 60 70 80 90 Load (N) Extension (mm) Load - Extension Curve Experiment 1 Experiment 2 Simulation 38 Looking at the load - extension curve of Trilinear CZM with thin cohesive layer in Figure 4.9 , the elastic part of the load - extension curve agre es with the experiment. The peak load is over - predicted. After the peak load, the instability is observed. The traction - separation curve is compared with the Bilinear C ZM law in Figure 4. 10 . The fracture toughness value s are the same for both CZM law s. The elastic slope of traction - separation curve for the Bilinear CZM to reach the peak traction is the same as the Trilinear CZM . With the same maximum traction, Trilinear CZM has a smaller ultimate displacement at failure (u) that eventually leads to a steeper slope of the softening process. The simulation results are shown in Figure 4. 11 . The orange and green curve s are corresponding to Bilinear and Trilinear CZM s . The elastic behavior of Trilinear CZM has the same trend as that of the Bilinear CZM. The non - linearity is not observed in load - extension curve s of both Bilinear and Trilinear CZM s . The peak load of the 0 10 20 30 40 50 60 0 10 20 30 40 50 60 70 80 90 Load (N) Extension (mm) Load - Extension Curve Experiment 1 Experiment 2 Simulation Figure 4.9 Load - extension curve of Trilinear CZM with thin cohesive layer 39 Trilinear CZM agrees with that of the Bilinear CZM. The instability of Trilinear CZM is observed in load - extension curve after the peak load . 0 0.002 0.004 0.006 0.008 0.01 0 0.1 0.2 Traction (GPa) Displacement (mm) Trilinear CZM Bilinear CZM Figure 4.10 Comparison of traction - separation curves for Bilinear CZM and Trilinear CZM 0 10 20 30 40 50 60 0 10 20 30 40 50 60 70 80 90 Load (N) Extension (mm) Load - Extension Curve Experiment 1 Experiment 2 Bilinear CZM Trilinear CZM Figure 4.11 Load - extension curve comparison of Bilinear CZM and Trilinear CZM 40 To find out the cause of the instability in Trilinear CZM, the cohesive element deletion in cohesive layer is investigated . Figure 4. 1 2 illustrates the cohesive element deletion inside the cohesive layer for the Bilinear CZ M . It takes 17ms for one row of red cohesive elements to be deleted and totally Figure 4.12 Cohesive element deletion of Bilinear CZM Figure 4.13 Cohesive element deletion of Trilinear CZM 41 66ms to delete two rows of cohesive elements. On the contrary , the Trilinear CZM has sudden cohesive element deletion in Figure 4.1 3 . Over the time interval of 4 ms, one row of cohesive elements has been deleted because of the element failure. In 53ms, two rows of cohesive elements have been deleted for the Trilinear CZM. The relatively large number of cohesive elements deletion in short amount of time leads to a sudden internal energy transformation which result s in instability for Trilinear CZM. 4.1.3 Q3D Composite The numerical prediction results for the UD composite reveals that the Bilinear C ZM is more stable than the Trilinear Cohesive Zone Model in crack propagation . The Bilinear CZM is more viable to simulate Mode I interlaminar fracture toughness test. Higher fracture toughness value s are assigned to 60° bias tow cohesive elements which cause s several load - drops in load - extension curve . This load - drop phenomenon is also refer red - in literature paper s [ 41 - 42 ]. Two Figure 4.14 Load - Extension curve of Q3D composite laminate with thick cohesive layer 0 10 20 30 40 50 60 0 10 20 30 40 50 60 70 80 90 Load (N) Extension (mm) Load - Extension Curve Experiment 1 Experiment 2 Simulation 42 cases for cohesive layer thickness are also investigated here. Look ing at numerical prediction of Q3D composite laminate with thick cohesive layer in Figure 4.1 4 , t he increase of the elastic load before crack initiation is consistent with the load - extension curve of experiment 2 . The numerical prediction of peak load - experimental result. Each peak in the experimental load - extension curve represents a fiber tow bridging event . S everal load drops are observed in the experimental load - extension curve. This phenomenon occurs when a crack starts to propagate and reaches bias bridging tows . As the bias tow has a higher delamination resistance , it requires a higher load to break the fiber tow for crack growth to continue . The general trend of crack growth at the interface agrees with the experiment al load - extension curve. However, the numerical simulation is not sufficiently accurate to mimic stick - slip in the load - extension curve after the peak load . 0 10 20 30 40 50 60 0 10 20 30 40 50 60 70 80 90 Load (N) Extension (mm) Load - Extension Curve Experiment 1 Experiment 2 Simulation Figure 4.15 Load - Extension curve of Q3D composite laminate with thin cohesive layer 43 The load - extension curve of the Q3D composite laminate with thin cohesive layer show consistency of elastic process compared with experiment 2. The peak load is slightly over - - ly predicted after the peak load. 4.2 Mode II Results 4.2.1 UD Composite with Bilinear CZM The numerical prediction of the UD compo site with the Bilinear C ZM is investigated here . The green curve in Figure 4.1 6 represents the numerical prediction of bilinear CZM for the UD composite laminate . The e lastic part of the load - extension curve before peak load is consistent with experimental results . The model with the Bilinear CZM precisely predict the load to initia te the crack. However, the peak load is under - predict ed . The stiffness of the numerical pr ediction is Figure 4.16 Comparison of the ENF experimental Load - deflection curve of UD composite with simulation using the Bilinear CZM. 0 50 100 150 200 250 300 350 0 1 2 3 4 5 6 Load (N) Extension (mm) Load - Extension Curve Experiment 1 Experiment 2 Simulation 44 reduced because of the crack growth . The sudden load drop after the peak traction of experimental curve does not appear in numerical prediction. To investigate whether the increase of in - plane Figure 4.17 CZM comparison with different stiffness 0 50 100 150 200 250 300 350 0 1 2 3 4 5 6 7 8 Load (N) Extension (mm) Load - Extension Curve Bilinear CZM Bilinear CZM 2 Figure 4.18 Load - deflection curve comparison with two CZMs 45 stiffness in CZM will stiffen the bending behavior of the composite beams , CZM 2 is introduced with an increased in - plane stiffness (ET). The traction separation curve is presented in Figure 4.1 7 . The stiffness (ET) change results in the sh ift of the displacement at the peak traction (L). However, t he load - extension curves generated in simulations with two CZMs are overlapping each other , Figure 4.1 8 . This shows that t h e increases of in - plane stiffness (ET) does not affect the numerical prediction results. 4.2.2 UD Composite with Trilinear CZM Since the Bilinear CZM demonstrates significant stiffness reduc tion in the elastic part of the load - extension curve and no discrepancy due to load drop after the peak traction, the Trilinear CZM is investigated to see if it can produce different results. The Trilinear CZM result is provided in Figure 4.1 9 . The red and blue curve s represent the UD experiment result s . The load to initial the Figure 4.19 Comparison of ENF experimental Load - deflection curve with simulation with a Trilinear CZM. 0 50 100 150 200 250 300 350 0 1 2 3 4 5 6 Load (N) Extension (mm ) Load - Extension Curve Experiment 1 Experiment 2 Simulation 46 crack agrees well with experiment. The decrease of bending stiffness due to crack growth is small compar ed with the Bilinear CZM. The predicted peak load is slightly smaller than the experimental peak load. The discrepancy after the peak load is observed in Figure 4.1 9 . I nstability after the peak 0 0.002 0.004 0.006 0.008 0.01 0 0.1 0.2 Traction (GPa) Displacement (mm) Trilinear CZM Bilinear CZM Figure 4.20 Traction - separation curves comparison for Bilinear and Trilinear CZMs 0 50 100 150 200 250 300 350 0 1 2 3 4 5 6 7 8 Load (N) Extension (mm) Load - Extension Curve Experiment 1 Experiment 2 Bilinear CZM Trilinear CZM Figure 4.21 Load - extension curves comparison for Bilinear and Trilinear CZMs 47 load for Trilinear CZM is also detected in ENF simulation due to sudden cohesive elements deletion . Overall, the Trilinear CZM provides relative ly accurate numerical prediction for the delamination behavior of UD composite s . The traction - separation curves for Bilinear and Trilinear CZMs are compared in Figure 4.20, the elastic load increasing to reach the peak traction is the same for both CZMs. The fracture toughness value for both CZMs are the same. The slope of the softening process is steeper for Trilinear CZM. The load - extension curve in Figure 4.21 shows that the elastic load increase and softening behavior of both CZMs are consistent with each other . The Trilinear CZM predicts a higher peak load for crack to propagates and has less stiffness reduction compared with the Bilinear CZM. Figure 4.22 Load - extension curve of Bilinear CZM for Mode II 0 50 100 150 200 250 300 350 400 450 0 1 2 3 4 5 6 7 Load (N) Extension (mm) Load - Extension Curve Experiment 1 Experiment 2 Simulation 48 4.2.3 Q3D Composite Contrary to Mode I DCB test simulation s , the model with a Trilinear CZM generated accurate numerical prediction of UD composite under Mode II ENF test . To numerically predict the delamination behavior of Q3D composite, the model with Bilinear CZM was tested first to check the prediction accuracy . Higher fracture toughness values are used for bias tow cohesive elements (brown elements) as shown in Figure 3.4. The load - deflection curve predicted with Bilinear CZM is shown in Figure 4. 22 . The result reveal ed an under - prediction of the peak load coma pared with experimental value ( green curve) . The Q3D composite became les s stiff after the crack starts to propagate. The trend in discrepency between the experimental curve and predict ed response in Q3D ENF simulations with the Bili nea r CZM is consistant with UD ENF simulations. Because the B ilinear CZM fail ed to accurately predict delamination behavior of Q3D composite laminate, Trlinear CZM wa s then tested for cohesive layer. A more accurate numerical prediction Figure 4.23 Load - extension curve of Trilinear CZM for Mode II 0 50 100 150 200 250 300 350 400 450 0 1 2 3 4 5 6 Load (N) Extension (mm) Load - Extension Curve Experiment 1 Experiment 2 Simulation 49 is obtained , as shown in Figure 4. 23 . The bending stiffness reduction is improved compar ed with Bilinear CZM prediction results in Figure 4. 22 . The under - prediction of peak load wa s still seen for Trilinear CZM. The discrepancy of the load drop is more distinct comparing to experimental result. In general, the Trilinear CZM produces more accurate numerical prediction of delamiantion behavior for Q3D composite laminate than the Bilienar CZM does. 50 Chapter 5. CONCLUSION 5.1 Summary and Conclusion In this research, the delamination behavior of UD and Q3D composites are simulated using Cohesive Zone Model (CZM). The interface in UD composite was simulated with one set of CZM parameters. To simulate Q3D composite, two sets of CZM parameters are used. A set of CZM para meter with higher fracture toughness value is assigned to bias tow cohesive elements which - Bilinear CZM and Trilinear CZM are used to simulate Mode I and Mode II experiments of UD and Q3D composite s . The n umerical simulation results of Bilinear CZM and Trilinear CZM are compared with the experimental result s, Table 5.1 and 5.2 shows that Bilinear CZM has better accuracy when predict ing the peak load of UD and Q3D composites under Mode I condition . The Trilinear CZM gives an uns t able softening process of load - extension curve because of sudden cohesive elements deletion in cohesive layer. The load to initiate crack propagation is consistent with experiment. Howev er, numerical simulation s could not precisely predict the instability of the crack propagation. The linearity is observed in the elastic part of the load - extension curve. The two different definitions of shell thickness reference plane result in thick and thin cohesive layer thickness. The results in Table 5.1 - 5.4 shows that a thick cohesive layer better predicts the peak load and displacement at the peak load for both UD and Q3D composites. To accurately simulate Mode I experiment of UD and Q3D composites, the Bilinear CZM with a thick cohesive layer would be the optimal method. In this thesis, t he effect of the facture toughness value on the numerical simulation results is also investigated , the initial fracture toughness value is used in 51 Bilinear CZM for the UD composite laminate under Mode I condition. The results show that using the initial fracture toughness value better predicts the load to initiate the crack propagation. On th e contrary, using the average fracture toughness value gives a better numerical prediction of crack propagation. A h igher fracture toughness value results in a higher peak load of the load - extension curve. The optimal method to numerically predict t he delamination behavior of UD composites is to use both and in CZM. In the Mode II simulation of UD and Q3D composite laminates, both the Bilinear and Trilinear CZMs under - predict the peak load. However, the Trilinear CZM clearly predicts a h igher peak load for crack to propagates. The elastic and softening behavior for both CZMs are consistent. Whereas, t he stiffness reduction in elastic load increasing process is less for Trilinear CZM. The changing stiffness of CZM does not have significantly influence on simulation results. The instability and discrepancy of the load drop in softening process of load - extension is spotted for Trilinear CZM. In general , Table 5.1 - 5.4 show that the Trilinear C ZM, when simulating UD and Q3D composites under Mode I and Mode II conditions, demonstrates a better prediction accuracy of the peak load and displacement at the peak load compared with Bilinear CZM . 52 Table 5.1 Comparison of the predicted and measured peak load with Bilinear CZM Mode I with Bilinear CZM Mode II with Bilinear CZM UD (Thick) UD (Thin) Q3D (Thick) Q3D (Thin) UD Q3D Predicted peak load (N) 49.3 51.6 46.5 48.38 251 25 2 Measured peak load (N) 46.598 46.598 45.23 45.23 317.25 358.9 Error (%) 5.5 9.7 2.7 6.5 - 2 6.4 - 42.4 Table 5. 2 Comparison of the predicted and measured peak load with Trilinear CZM Mode I with Trilinear CZM Mode II with Trilinear CZM UD (Thick) UD (Thin) UD Q3D Predicted peak load (N) 52 52.7 28 7 31 1 Measured peak load (N) 46.598 46.598 317.25 358.9 Error (%) 10.4 11.6 - 10.5 - 1 5.4 53 Table 5. 3 C omparison of the predicted and measured displacement at the peak load with Bilinear CZM Mode I with Bilinear CZM Mode II with Bilinear CZM UD (Thick) UD (Thin) Q3D (Thick) Q3D (Thin) UD Q3D Predicted displacement at peak load (mm) 13.44 13.2 29.2 28.4 3.32 3.35 Measured displacement at peak load (mm) 19.41 19.41 29.75 29.75 3.57 4.27 Error (%) - 44.4 - 47.0 - 1.9 - 4.8 - 7.5 - 27.5 Table 5. 4 C omparison of the predicted and measured displacement at the peak load with Trilinear CZM Mode I with Trilinear CZM Mode II with Trilinear CZM UD (Thick) UD (Thin) UD Q3D Predicted displacement at peak load ( mm ) 14.2 14 3.49 3.74 Measured displacement at peak load ( mm ) 19.41 19.41 3.57 4.27 Error (%) - 36.7 - 38.6 - 2.3 - 1 4.2 54 5.2 Future Work To further improve the numerical model, the following works are suggested: (1) v alidate the numerical model by simulating Mode III and mixed - mode fracture toughness test , a nd (2) examine the numerical model by simulating impact test. In addition, a systemati c characterization of the mechanical properties of Q3D composite is recommended. This would give accurate i nput of composite properties in LS - DYNA which could improve the numerical prediction accuracy on softening process of load - extension curve . Accurately i dentifying the fiber bridging area at the interface of Q3D composite laminates using microscope could also improve the accuracy of prediction on softening process of load - extension curve. In the Mode II simulation of Q3D composites , the instability o f softening process is observed for Trilinear CZM . Other cohesive zone law could implement in cohesive zone model such as arbitrary normalized traction - separation law. Instead of having discontinuities of traction - separation law for Bilinear and Trilinear CZMs [ 43 ] , arbitrary normalized traction - separation law is continuous throughout the curve which might reduce the instability of the softening process. 55 APPENDIX 56 Table A.1 Different Cohesive Zone Laws MAT_138 MAT_185 MAT_186 MAT_041_050 User defined material model, toughness and peak can be defined with different values 57 BIBLIOGRAPHY 58 BIBLIOGRAPHY [1] Rosario, Kirit, and Dahsin - Three - Dimensional Composites - with vol. 44, no. 25, Dec. 2010, pp. 2953 2973. [2 ] A.P. Mouritz, Compression properties of z - pinned composite laminates, Composites Science and Technology, Volume 67, Issues 15 16, 2007, Pages 3110 - 3120, ISSN 0266 - 3538, https://doi.org/10.1016/j.compscitech.2007.04.017. [3] Dickinson LC, Farley GL, Hinders MK. Prediction of effective three - dimensional elastic constants of translaminar reinforced composites. J Comp Materials 1999; 33:1002 29. [4] Grassi, M & Zhang, Xiang & Meo, Michele. (2002). Prediction of stiffness and stresses in z - fibre reinforced composite laminates. Composites Pa rt A: Applied Science and Manufacturing. 33. 1653 - 1664. 10.1016/S1359 - 835X (02)00137 - 9. [5] Tan, K. T., Watanabe, N., Sano, M., Iwahori, Y., & Hoshi, H. (2010). Interlaminar Fracture Toughness of Vectran - stitched Composites - Experimental and Computational Analysis. Journal of Composite Materials , 44 (26), 3203 3229. https://doi.org/10.1177/0021998310369581 [6] Tony Wente, Xinyu Mao, Danielle Zeng, Jeff Dahl, Xinran Xiao, Interlaminar fracture toughness of a quasi 3D braided composite, ASC conference paper, 2019 [7] W. Zhou, Peridynamic modeling and impact testing of dynamic damage, fracture, and failure process in fiber - reinforced composite materials. PhD thesis, Michigan State University, 2018. [8] A.B. Nair, R. Joseph, Eco - friendly bio - composites using nat ural rubber (NR) matrices and natural fiber reinforcements, Editor(s): Shinzo Kohjiya, Yuko Ikeda, Chemistry, Manufacture and Applications of Natural Rubber, Woodhead Publishing, 2014, Pages 249 - 283, ISBN 9780857096838. [9] H Mao, S Mahadevan, Fatigue dama ge modelling of composite materials, Composite Structures, Volume 58, Issue 4, 2002, Pages 405 - 410, ISSN 0263 - 8223. 59 [10] Philippe H Geubelle, Jeffrey S Baylor, Impact - induced delamination of composites: a 2D simulation, Composites Part B: Engineering, Vol ume 29, Issue 5, 1998, Pages 589 - 602, ISSN 1359 - 8368. [11] Rosario, Kirit & Liu, Dahsin. (2010). Assessment of Quasi - Three - Dimensional Composites - with Discussions on Fiber Straining and Weaving Effectiveness. Journal of Composite Materials - J COMPOS MAT ER. 44. 2953 - 2973. [12] T - A Robust Cohesive Zone Model for Cyclic Loading . [13] Dominik Bender, Jens Schuster, Dirk Heider, Flow rate control during vacuum - assisted resin transfer molding (VARTM) processing, Composites Science and Technology, Volume 66, Issue 13, 2006, Pages 2265 - 2271, ISSN 0266 - 3538. [14] Roberts, Gary, D.; Goldberg, Robert, K.; Binienda, Wieslaw, K.; Arnold, William, A.; Littell, Justin, D.; Kohlman, Lee, [15] Liu Sheng, Quasi - impact damage initiation and growth of thick - section and toughened composite materials, International Journal of Solids and St ructures, Volume 31, Issue 22, 1994, Pages 3079 - 3098, ISSN 0020 - 7683, https://doi.org/10.1016/0020 - 7683(94)90042 - 6. [16] Zou, Z., Reid, S. R., Li, S., & Soden, P. D. (2002). Modelling Interlaminar and Intralaminar Damage in Filament - Wound Pipes under Quasi - Static Indentation. Journal of Composite Materials , 36 (4), 477 499. https://doi.org/10.1177/0021998302036004539 . NASA/CR - 2002 - 211628, ICASE Report No. 200 2 - 10, April 2012. [18] E.F. Rybicki, M.F. Kanninen, A finite element calculation of stress intensity factors by a modified crack closure integral, Engineering Fracture Mechanics, Volume 9, Issue 4, 1977, Pages 931 - 938, ISSN 0013 - 7944, https://doi.org/10.10 16/0013 - 7944(77)90013 - 3. [19] I.S. Raju, Calculation of strain - energy release rates with higher order and singular finite elements, Engineering Fracture Mechanics, Volume 28, Issue 3, 1987, Pages 251 - 274, ISSN 0013 - 7944, https://doi.org/10.1016/0013 - 7944(8 7)90220 - 7. 60 [20] Viggo Tvergaard, John W. Hutchinson, The relation between crack growth resistance and fracture process parameters in elastic - plastic solids, Journal of the Mechanics and Physics of Solids, Volume 40, Issue 6, 1992, Pages 1377 - 1397, ISSN 002 2 - 5096. [21] Needleman AA. A Continuum Model for Void Nucleation by Inclusion Debonding. ASME. J. Appl. Mech. 1987;54(3):525 - 531. doi:10.1115/1.3173064. [22] G.T. Camacho, M. Ortiz, Computational modelling of impact damage in brittle materials, Internation al Journal of Solids and Structures, Volume 33, Issues 20 22, 1996, Pages 2899 - 2938, ISSN 0020 - 7683, https://doi.org/10.1016/0020 - 7683(95)00255 - 3. [23] Ku - Hyun Jung, Do - Hyoung Kim, Hee - June Kim, Seong - Hyun Park, Kyung - Young Jhang, Hak - Sung Kim, Finite element analysis of a low - velocity impact test for glass fiber - reinforced polypropylene composites considering mixed - mode interlaminar fracture toughness, Composite Structures, Volume 160, 2017, Pages 446 - 456, ISSN 0263 - 8223. [24] A. Matzenmiller, J . Lubliner, R. Taylor. A constitutive model for anisotropic damage in fiber - composites Mech Mater, 20 (1995), pp. 125 - 152. [25] LS - [26] D. M. Liljedahl, C & D. Crocombe, A & Abdel Wahab, Magd & Ashcroft, Ian. (2006). Damage modelling of adhesively bonded joints. International Journal of Fracture. 141. 147 - 161. 10.1007/s10704 - 006 - 0072 - 9. [27] A. Turon, P.P. Camanho, J. Costa, J. Renart, Accurate simulation of delamination growth under mixed - mode l oading using cohesive elements: Definition of interlaminar strengths and elastic stiffness, Composite Structures, Volume 92, Issue 8, 2010, Pages 1857 - 1864, ISSN 0263 - 8223, https://doi.org/10.1016/j.compstruct.2010.01.012. merical Simulation of Mixed - Mode Progressive Delamination in 1415 1438. [29] Seung J. Song and Anthony M. Waas . "Energy - based mechanical model for mixed mode failure of laminated composites", AIAA Journal, Vol. 33, No. 4 (1995), pp. 739 - 745. https://doi.org/10.2514/3.12639. 61 Journal o f Composite Materials , vol. 40, no. 22, Nov. 2006, pp. 2025 2046, doi:10.1177/0021998306061320. [31] A.P. Mouritz, Review of z - pinned composite laminates, Composites Part A: Applied Science and Manufacturing, Volume 38, Issue 12, 2007, Pages 2383 - 2397, ISS N 1359 - 835X, https://doi.org/10.1016/j.compositesa.2007.08.016. [32] Cox, B & Narayanaswamy, Sridhar. (2002). A Traction Law for Inclined Fiber Tows Bridging Mixed - Mode Cracks. Mechanics of Advanced Materials and Structures. 9. 299 - 331. 10.1080/15376490290 096973. [33] Chen, Leishan & G. Ifju, Peter & Sankar, Bhavani. (2005). Analysis of Mode I and Mode II Tests for Composites with Translaminar Reinforcements. Journal of Composite Materials - J COMPOS MATER. 39. 1311 - 1333. 10.1177/0021998305050425. [34] ASTM D5528 Standard test method for mode I interlaminar fracture toughness of unidirectional fiber - reinforced polymer matrix composites. In: Annual book of ASTM standards vol. 15.06. USA: ASTM;2002. [35] ASTM, "Standard Test Method for Determination of the Mode II Interlaminar Fracture Toughness of Unidirectional Fiber - Reinforced Polymer Matrix Composites," in ASTM D7905/D7905M - 14, ed. West Conshohocken, PA: ASTM International 2014. [36] - DYNA Support Site, www.dynasupport.com/tutorial/contact - modeling - in - ls - dyna/contact - types. [37] LS - - DYNA R7.1, May 26, 2014. Epox [39] Gareth W. Beckermann, Kim L. Pickering, Mode I and Mode II interlaminar fracture toughness of composite laminates interleaved with electrospun nanofibre veils, Composites Part A: Applied Science and Manufacturing, Volume 72, 2015, Page s 11 - 21, ISSN 1359 - 835X. 62 [40] I.S. Floros, K.I. Tserpes, T. Löbel, Mode - I, mode - II and mixed - mode I+II fracture behavior of composite bonded joints: Experimental characterization and numerical simulation, Composites Part B: Engineering, Volume 78, 2015, Pages 459 - 468, ISSN 1359 - 8368. [41] Yasuyo Tanzawa, Naoyuki Watanabe, Takashi Ishikawa, Interlaminar fracture toughness of 3 - D orthogonal interlocked fabric composites, Composites Science and Technology, Volume 59, Issue 8, 1999, Pages 1261 - 1270, htt ps://doi.org/10.1016/S0266 - 3538(98)00167 - 5. [42] Kimberley A. Dransfield, Laut K. Jain, Yiu - Wing Mai, On the effects of stitching in CFRPs I. mode I delamination toughness, Composites Science and Technology, Volume 58, Issue 6, 1998, Pages 815 - 827, ISSN 02 66 - 3538, https://doi.org/10.1016/S0266 - 3538(97)00229 - 7. [43] A. Turon, P.P. Camanho, J. Costa, C.G. Dávila, A damage model for the simulation of delamination in advanced composites under variable - mode loading, Mechanics of Materials, Volume 38, Issue 11, 2 006, Pages 1072 - 1089, ISSN 0167 - 6636, https://doi.org/10.1016/j.mechmat.2005.10.003.