INTEGRATION OF PROCESS SIMULATION WITH DAMAGE MODELLING OF COMPOSITE LAMINATES USING LAYERWISE ELEMENTS by Hamidreza Bakhtiarizadeh B.Sc., Sharif University of Technology, 2011 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Civil Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) August 2015 © Hamidreza Bakhtiarizadeh, 2015 ii Abstract During the last decades, several simulation models have been introduced in numerical modelling of composite structures. Research works in this area are focused on either manufacturing process or in-service behaviour analysis of composite materials. Although, these two phenomena are not separable and the effects of curing process on in-service behaviour of composite material is inevitable in practice, process-induced residual stresses are rarely considered in the in-service analysis of materials. The main objective of this research is to combine the process simulation and in-service analysis of composite materials within a common finite element framework. This would make it possible to model a complete life cycle of composite material from the beginning of the curing process all the way to its failure. In this approach the process-induced residual stresses are carried over to the in-service analysis and failure of the material. The model is implemented in an existing framework for process modelling developed at UBC Composites Group. Minimizing computational cost and accuracy are two important objectives in modelling of composite structures and usually there is a trade-off between these two objectives. In this research with the objective of minimizing computational cost and having the capability of capturing through-thickness stresses accurately, the Layerwise element is selected for the finite element modelling framework. The performance of the integrated process/damage simulation framework is tested through modelling of flat curved composite laminates that undergo various processing (cure) cycles and are subsequently subjected to mechanical loads. iii Table of Contents Abstract .......................................................................................................................................... ii Table of Contents ......................................................................................................................... iii List of Tables ............................................................................................................................... vii List of Figures ............................................................................................................................. viii List of Symbols .............................................................................................................................xv Acknowledgements .................................................................................................................. xviii Chapter 1: Introduction ................................................................................................................1 Chapter 2: Literature Review .......................................................................................................4 2.1 Composite Laminated Plate Theories ............................................................................. 4 2.1.1 Classical Laminated Plate Theory .............................................................................. 4 2.1.2 First-Order Laminated Plate Theory ........................................................................... 5 2.1.3 Third-Order Plate Theory ........................................................................................... 6 2.1.4 Layerwise Theory ....................................................................................................... 8 2.1.5 Displacement Field Assumption ................................................................................. 9 2.1.6 Strain-Displacement Relation ................................................................................... 12 2.2 Finite Element Formulation of Layerwise Theory........................................................ 13 2.2.1 Displacement Domain Discretization ....................................................................... 13 2.2.2 Strain-Displacement Relation ................................................................................... 14 2.3 Composite Material Process Modelling ........................................................................ 15 2.3.1 Process Modelling Approaches................................................................................. 15 2.3.2 COMPRO Component Architecture (CCA) ............................................................. 18 2.3.3 Thermal-Stress Analysis ........................................................................................... 19 iv 2.4 Figures........................................................................................................................... 20 Chapter 3: Methodology..............................................................................................................24 3.1 Objectives ..................................................................................................................... 24 3.2 Process Modelling Formulation .................................................................................... 24 3.2.1 Element Type ............................................................................................................ 25 3.2.2 Thermochemical Module .......................................................................................... 25 3.2.3 Stress-Deformation Module ...................................................................................... 26 3.3 In-Plane Progressive Damage Modelling ..................................................................... 28 3.4 Linking Process Simulation with In-Plane Damage Modelling ................................... 33 3.4.1 COMPRO and CODAM Inconsistencies.................................................................. 33 3.4.2 Combined Modelling Approach ............................................................................... 34 3.4.3 Algorithm .................................................................................................................. 36 3.4.3.1 Pseudo-Code ..................................................................................................... 37 3.6 Tables ............................................................................................................................ 39 3.8 Figures........................................................................................................................... 40 Chapter 4: Verification................................................................................................................41 4.1 Modelling of the Elastic Plate ....................................................................................... 41 4.1.1 Isotropic Homogeneous Material .............................................................................. 42 4.1.1.1 Results ............................................................................................................... 43 4.1.2 Composite Plate ........................................................................................................ 45 4.1.2.1 Results ............................................................................................................... 46 4.1.3 Laminated Plate under Tension: Free-edge Effect .................................................... 47 4.1.3.1 Results ............................................................................................................... 47 v 4.2 Process Modelling ......................................................................................................... 48 4.2.1 Flat Geometry ........................................................................................................... 48 4.2.1.1 Problem Description ......................................................................................... 49 4.2.1.2 Results ............................................................................................................... 50 4.2.1.3 Deformation ...................................................................................................... 50 4.2.1.4 Transverse Normal Stress Comparison for Different Cure Cycles ................... 51 4.2.1.5 Free Edge Stresses Before and After Tool-Removal ........................................ 51 4.2.1.6 Residual stress development during the curing process ................................... 53 4.2.2 Curved Geometry ...................................................................................................... 55 4.2.2.1 Results ............................................................................................................... 55 4.3 Progressive Damage Modelling .................................................................................... 57 4.3.1 Unidirectional laminate ............................................................................................. 59 4.3.2 Quasi-isotropic Laminate .......................................................................................... 62 4.3.2.1 Saturation Strain Calculation ............................................................................ 62 4.3.2.2 Simulation Results ............................................................................................ 63 4.4 Tables ............................................................................................................................ 65 4.6 Figures........................................................................................................................... 66 Chapter 5: Numerical Examples ................................................................................................93 5.1 Effect of Residual Stresses on Damage Initiation of Laminated Plates ........................ 93 5.1.1 Problem Description ................................................................................................. 93 5.1.2 Results ....................................................................................................................... 94 5.2 Curved-Part Delamination ............................................................................................ 94 5.2.1 Problem Description ................................................................................................. 95 vi 5.2.2 Results ....................................................................................................................... 96 5.3 Figures........................................................................................................................... 98 Chapter 6: Conclusion ...............................................................................................................104 Bibliography ...............................................................................................................................105 vii List of Tables Table 3.1 Inputs and outputs of CCA module .............................................................................. 39 Table 4.1 Mechanical properties of Steel ..................................................................................... 65 Table 4.2 Mechanical properties of fully cured AS4/8552 FRP ................................................... 65 viii List of Figures Figure 2.1 Stresses equilibrium on the layers interface of a composite plate. Transverse stresses are continuous across the interface of adjacent layers .................................................................. 20 Figure 2.2 Through-thickness shape functions of layerwise element for (a) linear and (b) quadratic interplation through the thickness. The shape functions are shown for an element with 3 numerical layers through its thickness. ...................................................................................... 21 Figure 2.3 Nodes notation for a layerwise element with 2 numerical layers, a quadratic shape function through thickness and a linear shape function in-plane .................................................. 22 Figure 2.4 Flow chart showing the integrated sub-model approach for process modelling [Arafath 2007] ............................................................................................................................... 22 Figure 2.5 Schematic showing the sub-model structure used in COMPRO [Johnston (1997)] ... 23 Figure 3.1 The flowchart of the combined approach analysis ...................................................... 40 Figure 4.1 Schematic of the flat homogenous, isotropic plate with material properties of steel under end shear load ..................................................................................................................... 66 Figure 4.2 Deflection of the bottom plane of the plate shown in Figure 4.1 calculated using the LT element with two different meshes compared with the analysis results using the ABAQUS 3D continuum elements result ...................................................................................................... 66 Figure 4.3 Through-thickness distribution ofσx obtained at the centre of the plate (x=25 mm, y=5mm) using the LT element with two different meshes compared with the analysis results using the ABAQUS 3D continuum elements result ..................................................................... 67 Figure 4.4 Through-thickness distribution of the transverse shear stress τxz obtained at all in-plane integration points of the LT elements with (a) full and (b) reduced integration scheme. ix Also, in-plane integration points numbering for a typical LT element is shown. The average values compared with the results of ABAQUS 3D continuum elements modelling for the plate-bending problem shown in Figure 4.1. The average values agree with the 3D modelling of the problem. ........................................................................................................................................ 68 Figure 4.5Transverseshearstressτxz distribution through the thickness of the plate shown in Figure 4.1 calculated at the centre point using both coarse and fine mesh of LT elements with reduced integration scheme compared to the corresponding ABAQUS 3D element modelling result. ............................................................................................................................................. 69 Figure 4.6 Schematic of a laminated plate with [0/90/0] layup under an end shear load ............. 69 Figure 4.7 Deflection of the bottom surface of the laminated plate shown in Figure 4.6 using LT elements and compared with ABAQUS 3D continuum elements ................................................ 70 Figure 4.8 Through-thicknessdistributionoftheaxialstressσx measured at the centre of the plate (x=25mm, y=5mm) .............................................................................................................. 70 Figure 4.9 Transverse shear stress distribution through the thickness of the laminated plate calculated at the centre point of the plate using reduced integrated LT element and compared with ABAQUS 3D element modelling results .............................................................................. 71 Figure 4.10 Transverse shear strain distribution through the thickness of the laminated plate calculated at the centre point of the plate using reduced integrated LT element and compared with ABAQUS 3D element modelling results .............................................................................. 71 Figure 4.11 Schematic of a laminated plate with (a) [0/90/0] and (b)[90/0/90] layup under axial load.Thefreeedgepeelstressσz is calculated near the free edge using LT element and ABAQUS 3D continuum element. ............................................................................................... 72 x Figure 4.12 Transverse normal inter-laminar stress distribution at an interface between 0° and 90° layers near the free edge calculated using the LT element and compared with ABAQUS 3D continuum element modelling for (a)[0/90/0] and (b)[90/0/90] layups ........................................ 73 Figure 4.13 Schematic of a cross-ply composite plate being cured on a rigid invar tool ............. 74 Figure 4.14 Two different cure cycles used for curing of the laminate shown in Figure 4.13 ..... 74 Figure 4.15 Deformation profile after tool-removal for the laminate shown in Figure 4.13 undergoing Cure Cycle number 1 shown in Figure 4.14 .............................................................. 75 Figure 4.16 Transverse normal stress calculated at the interface between layer 1 (0°) and layer 2 (90°) after tool-removal and subsequent uniaxial loading for two different cure cycles. Also the results for the case without cure modelling are shown. ................................................................ 75 Figure 4.17 Average stress on layer 1 (bottom 0° ply) and layer 2 (90° ply) distribution along the path near the free edge before and after tool-removal, compared with COMPRO 3D modelling 76 Figure 4.18Averageσy stress evaluated the bottom layer after tool-removal calculated using COMPRO 3D element and LT element modelling for the plate shown in Figure 4.13 under Cure Cycle 1 .......................................................................................................................................... 76 Figure 4.19 Inter-laminartransverseshearstressτyz at the interface between layers 1 (0°) and 2 (90°) y-axis (with the origin at the free edge) before and after tool-removal, compared with COMPRO 3D results .................................................................................................................... 77 Figure 4.20Transversenormalstress(σz ) at the interface between layers 1 (0°) and 2 (90°) y-axis (with the origin at the free edge) before and after tool-removal, compared with COMPRO 3D results ...................................................................................................................................... 78 xi Figure 4.21 Transverse shear stress distribution through the thickness of laminate on the free edge and at a point far from the free edge before and after tool-removal.σz is increased on the free edge and also it is reduced after tool-removal at the bottom of the plate. ............................. 79 Figure 4.22 Residual stresses evolution during the cure cycle, compared with COMPRO 3D modelling ...................................................................................................................................... 80 Figure 4.23 Schematic of the curved part with layup of [0/90]s made from AS4/8552 CFRP cured on a rigid tool. ............................................................................................................................... 80 Figure 4.24 The one hold cure cycle is used for the curved part curing and the resin degree of cure during the processing simulation in COMPRO and Layerwise element .............................. 81 Figure 4.25 Curved part deformation after tool-removal in r-θplane.Tomakethedeformationvisible in the figure, displacements are scaled by a factor of 10 .................................................. 81 Figure 4.26 (a)Radial and (b)tangential displacement of the curved part at r = 50mm calculated from the LT element modelling compared with COMPRO 3D modelling for a unidirectional part ([0] lay-up) .................................................................................................................................... 82 Figure 4.27 (a)Radial and (b)tangential displacement of the curved part at r = 50mm calculated from layerwise element modelling compared with COMPRO 3D modelling for a cross-ply part ([0/90]s lay-up) .............................................................................................................................. 83 Figure 4.28 In-plane stress distribution through the thickness of the curved [0°,90°]s part at θ=45°beforeandaftertool-removal calculated using layerwise element modelling compared with COMPRO 3D modelling ...................................................................................................... 84 Figure 4.29 Transverse normal (or radial) stress distribution through the thickness of the curved [0°,90°]s partatθ=45°atdifferentcuringstagescalculatedfromlayerwiseelementmodellingcompared with COMPRO 3D modelling ...................................................................................... 84 xii Figure 4.30 Inter-laminar transverse normal stress variation at the interface between the top 0º and90ºpliesatθ=45°sectionofthecurvedpartduringthecurecycleobtainedfrommodellingwith layerwise element and compared with COMPRO 3D modelling ......................................... 85 Figure 4.31 Schematic figure of the OCT test modelled using layerwise element with damage modelling capability in ABAQUS implicit................................................................................... 86 Figure 4.32 A rectangular area (51mm×30mm) around the notch tip has a uniform fine mesh with element size of 1mm×1mm .................................................................................................. 86 Figure 4.33 The FE model of an OCT test on a 90° unidirectional lamina made of AS4/8552 CFRP. The modelling with higher order layerwise element failed to converge. The red circle shows the area with damaged elements. ....................................................................................... 87 Figure 4.34 Simulation of OCT where the damage zone is limited to one row of elements ........ 87 Figure 4.35 OCT test modelled using CODAM in ABAQUS implicit wit 2D (a) linear (4-noded) and (b) quadratic (8-noded) elements. The linear element modelling successfully converged to the end of complete specimen failure however the modelling with the higher order elements failed to converge after damage occurred in the first element ...................................................... 88 Figure 4.36 The middle node of damaged elements show an unexpected behaviour. The highest residual force (which makes the results diverge) occurs on the middle node of the element which is about to get damaged because the stiffness at both integration points above and below the middle node goes to zero. ............................................................................................................. 89 Figure 4.37 The symmetry of the problem is applied to the model by fixing the middle node of the crack-band elements. The solution is converges successfully until complete failure of the unidirectional (90º) lamina............................................................................................................ 89 xiii Figure 4.38 The load applied load versus CMOD in the OCT test simulation of a unidirectional (90º) composite plate in ABAQUS implicit using the higher order layerwise element. .............. 90 Figure 4.39 Stress vs. strain from unidirectional strain simulation of one-element to find the fracture energy for different values of fibre and matrix saturation strains. .................................. 90 Figure 4.40 Energy density vs. damage saturation strain plotted to find the fibre and matrix damagesaturationstrainswhichcorrespondstotheenergydensityofγs=110 MPa ................... 91 Figure 4.41 The OCT test simulation of [90/45/0/-45]s laminate made from Hexcel AS4/8552 CFRP using layerwise element in ABAQUS implicit comparing with simulation of the same laminate using an explicit solver in LS-DYNA. ........................................................................... 91 Figure 4.42 The OCT test simulation of [90/45/0/-45]s laminate made from AS4/8552 using layerwise element in ABAQUS implicit with to different damage saturation strain parameter for CODAM. It shows the behaviour of the model with higher saturation strain is smoother and consequently its chance of convergence is higher than the one with lower damage saturation strain. ............................................................................................................................................. 92 Figure 5.1 Schematic of a [0/90]s laminate made from AS4/8552 CFRP : (a) cured on a rigid invar tool in autoclave using two different cure cycles and (b) the laminate is subsequently removed from the tool and subjected to an external uniform in-plane displacement. .................. 98 Figure 5.2 Two different cure cycles used for curing of the laminate in Figure 5.1 .................... 99 Figure 5.3 Laminate stiffness vs. applied displacement for two different cure cycles and without manufacturing simulation. Matrix damage initiation occurs with lower applied displacement for the ones with cure modelling because of the residual stresses in the matrix. ............................... 99 Figure 5.4 Schematic figure of the curved-part modelled for detecting delamination ............... 100 xiv Figure 5.5 A displacement applied to both ends of the curved-part until delamination occurs in one of the plies interfaces ........................................................................................................... 100 Figure 5.6 Transverse normal stress distribution in the curved part when a displacement of 1mm is applied to both ends without process modelling. The stress has not reached the critical value fordelaminationσcr=79.4MPa (which corresponds 0.8% strain in matrix). The highest stress occurs atInterface1andatsectionθ=45º. .................................................................................. 101 Figure 5.7 Transverse normal stress on 3 potential delamination interfaces when a displacement of 1mm is applied to both ends of the curved part without process modelling. ......................... 101 Figure 5.8 The temperature cycle applied to the cuved part for process modelling ................... 102 Figure 5.9 The curved-part is modelled with two different (a) convex and (b) concave rigid tools...................................................................................................................................................... 102 Figure 5.10TransversenormalstressatInterface1andatsectionθ=45ºasafunctionoftheapplied displacement. While the figure illustrates the modelling with curing simulation shows lower value of the residual stress in the material, there is no significant difference between the results for convex and concave tools. ......................................................................................... 103 xv List of Symbols 𝑨 Laminate’sin-plane stiffness matrix in 2D formulation 𝑩 Strain-displacement matrix of the element C 3D stiffness matrix 𝐸1 LongitudinalYoung’smodulus of a composite lamina 𝐸2 Transverse Young’smodulus of a composite lamina 𝐸3 TransverseYoung’s modulus of a composite lamina 𝑭 External load vector 𝐺12 In-plane shear modulus of a composite lamina G13 In-plane shear modulus of a composite lamina 𝐺23 Transverse shear modulus of a composite lamina ℎ Plate thickness ℎ𝑒 Element height 𝑲 Stiffness matrix 𝑲𝑒𝑙𝑒 Stiffness matrix of the element 𝑵𝑱 In-plane shape function of the element for 𝐽𝑡ℎ node 𝑸 The plane stress in-plane stiffness matrix of a lamina 𝑅 Stiffness reduction coefficient 𝑹 Internal load vector 𝑹𝑒𝑙𝑒 Internal load vector 𝑻 Transformation matrix for strain 𝚫𝒖 Incremental displacement vector xvi 𝒖 Displacement vector 𝑢, 𝑣, 𝑤 Displacement in x, y and z 𝑈, 𝑉,𝑊 Nodal displacement in x, y and z 𝑢0, 𝑣0, 𝑤0 Mid-plane displacement of the plate in x, y and z 𝑤𝑖 Weight factors for the numerical integration 𝛼 Resin degree of cure 𝛾𝑠 Energy density of the element 𝝐𝑡𝑜𝑡𝑎𝑙 Total strain vector 𝝐𝑚𝑒𝑐ℎ Mechanical strain vector 𝝐𝑓𝑟𝑒𝑒 Free strain vector ?̂? Effective mechanical strain 𝜖𝑖 Damage initiation strain 𝜖𝑠 Damage saturation strain 𝝈 Stress vector Φ Trough-thickness shape function of the layerwise element for in-plane displacement Ψ Trough-thickness shape function of the layerwise element for transverse displacement 𝜔 Damage parameter Subscript 𝑚 Matrix 𝑓 Fibre xvii Superscript 𝑖 Initiation 𝑠 Saturation 𝑑 Damaged 0 Undamaged xviii Acknowledgements ByGod’swill,aidandsupportthecompletionofthisthesishasbecomeareality. I would like to acknowledge the help of those without whom I could not have finished this work. Firstly, my deepest gratitude goes to my supervisor, Dr. Reza Vaziri, for his guidance, encouragement and patience throughout this work. His valuable suggestions and constructive criticism have made this thesis a better product. I would also thank Dr. Anoush Poursartip for his support during the course of my studies. I would like to thank Dr. Alireza Forghani for his encouragement and providing constructive feedbacks during all stages of my studies. I would like to thank former and present members of Composite Research Network, Dr. Ahmed Arafath, Kamyar Gordnian, Dr. Navid Zobeiry, Dr. Sardar Malek Mohammadi, Mina Shahbazi and Sina Amini for their helpful discussions and joyful company throughout this work. I am also indebted to Payam Saffari for answering my numerical questions. Finally, my special thanks are due to my mother, father and sister for all their supports. Most of all, I would like to thank my lovely wife, Ghazal, for her patience and sacrifice throughout my studies as her words of kindness and support were my constant source of motivation. 1 Chapter 1: Introduction Fibre-Reinforced Polymer (FRP) materials are replacing traditional materials in advanced structural applications due to their high strength and stiffness, low weight and durability. They are being used in a several structures such as aircrafts, cars, storage tanks, pressure vessels and sport equipment. They also have applications in civil engineering such as retrofit of existing structures. The high-performance structural components made of advanced thermoset matrix composites, which is of interest in this research, are commonly manufactured using an autoclave process. This process starts from stacking pre-impregnated sheets of unidirectional fibres (prepreg) at different orientations followed by subjecting the whole part to a controlled cycle of temperature and pressure inside an autoclave. During the curing process, due to the cure shrinkage of matrix, mismatch between mechanical properties of fibres and matrix and geometric features, internal stresses are developed within the material, which leads to deformations in the structure after tool-removal. Although sometimes these stresses are relieved due to the deformation, they usually do not vanish completely. These residual stresses along with other outcomes of processing such as voids and micro-cracks affect the in-service behaviour of the materials and its strength. Therefore an accurate simulation of composite structures in service and predicting their failure requires an accurate incorporation of manufacturing-induced effects. 2 During the past decades, several numerical models have been developed to model the behaviour of composite materials both during processing and in-service. However, there have been very few attempts at integrating these two simulations. In other words, process-induced effects have scarcely been considered in damage modelling of composite materials. The goal of this research work is to develop a numerical modelling framework for composite materials that is capable of modelling a composite part from the beginning of the manufacturing to the end of its complete failure seamlessly within a common finite element mesh. To do so, the framework should consist of the following attributes: Efficiency: The numerical cost of this model should be minimized by using efficient numerical techniques and element type, and simplified models. Accuracy: The model should be capable of capturing both transverse (out-of-plane) and in-plane stresses accurately. Completeness of modes of failure: Damage modelling within this framework should comprehensive and include the ability to detect the initiation and propagation of the different modes of failure in composites such as delamination, matrix or fibre failure. Based on these objectives the thesis is organized as follows: Chapter 2: Literature review. In this chapter, a brief literature review of composite laminated plate theories suitable for 3D analysis is presented and the appropriate theory is selected for this research. The finite element formulation for this selected theory is then described. Finally a brief review of literature on process modelling is presented in the last section of this chapter. 3 Chapter 3: Methodology. In this chapter, the selected process modelling and damage modelling formulation is presented and then the incompatibilities of available process and damage modellings are investigated. Finally, a methodology for modification and linkage of these two modelling regimes is proposed. Chapter 4: Verification. To verify the capability of the proposed modelling framework, in this chapter, a few simple examples of composite material during the curing process and under different loading conditions are modelled using this framework and the proposed element. Then the result of analysis is compared with 3D analysis of the same examples using 3D continuum elements in ABAQUS. Chapter 5: Numerical case studies. In this chapter, two different composite structures are modelled through the different curing conditions and the effect of considering process-induced residual stresses on damage initiation point is investigated. In the first case study a flat composite pare is modelled using two different cure cycles to understand the effect of cure cycles in damage initiation in composite materials, and in the second case study the cure process of a curved composite part is modelled using two different convex and concave rigid tools. 4 Chapter 2: Literature Review 2.1 Composite Laminated Plate Theories Modelling of composite laminated plates is a difficult process due to the complex interaction of fiber and matrix, and layers bonding. It could be modelled through a 3D simulation of the plate or simplifying the plate to a 2D model. The high aspect ratio of plates makes their 3D modelling more difficult than other types of structures. In order to maintain a reasonable aspect ratio, the size of the elements should be small which leads to a high computational cost. Generally, composite laminates are modelled as 2D plates with equivalent properties calculated from Laminated Plates Theories. However, during the smearing process to simplify the problem to a 2D model, some of the important aspects that are needed to describe the through thickness behaviour will be lost. The objective of this section is to give an overview of plate theories with the aim of introducing a method suitable for effective modelling of through thickness effects. 2.1.1 Classical Laminated Plate Theory The classical laminated plate theory (CLPT) is an extension of the classical plate theory to composite laminates (Reissner 1961; Whitney1969; Reddy 1997). In this theory it is assumed that the Kirchhoff hypothesis holds: 1. Planes perpendicular to the mid-surface before deformation remain plane after deformation. 2. The transverse normals do not experience elongation. 3. The transverse normals remain perpendicular to the mid-surface after deformation. Kirchhoff hypothesis in plate bending, leads to the following equations for the displacement field in CLPT 5 { 𝑢(𝑥, 𝑦, 𝑧, 𝑡) = 𝑢0(𝑥, 𝑦, 𝑡) − 𝑧𝜕𝑤0𝜕𝑥𝑣(𝑥, 𝑦, 𝑧, 𝑡) = 𝑣0(𝑥, 𝑦, 𝑡) − 𝑧𝜕𝑤0𝜕𝑦𝑤(𝑥, 𝑦, 𝑧, 𝑡) = 𝑤0(𝑥, 𝑦, 𝑡) (2.1) where 𝑢, 𝑣 and 𝑤 are displacement field in each of the 𝑥, 𝑦 and 𝑧 andthesuffix‘0’returntothemid-plane values. Based on these equations, the transverse displacement is independent of the normal coordinate, and therefore, the derivative of the transverse displacement with respect to 𝑧 is zero, which implies that the transverse normal strain is also zero i.e. 𝜖𝑧𝑧 = 0. The third assumption results in zero transverse shear strains (𝜖𝑥𝑧 = 0, 𝜖𝑦𝑧 = 0). Therefore CLPT is not capable of predicting out-of-plane strains and hence stresses through the thickness of composite laminated plates. 2.1.2 First-Order Laminated Plate Theory The hypotheses in the first-order shear deformation laminated plate theory (FSDT) are the same as CLPT; except that the third assumption is relaxed (Whitney 1970; Reddy 1997). Accordingly, the plane sections remain plane after deformation but they do not remain perpendicular to the mid-surface. Therefore the displacement field could be written as {𝑢(𝑥, 𝑦, 𝑧, 𝑡) = 𝑢0(𝑥, 𝑦, 𝑡) − 𝑧𝜙𝑥(𝑥, 𝑦, 𝑡)𝑣(𝑥, 𝑦, 𝑧, 𝑡) = 𝑣0(𝑥, 𝑦, 𝑡) − 𝑧𝜙𝑦(𝑥, 𝑦, 𝑡)𝑤(𝑥, 𝑦, 𝑧, 𝑡) = 𝑤0(𝑥, 𝑦, 𝑡) (2.2) The transverse shear strain in this theory is given by 6 𝜖𝑥𝑧 =12(𝜕𝑢𝜕𝑧+𝜕𝑤𝜕𝑥)=12(−𝜙(𝑥, 𝑦, 𝑡) +𝜕𝑤0(𝑥, 𝑦, 𝑡)𝜕𝑥) (2.3) Equation (2.3) shows that although the transverse shear strains in FSDT is nonzero, it is not a function of the transverse coordinate, 𝑧. In other words 𝜖𝑥𝑧 is constant through the thickness. Therefore, FSDT is not capable of capturing the variation of transverse stresses through the thickness. 2.1.3 Third-Order Plate Theory While CLPT and FSDT theories adequately represent the plate’sbehavior,higher-order theories are capable of detecting transverse (out-of-plane) stresses more accurately (Reddy 1984; Reddy 1990). The idea in such theories is to expand the assumed in-plane displacement field of the plate to a certain power of the transverse coordinate, z. For example, the third-order plate theory displacement components are defined as {𝑢(𝑥, 𝑦, 𝑧, 𝑡) = 𝑢0 + 𝑧𝜙𝑥 + 𝑧2𝜃𝑥 + 𝑧3𝜆𝑥𝑣(𝑥, 𝑦, 𝑧, 𝑡) = 𝑣0 + 𝑧𝜙𝑦 + 𝑧2𝜃𝑦 + 𝑧3𝜆𝑦𝑤(𝑥, 𝑦, 𝑧, 𝑡) = 𝑤0 (2.4) in which the mid-plane displacement and notations are given by: 𝑢0 = 𝑢(𝑥, 𝑦, 0, 𝑡) 𝑣0 = 𝑣(𝑥, 𝑦, 0, 𝑡) 𝑤0 = 𝑤(𝑥, 𝑦, 0, 𝑡)𝜙𝑥 = (𝜕𝑢𝜕𝑧)𝑧=0𝜙𝑦 = (𝜕𝑣𝜕𝑧)𝑧=02𝜃𝑥 = (𝜕2𝑢𝜕𝑧2)𝑧=02𝜃𝑦 = (𝜕2𝑣𝜕𝑧2)𝑧=06𝜆𝑥 = (𝜕3𝑢𝜕𝑧3)𝑧=06𝜆𝑦 = (𝜕3𝑣𝜕𝑧3)𝑧=0 (2.5) 7 There are 9 dependent unknowns which can be reduced by imposing traction-free boundary conditions 𝜎𝑥𝑧(𝑥, 𝑦, ±ℎ/2, 𝑡) = 0 and 𝜎𝑦𝑧(𝑥, 𝑦, ±ℎ/2, 𝑡) = 0, where ℎ is the plate thickness. After applying boundary conditions on the displacement field and calculating strains, the transverse shear stress can be written as follows: {𝛾𝑥𝑧𝛾𝑦𝑧} = {𝛾𝑥𝑧(0)𝛾𝑦𝑧(0)} + 𝑧2 {𝛾𝑥𝑧(2)𝛾𝑦𝑧(2)} (2.6) Equation (2.6) shows that the transverse shear stress is a quadratic function with respect to the transverse coordinate, z. This calculation is more accurate than FSDT and CLPT, and it is usually sufficient for a single layer of plate as it can predict the exact transverse stress when it is a quadratic function through the thickness. However, when it comes to composite laminated plates, the quadratic distribution does not accurately represent the stress distribution estimation. For equilibrium of inter-laminar force the transverse stresses should be continuous through the thickness of the laminate while the in-plane stresses can be discontinuous. The following stress continuity and discontinuity conditions hold between the layers’ interfaces through the thickness. {𝜎𝑥𝑥𝜎𝑦𝑦𝜎𝑥𝑦}𝑡𝑘≠ {𝜎𝑥𝑥𝜎𝑦𝑦𝜎𝑥𝑦}𝑏𝑘+1 , {𝜎𝑥𝑧𝜎𝑦𝑧𝜎𝑧𝑧}𝑡𝑘= {𝜎𝑥𝑧𝜎𝑦𝑧𝜎𝑧𝑧}𝑏𝑘+1 (2.7) where 𝑡𝑘 refers to the top surface of the 𝑘𝑡ℎ layer and 𝑏𝑘+1 refers to the bottom surface of the (𝑘 + 1)𝑡ℎ layer, and (𝑘 + 1)𝑡ℎ is on top of the 𝑘𝑡ℎ layer as Figure 2.1 shows. Since different layers have different stiffness properties, the stress continuity on the layers’ boundaries, will cause discontinuity in strains. 8 {𝛾𝑥𝑧𝛾𝑦𝑧𝛾𝑧𝑧}𝑡𝑘≠ {𝛾𝑥𝑧𝛾𝑦𝑧𝛾𝑧𝑧}𝑏𝑘+1 (2.8) Note that the transverse shear strain obtained from the third-order theory is continuous at the layers’boundaries due to the continuous displacement field assumption through the thickness. However, the equilibrium equations result in discontinuous strain field through the thickness. It can be concluded that in order to have the capability of detecting accurate transverse stresses in composite laminates using the Equivalent Single Layer (ESL) theory, the first derivative of displacement field should be discontinuous across layer interfaces leading to discontinuous strains through the thickness of plate. The Layerwise Theory, introduced by Reddy (1990), and presented in the next section provided such capability for laminated plate analysis. 2.1.4 Layerwise Theory As shown in the previous section, ESL theories do not predict the transverse stresses accurately. Modelling of a plate using 3D FEM model is a way to overcome this modelling problem. Because a 3D model is capable of capturing transverse stresses accurately when its mesh is fine enough. A 3D model with one element through the thickness clearly cannot predict accurate stressed through the thickness as it suffers from the same deficiencies as the ESL theories, because its displacement field is continuous through the thickness thus requiring a mesh refinement is needed through the thickness. However, 3D finite element modelling with many elements through the thickness becomes computationally costly as the requirement for keeping the aspect ratio of elements reasonable leads to large number of elements in the model. Layerwise Theory (LT), developed by Reddy (1990), is a powerful approach for laminated plates 9 modelling. The layerwise element is capable of predicting transverse stresses accurately with one element through the thickness and its computational cost is less than a robust 3D model. Moreover, to capture through-thickness response, unlike using conventional 3D element, there is no need to re-mesh the whole structure. Instead, results can be refined by increasing the number of layers within the layerwise element. The strain continuity problem is addressed in layerwise theory through changing the displacement field assumptions of ESL theories. The displacement field assumption in LT exhibits only 𝐶0 continuity on the laminate layers boundaries through the thickness. In other words, while the transverse displacement is continuous, its derivative can be discontinuous through the plate thickness. This feature provides the ability to the LT to model discontinuous transverse strains (𝜖𝑥𝑧 , 𝜖𝑦𝑧 , 𝜖𝑧𝑧) at the boundaries of the layers. 2.1.5 Displacement Field Assumption In LT, the displacements in the plate are written as { 𝑢(𝑥, 𝑦, 𝑧, 𝑡) =∑𝑈𝐼(𝑥, 𝑦, 𝑡)Φ𝐼(𝑧)𝑚𝐼=1𝑣(𝑥, 𝑦, 𝑧, 𝑡) =∑𝑉𝐼(𝑥, 𝑦, 𝑡)Φ𝐼(𝑧)𝑚𝐼=1𝑤(𝑥, 𝑦, 𝑧, 𝑡) =∑𝑊𝐼(𝑥, 𝑦, 𝑡)Ψ𝐼(𝑧)𝑛𝐼=1 (2.9) where 𝑢, 𝑣 and 𝑤 denote the total displacement components in the x, y and z-directions, Φ𝐼(𝑧) and Ψ𝐼(𝑧) are continuous shape functions of the through thickness coordinate, 𝑚 is the number of nodes through the thickness of laminate for displacement in x and y coordinates, and 𝑛 is the 10 number of nodes through the thickness for transverse displacement interpolation. (𝑈𝐼 , 𝑉𝐼 ,𝑊𝐼) represents the displacements at the 𝐼𝑡ℎ plane at 𝑧 = 𝑧𝐼. In the LT, there are numerical layers through the thickness, which are used for displacement interpolation through the thickness. These layers are not necessarily same as the physical layers in the laminated composite materials. Each numerical layer can have more or less than one material layer inside. However, typically each numerical layer represents a material layer of the laminate. The function {Φ𝐼(𝑧)} is a shape function of the in-plane displacements’ interpolation through the thickness. This function is independent of x and y, and used for z-direction discretization. As shown in Figure 2.2 show, ΦI is a continuous function through the thickness of the plate, however its first derivative can be discontinuousonnumericallayers’boundaries. It has a value of 1 on 𝐼𝑡ℎ node through the thickness and 0 on all other nodes. The order of this interpolation function withineachnumericallayerdependsonthenumberofthelayer’snodes.For instance, as Figure 2.2–(a) shows, 𝛷 consists of linear functions within each layer when number of nodes on each layer is 2, whereas when there are 3 nodesoneachlayer(twoonthelayer’sboundaryand one in the middle), the function consists of quadratic functions on each numerical layer. The following formulation represents the 𝛷 function when there are 2 nodes through the thickness within each layer, hence the function would be a piecewise linear function (Figure 2.2): 11 { Φ1(𝑧) = 𝜓11(𝑧), 𝑧1 ≤ 𝑧 ≤ 𝑧2Φ𝐼(𝑧) = {𝜓2𝐼−1(𝑧),𝜓1𝐼(𝑧),𝑧𝐼−1 ≤ 𝑧 ≤ 𝑧𝐼+1𝑧𝐼 ≤ 𝑧 ≤ 𝑧𝐼+2Φ1(𝑧) = 𝜓2𝑁(𝑧), 𝑧𝑁−1 ≤ 𝑧 ≤ 𝑧𝑁 {𝜓1𝑘 =𝑧 − 𝑧𝑘+1𝑧𝑘 − 𝑧𝑘+1𝜓2𝑘 =𝑧 − 𝑧𝑘𝑧𝑘+1 − 𝑧𝑘 (2.10) where z is the through thickness coordinate, 𝑧𝑘 denotes the z-coordinate of bottom surface of the 𝑘𝑡ℎ numerical layer (which is equivalent to the location of the 𝑘𝑡ℎ node), and N is the number of numerical layers through the plate thickness. Figure 2.2 shows profiles of through thickness shape functions (Φ) in 3 layers element for linear (2 nodes within each numerical layer) and quadratic interpolation (3 nodes within each numerical layer) through the element thickness. Due to the fact that the in-plane displacement fields in the z-direction are constructed from piecewise interpolation functions, the displacement derivatives are continuous within each layer but are discontinuous across the layer interfaces. Ψ𝐽(𝑧) functions are shape functions in the z-direction for transverse displacement discretization. They can be the same as, or different from, Φ functions (in-plane displacement shape functions). The order of this functions depends on the number of nodes through the thickness for transverse displacement interpolation. For example, when the material is assumed to be incompressible the 12 transverse displacement (𝑤) will be assumed to be constant through the thickness. In such case Ψ(𝑧) = 1. When the number of nodes for both the transverse and in-plane displacements are same, an identical interpolation functions are used for Φ and Ψ. In this research as the ultimate goal is to model delamination in composite laminated plates, it is needed to predict transverse displacement variation through the thickness, so that more than one node through the thickness is needed for transverse displacement interpolation. Therefore for simplification of the numerical analysis formulation the same nodes and shape functions are used for interpolation of transverse(𝑤) and in-plane (𝑢 and 𝑣) displacements, hence the through thickness interpolation functions for in-plane and transverse displacement is identical (𝛹 = 𝛷). 2.1.6 Strain-Displacement Relation Strains can be written as follows, using displacement fields in Equation (2.9): 𝜖𝑥𝑥 =∑𝜕𝑈𝐼𝜕𝑥Φ𝐼𝑁𝐼=1𝜖𝑦𝑦 =∑𝜕𝑉𝐼𝜕𝑦Φ𝐼𝑁𝐼=1𝜖𝑧𝑧 =∑𝑊𝐼𝜕Φ𝐼𝜕𝑧𝑁𝐼=1𝜖𝑥𝑦 =∑𝜕𝑈𝐼𝜕𝑦Φ𝐼𝑁𝐼=1+∑𝜕𝑉𝐼𝜕𝑥Φ𝐼𝑁𝐼=1𝜖𝑥𝑧 =∑𝑈𝐼𝜕Φ𝐼𝜕𝑧𝑁𝐼=1+∑𝜕𝑊𝐼𝜕𝑥Φ𝐼𝑁𝐼=1𝜖𝑥𝑦 =∑𝑉𝐼𝜕Φ𝐼𝜕𝑧𝑁𝐼=1+∑𝜕𝑊𝐼𝜕𝑦Φ𝐼𝑁𝐼=1 (2.11) 13 In the above equations the number of nodes is assumed to be identical for both the in-plane and transverse displacements. Equation (2.11) shows that the strain components 𝜖𝑧𝑧, 𝜖𝑥𝑧 and 𝜖𝑦𝑧 contain a term that has derivatives of Φ, and as discussed, this term is discontinuous through the thickness across layer boundaries. Therefore the Layerwise displacement field assumption can have discontinuous transverse shear stress on boundaries. 2.2 Finite Element Formulation of Layerwise Theory 2.2.1 Displacement Domain Discretization Part of the displacement field discretization for finite element modelling is done in Equation (2.9), where the displacement field is discretized in the z-direction. To complete the discretization in terms of x and y coordinates using 2D continuum elements. The function (𝑈𝐼 , 𝑉𝐼,𝑊𝐼) in Equation (2.9) should be written terms of 2D shape functions in x-y domain. Equation (2.12) represents the discretized displacement field in the x-y plane. { 𝑈𝐼(𝑥, 𝑦, 𝑡) =∑𝑈𝐼𝐽(𝑡)𝑁𝐽(𝑥, 𝑦)𝑝𝐽=1𝑉𝐼(𝑥, 𝑦, 𝑡) =∑𝑉𝐼𝐽(𝑡)𝑁𝐽(𝑥, 𝑦)𝑝𝐽=1𝑊𝐼(𝑥, 𝑦, 𝑡) =∑𝑊𝐼𝐽(𝑡)𝑁𝐽(𝑥, 𝑦)𝑝𝐽=1 (2.12) in which 𝑁𝐽(𝑥, 𝑦) are 2D Lagrangian interpolation functions and 𝑝 is number of nodes in the 2D continuum element. For example, the element can be a 4-noded (𝑝 = 4) rectangular or 8-noded (𝑝 = 4) quadratic (Serendipity) element. 14 The displacement field is discretized in two steps based on LT assumptions. In the first step (Equation (2.9)) it was discretized in the z-direction by different shape functions for the in-plane and transverse displacements, and in the second step, discretization was done in the in-plane coordinates using the same shape functions 𝑁(𝑥, 𝑦) in all directions. By combining Equation (2.9) and Equation (2.12) the fully discretized form of the displacement field in LT can be obtained as follows: { 𝑢(𝑥, 𝑦, 𝑧, 𝑡) =∑∑𝑁(𝑥, 𝑦)Φ𝐼(𝑧)𝑈𝐼𝐽(𝑡)𝑝𝐽=1𝑚𝐼=1𝑣(𝑥, 𝑦, 𝑧, 𝑡) =∑∑𝑁(𝑥, 𝑦)Φ𝐼(𝑧)𝑉𝐼𝐽(𝑡)𝑝𝐽=1𝑚𝐼=1𝑤(𝑥, 𝑦, 𝑧, 𝑡) =∑∑𝑁(𝑥, 𝑦)Φ𝐼(𝑧)𝑊𝐼𝐽(𝑡)𝑝𝐽=1𝑚𝐼=1 (2.13) in which 𝑁 and Φ are the in-plane and through thickness shape functions, respectively. 2.2.2 Strain-Displacement Relation By assuming small deformation in the plate, the displacement gradients will be small in all directions. Therefore, the square (and higher powers), and products of displacement gradients are negligible in the Green-Lagrange strain, and as a result, linear strain-displacement equations will be used in finite element modelling. By substituting the discretized displacement fields (Equation (2.13)) in the linear strain-displacement relation (Equation (2.11)), the strain components can be written in the following discretized form: 15 𝜖𝑥𝑥 =∑∑𝜕𝑁𝐽𝜕𝑥Φ𝐼𝑈𝐼𝐽𝑝𝐽=1𝑁𝐼=1𝜖𝑦𝑦 =∑∑𝜕𝑁𝐽𝜕𝑦Φ𝐼𝑉𝐼𝐽𝑝𝐽=1𝑁𝐼=1𝜖𝑧𝑧 =∑∑𝑁𝐽𝜕Φ𝐼𝜕𝑧𝑊𝐼𝐽𝑝𝐽=1𝑁𝐼=1𝜖𝑥𝑦 =∑∑𝜕𝑁𝐽𝜕𝑦Φ𝐼𝑈𝐼𝐽𝑝𝐽=1𝑁𝐼=1+∑∑𝜕𝑁𝐽𝜕𝑥Φ𝐼𝑉𝐼𝐽𝑝𝐽=1𝑁𝐼=1𝜖𝑥𝑧 =∑∑𝑁𝐽𝜕Φ𝐼𝜕𝑧𝑈𝐼𝐽𝑝𝐽=1𝑁𝐼=1+∑∑𝜕𝑁𝐽𝜕𝑥Φ𝐼𝑊𝐼𝐽𝑝𝐽=1𝑁𝐼=1𝜖𝑦𝑧 =∑∑𝑁𝐽𝜕Φ𝐼𝜕𝑧𝑉𝐼𝐽𝑝𝐽=1𝑁𝐼=1+∑∑𝜕𝑁𝐽𝜕𝑦Φ𝐼𝑊𝐼𝐽𝑝𝐽=1𝑁𝐼=1 (2.14) in which (𝑈𝐼𝐽, 𝑉𝐼𝐽,𝑊𝐼𝐽) is the displacement vector annotated with the node number (𝐼, 𝐽). Node (𝐼, 𝐽) corresponds to the 𝐼𝑡ℎ node through the thickness and 𝐽𝑡ℎ node of the 2D element. Figure 2.3 shows the nodes notation for a typical layerwise element with 2 numerical layers, a quadratic shape function through the thickness and linear shape function in-plane. 2.3 Composite Material Process Modelling 2.3.1 Process Modelling Approaches Many complex physical phenomena such as heat transfer, resin flow, and stress development are present in manufacturing process of composite materials. Achieving an accurate simulation of composite manufacturing process depends on accuracy in modelling these phenomena. “Integratedsub-model”approach,firstintroducedbyLoosandSpringer(1983), is a widespread approach in process modelling of composite materials. This method divides the process modelling into three simpler sub-models, which could be studied more or less independently 16 from each other. Figure 2.4 shows this approach, which has been used by many researchers such as Bogetti and Gillespie (1991,1992), White and Hahn (1992), Johnston et al. (2001), and Zhu et al (2003). The governing equation, the material behavior and the boundary conditions are the three different but related components of each sub-model. The method used to solve the governing equations, the material constitutive behavior approximation, and representation of the physical boundary conditions, all directly influence the accuracy of the performance of aforementioned sub-models. One of the approximations in solving the governing equations is the analysis dimension. The dimension of the governing equations solution is chosen to be 1D, 2D or 3D [Loos and Springer (1983), Bogetti and Gillespie, (1991, 1992), Johnston et al. (2001), Zhu et al, (2003)]. The other approximation is the numerical technique used in solving the governing equations. Different methods have been used in the literature for this purpose, such as laminate plate theory [Hahn and Pagano (1975), Loos and Springer (1983), White and Hahn (1992) and Bogetti and Gillespie (1992)], finite difference method [Bogetti and Gillespie (1991)] and the finite element technique [Chen and Ramkumar (1988), Johnston et al. (2001), and Zhu et al (2003)]. The state of resin in Fibre Reinforced Polymers (FRP) transitions from viscous fluid to the elastic solid during the manufacturing process in an autoclave. The resin behaviour can be viscous fluid, viscoelastic solid or elastic solid in different stages of material curing process. The accuracy of constitutive models for simulation of the resin mechanical behaviour during the cure 17 cycle, determine the accuracy of the stress development model. Many constitutive models have been used in the literature for calculation of residual stresses developed during the manufacturing such as purely elastic [Loos and Springer (1983), Nelson and Cairns (1989)], Cure Hardening Instantaneously Linear Elastic or CHILE [Bogetti and Gillespie (1992), Johnston et al. (2001)], pseudo-viscoelastic [Zobeiry et al. (2006a)] and fully viscoelastic [Kim and White (1996), Zhu et al (2003), Zobeiry et al. (2006b)]. Among all processing simulation constitutive models, CHILE models are popular in the literature because it is easy to characterize and implement. CHILE models are essentially elastic at each time increment, however the elasticity modulus changes as a function of temperature and degree of cure. In other words: 𝛥𝜎 = 𝐸(𝑇, 𝛼)𝛥𝜖 (2.15) in which 𝛥𝜎 is the incremental mechanical stress, 𝐸(𝑇, 𝛼) is the instantaneously elasticity modulus as a function of temperature (𝑇) and degree of cure (𝛼), and 𝛥𝜖 is the incremental strain. Also, it can be written in the integral form as a follow: 𝜎(𝑡) = ∫𝐸(𝑇, 𝛼)𝑑𝜖𝑑𝜏𝑑𝜏𝑡0 (2.16) in which, 𝜎 is the stress, t is the current time, and 𝜏 is the time integration variable. According to the nature of the material behaviour, it is expected that only fully viscoelastic formulations for polymers predict accurate deformations and stresses, however the CHILE constitutive models have been shown to provide good predictions in analysis of the thermoset polymers behaviour (Fernlund et al. 2003). 18 As mentioned above, various numerical techniques combined with different material constitutive models in a specific dimension can be used for process modelling of composite structures. The process models available in the literature mostly consist of a certain numerical technique and a specific material constitutive model. The integrated sub-model approach is the one used in development of COMPRO at UBC, in which the communication among these sub-models is achieved through three different state variables, namely, temperature, degree of cure of the resin, and the fibre volume fraction. The finite element method (FEM) is the numerical technique used in this approach, followed by the CHILE material constitutive model. COMPRO has the capability of implementing the material behavior in a very modular fashion. As discussed before, one of the most common techniques in numerical modelling of physical problems is FEM. It is also widely used in modelling of composite material processing. The choice of the types of elements are used for modelling depends on the specifications of the problem. An efficient and verified finite element process model should have the capability of using different types of elements and material models. 2.3.2 COMPRO Component Architecture (CCA) As Figure 2.5 shows, modelling in CCA is performed using different material subroutines, namely, cure kinetics, resin properties, fibre properties, thermal expansion, cure shrinkage, and micromechanics modules. Different material models for each subroutine have been implemented in CCA for different composite materials. In order to model a new material using CCA the appropriate material models need to be addedtothe“Modulebank”.Theinterfacemoduleisalink between the material model modules and commercial finite element codes. Due to the fact 19 that separate modules are used for material models that interact with the FE codes, CCA can be linked to any other commercial FE code by modifying only the interface module while the rest of modules remain unchanged. 2.3.3 Thermal-Stress Analysis In this study, the thermochemical and the stress modules are the only modules that are exercised. A uniform fibre volume fraction without any resin flow is assumed. When the stress/displacement solution depends on the temperature field and there is no inverse dependency, a sequentially coupled thermal-stress analysis can be used which is performed through two steps; the first step is to solve the pure heat transfer problem to obtain the temperature solution. This solution is then fed into stress analysis as a predefined temperature field. Temperature, as a pre-defined variable read into the FE code (ABAQUS) at the nodes, and then interpolated to the desired points within elements as needed. According to the brief literature review of composite laminated plate analysis, presented in this chapter, layerwise element is selected as the FEM element for modelling of laminated plates in this study. In addition, according to the presented literature on composite material process modelling, CHILE approach, which has already been implemented in the CCA, is selected as the process modelling framework. In the next chapter the selected damage modelling framework and its linkage to process modelling will be introduced. 20 2.4 Figures Figure 2.1 Stresses equilibrium on the layers interface of a composite plate. Transverse stresses are continuous across the interface of adjacent layers 12…kk+1…n21 (a) (b) Figure 2.2 Through-thickness shape functions of layerwise element for (a) linear and (b) quadratic interplation through the thickness. The shape functions are shown for an element with 3 numerical layers through its thickness. 010 1Through Thickness Coordinate (z/h)Φ(z)010 1Through Thickness Coordinate (z/h)Φ(z)010 1Through Thickness Coordinate (z/h)Φ(z)010 1Through Thickness Coordinate (z/h)Φ(z)01-0.5 0 0.5 1Through Thickness Coordinate (z/h)Φ1(z)01-0.5 0 0.5 1Through Thickness Coordinate (z/h)Φ2(z)01-0.5 0 0.5 1Through Thickness Coordinate (z/h)Φ3(z)01-0.5 0 0.5 1Through Thickness Coordinate (z/h)Φ4(z)01-0.5 0 0.5 1Through Thickness Coordinate (z/h)Φ5(z)01-0.5 0 0.5 1Through Thickness Coordinate (z/h)Φ6(z)01-0.5 0 0.5 1Through Thickness Coordinate (z/h)Φ7(z)22 Figure 2.3 Nodes notation for a layerwise element with 2 numerical layers, a quadratic shape function through thickness and a linear shape function in-plane Figure 2.4 Flow chart showing the integrated sub-model approach for process modelling (Arafath 2007) (1,1)(2,1)(3,1)(4,1)(5,1)(5,4)(1,3)(2,3)(3,3)(4,3)(5,3)(1,2)(2,2)(3,2)(4,2)(5,2)Governing EquationMaterial BehaviourBoundary ConditionTHERMOCHEMICAL MODELResin Degree of CureTemperatureViscosityFLOW MODELResin PressureFibre Volume FractionSTRESS MODELStressesStrainsDisplacementsGoverning EquationMaterial BehaviourBoundary ConditionGoverning EquationMaterial BehaviourBoundary ConditionPROCESS VARABLESTemperaturePressure23 Figure 2.5 Schematic showing the sub-model structure used in COMPRO (Johnston 1997) State variablesFlowModuleStressModuleThermochemicalModuleOutputModuleInputModuleAutoclaveSimulationModuleMaterial Module24 Chapter 3: Methodology 3.1 Objectives Typically two different finite element modelling frameworks are available for manufacturing simulation and in-service structural analysis of composite structures. The objective of this research is to expand the ability of composite material modelling by combining these two frameworks into one integrated framework. To do so, a progressive in-plane damage modelling approach (CODAM) and COMPRO Component Architecture (CCA) for process modelling was implemented within a user element in the commercial FE code, ABAQUS Standard (or implicit). To improve the capability of this element in predicting both in-plane and out-of-plane damages (delamination) layerwise element is chosen, which has the capability of capturing both in-plane and out-of-plane (transverse) stresses accurately. 3.2 Process Modelling Formulation A modular approach is used here for process modelling of composite materials. CCA approach has been implemented as a user subroutine linked with ABAQUS finite element code (Arafath, 2007). Although, as explained in Section 2.3, different material models are available in CCA for numerous composite materials, only the materials model of HEXCEL AS4/8552 Carbon Fibre Reinforced Polymer (CFRP) has been implemented within the ABAQUS User Element (UEL). However, the capability for modelling different materials can be added to the UEL by simply changing the material model subroutines since the element formulation subroutine is different from the material model subroutines. The material model modules of HEXCEL AS4/8552 CFRP are availableinArafath’sDoctoralDissertation,AppendixA(Arafath, 2007). 25 3.2.1 Element Type Based on the discussion in Section 2.1 only the 3D elements are capable of predicting the through-thickness stresses accurately albeit in a computationally costly manner. Therefore a more efficient solid-like shell element with such capability, e.g. an element based on Layerwise Theory (LT) has been implemented as a User Element (UEL) in ABAQUS implicit by (Arafath, 2007). This type of element is not available as a built-in internal element of ABAQUS since it has multiple nodes through the thickness and its node numbers vary based on the number of numerical layers within the element. 3.2.2 Thermochemical Module The temperature and degree of cure of the resin are two major state variables which are used in calculation of other material properties. These two parameters are handled with thermochemical module. The thermochemical module consists of two different analysis for heat transfer and resin kinetics analysis. Table 3.1 shows the input and outputs of the thermochemical module. As explained in Section 2.3.3, the heat transfer analysis is done using a sequentially coupled method; whereas it is assumed that the stress/displacement solution does not affect the heat transfer analysis. Therefore in this study no heat transfer analysis is performed, instead the cure cycle temperature is applied to the whole material uniformly as a pre-defined variable. Using the temperature variable passed to the thermochemical module, the cure kinetics analysis is performed at all integration points of the element and it calculates the degree of cure rate of corresponding point as a function of material properties and temperature. In this study, the 26 AS4/8552 SFRP cure kinetics model is used, which was developed by Johnston and Hubert (1995). 3.2.3 Stress-Deformation Module The stress module is responsible for predicting strains and stresses in the material based on the given displacement. This module is based on the cure-hardening instantaneously linear elastic (CHILE) model introduced by (Johnston, 1997). In this model the curing process is modelled within smaller time increments. In each increment the resin properties are upgraded and then using the calculated material properties, the constitutive equations are solved with the assumption of linear elastic behaviour. The material response is assumed to be elastic in each time-increment, while the resin properties evolve during the cure process. In the implemented stress-deformation module in ABAQUS, there are two major subroutines: UEXPAN and UMAT. These subroutines are called at each integration point. UEXPAN is responsible to calculate incremental free strains due to CTE and cure shrinkage of the composite. The UMAT subroutine is used to define material (fibre and resin) properties at each time increment. Table 3.1 shows the input and outputs of the stress-deformation and micro-mechanics module. The fibre and resin properties are passed to the micromechanics module to calculate stiffness of the material at each integration point. Then the stiffness matrix (𝑪) is used to calculate the incremental stress at the integration point. 27 𝛥𝝈 = 𝑪 𝛥𝝐𝑚𝑒𝑐ℎ = 𝑪 (𝛥𝝐𝑡𝑜𝑡𝑎𝑙 − 𝛥𝝐𝑓𝑟𝑒𝑒) (3.1) where 𝛥𝝐𝑚𝑒𝑐ℎ is the incremental mechanical strain, 𝑪 is the instantaneous stiffness matrix and 𝛥𝝈 is the incremental stress at a given integration point. In the time step analysis, the total stress at the current time 𝑡𝑛, is calculated by adding the incremental stress calculated in the current increment to the total stress calculated in the previous step. 𝝈(𝑡+Δ𝑡) = 𝝈(𝑡) + 𝛥𝝈 (3.2) Subsequently the nodal forces of the element is calculated: 𝑹𝑒𝑙𝑒 = − ∫ 𝑩𝑇𝝈𝑑𝑉𝑒𝑙𝑒𝑽𝑒𝑙𝑒 (3.3) Where 𝑹𝑒𝑙𝑒 is the residual force vector of the element. The material stiffness matrix (𝑫) is also used for calculation of the element stiffness matrix: 𝑲𝑒𝑙𝑒 = ∫ 𝑩𝑇𝑪𝑩𝑑𝑉𝑒𝑙𝑒𝑽𝑒𝑙𝑒 (3.4) The governing finite element equation for the stress module, which is passed to the ABAQUS solver, is as follows: ∑𝑲𝑒𝑙𝑒𝒖𝒏𝒆𝒍𝒆= ∑𝑹𝑒𝑙𝑒𝒏𝒆𝒍𝒆 (3.5) Generally,ABAQUSsolvestheoverallsystemofequationsusingtheNewton’smethod: Solve: 𝑲𝜟𝒖 = 𝜳 Set: 𝒖 = 𝒖 + 𝜟𝒖 Iterate (3.6) where 𝜳 is the residual load vector. 28 3.3 In-Plane Progressive Damage Modelling The first generation of UBC Composite Damage Model (CODAM) was introduced by Williams et al. (Williams et al. 2003) for modelling damage in composite laminated plates at the sub-laminate level. The objective of composite damage models at the sub-laminate level is to predict the effective and overall behaviour of large-scale laminates. A simple model that represents the overall behaviour of the material in terms of macro-scale parameters, such as crack patterns and absorbed energy, can predict the response of actual scale models fairly well. McClennan (McClennan 2004) showed that a simple scalar damage model is able to predict the behaviour of notched specimens accurately, without predicting the micro-scale properties of the material such as formation of microcracks, matrix cracking and fibre breakage, while all of them contribute in the overall behaviour of the material. Since the overall load-displacement response of the material is mainly controlled by the laminate peak stress and the absorbed energy. In this model the laminates stiffness is degraded due to the effect of matrix cracking and fibre breakage.Thereduced laminate’sstiffness iswrittenasa functionof threedamageparametersthat represent the reduction of elasticity moduli in two in-plane principal directions and the shear modulus in the same coordinate system. The original CODAM formulation shows a dependency on the choice of principal coordinate system. Due to the fact that damage parameters were written in the in-plane principal directions of the laminate, these to two perpendicular directions become the proffered direction in damage 29 growth in this model. To resolve the objectivity problem of the original version, Forghani (2011, 2013) extended CODAM to include an orthotropic non-local formulation. The new formulation called CODAM2 was developed based on the components of strain in the lamina principal directions. CODAM2 reduces the material stiffness due to matrix cracking and fibre breakage in each layer and then the equivalent in-plane stiffness matrix of the laminate is captured based on integrating the layer’s stiffnesses. The model was implemented in a 2D plane stress formulation in the open-source FE code (OOFEM) as well as the commercial FE code, LS-DYNA. Only one integration point through the thickness of the element was used to evaluate the equivalent stiffness matrix of the laminate. The in-plane stiffness matrix of the laminate can be calculated according to the following equation: 𝑨𝑑 =∑𝑻𝑘𝑇 𝑸𝑘𝑑 𝑻𝑘𝑡𝑘𝑛𝑘=1 (3.7) where 𝑨𝑑 is the damaged in-plane stiffness matrix of the laminate, 𝑻𝑘 is the transformation matrix for stresses and strains, and 𝑡𝑘 is the thickness of the 𝑘𝑡ℎ ply of the laminate. According to the CODAM formulation 𝑸𝑘𝑑, the damaged stiffness matrix of 𝑘𝑡ℎ ply, is written as 30 𝑸𝑘𝑑 =[ 𝑅𝑓𝐸11 − 𝑅𝑓𝑅𝑚𝜈12𝜈21𝑅𝑓𝑅𝑚𝐸21 − 𝑅𝑓𝑅𝑚𝜈12𝜈210𝑅𝑚𝐸21 − 𝑅𝑓𝑅𝑚𝜈12𝜈210𝑆𝑌𝑀 𝑅𝑚𝐺12] (3.8) where 𝑅𝑓 and 𝑅𝑚 are respectively the fibre and matrix damage parameters. 𝐸1, 𝐸2 and 𝐺12 represent the longitudinal, transverse and shear moduli and, 𝜈12 and 𝜈21 are the major and minor Poisson’sratios. The fibre and matrix stiffness reduction factors (𝑅𝑓 and 𝑅𝑚) vary linearly with fibre and matrix damage parameters ranging from 1 to 0. 𝑅𝑓 = 1 − 𝜔𝑓𝑅𝑚 = 1 − 𝜔𝑚 (3.9) Where ωf and ωm are respectively the fibre and matrix damage parameters. In CODAM damage parameters are calculated based on the equivalent strains functions. The equivalent strain function that is used to update matrix damage parameter is a function of the transverse (𝜖22) and shear (𝛾12) components of the local strain in the composite layer. The equivalent strain function that affects fibre damage parameter is calculated based on the longitudinal normal strain of the composite layer (𝜖11). Following equations show the formulation of equivalent strain functions. 𝜖𝑚 = 𝑠𝑖𝑔𝑛(𝜖22)√(𝜖22)2 + (𝛾122)2 𝜖𝑓 = 𝜖11 (3.10) Then the damage parameters are defined according to the following equation 31 𝜔𝑚 =(𝜖𝑚−𝜖𝑚𝑖 )(𝜖𝑚𝑠 −𝜖𝑚𝑖 )𝜖𝑚𝑠𝜖𝑚 𝜔𝑓 =(𝜖𝑓−𝜖𝑓𝑖 )(𝜖𝑓𝑠−𝜖𝑓𝑖 )𝜖𝑓𝑠𝜖𝑓 for (𝜖𝑚 − 𝜖𝑚𝑖 ) > 0 (3.11) where 𝜖𝑓𝑖 and 𝜖𝑚𝑖 are damage initiation strains of fibre and matrix, and 𝜖𝑓𝑠 and 𝜖𝑚𝑠 are damage saturation strains of fibre and matrix, respectively. Although the original formulation of CODAM was designed for 2D plane stress element, the model is upgraded in this study to become capable of being implemented 3D layered element. In the 2D formulation of CODAM, in-plane elasticity moduli (𝐸1 and 𝐸2), shear modulus (𝐺12) and Poisson’s ratios (𝜈12 and 𝜈21) are updated using fibre and matrix damage parameters (first four equations of Equation (3.13). In the 3D formulation 𝜈13 and 𝜈23 are assumed to be reduced by the fibre and matrix reduction factors because of fibre breakage and matrix cracks, respectively. It is also assumed that there is no matrix cracking in the out-of-plane direction of the composite laminate, therefore 𝐸3, 𝐺13, 𝐺23, 𝜈31 and 𝜈32 will remain intact during the in-service behaviour of the material. Damaged transverse elastic properties (𝜈13𝑑 , 𝜈23𝑑 ) updated as such to satisfy reciprocity. The undamaged material elastic properties also satisfy the reciprocity. So that the following equation should be satisfied in the damaged material: 𝐸1𝑑𝜈13𝑑 =𝐸3𝜈31=𝐸1𝜈13 𝐸2𝑑𝜈23𝑑 =𝐸3𝜈32=𝐸2𝜈23 (3.12) 32 Therefore 𝐸1 and 𝜈13 damage parameters should be identical. Similarly, the damage parameters of 𝐸2 and 𝜈23 should be identical. Updated material properties in the 3D formulation of the CODAM is defined as follows: 𝐸1𝑑 = 𝑅𝑓𝐸1𝐸2𝑑 = 𝑅𝑚𝐸2𝜈12𝑑 = 𝑅𝑓𝜈12𝐺12𝑑 = 𝑅𝑚𝐺12𝜈21𝑑 = 𝑅𝑚𝜈21𝜈13𝑑 = 𝑅𝑓𝜈13𝜈23𝑑 = 𝑅𝑚𝜈23 (3.13) Although, the upgraded CODAM model, presented in this study has not been validated with experiments, its plane-stress condition is verified for a single element with the original 2D formulation of CODAM implemented in LS-DYNA. In damage modelling within the 3D layerwise element, each integration point belongs to one layer of the material. Therefore, the material stiffness matrix at each integration point in the global coordinate system is obtained from the following equation: 𝑪𝐺𝑑 = 𝑻𝑇 𝑪𝐿𝑑 𝑻 (3.14) where 𝑪𝑮𝒅 and 𝑪𝑳𝒅 are the damaged stiffness matrix in the global and local coordinate system, respectively. 𝑻 is the standard transformation (rotation) matrix for the strain vector. The calculation of 3D stiffness matrix in the local coordinate system is presented in the following equation: 33 𝑪𝐿𝑑 =[ (1 − 𝜈23𝑑 𝜈32)𝐴∗. 𝐸2𝑑𝐸3(𝜈21𝑑 + 𝜈23𝑑 𝜈31)𝐴∗. 𝐸2𝑑𝐸3(𝜈31 + 𝜈21𝑑 𝜈32)𝐴∗. 𝐸2𝑑𝐸30 0 0(1 − 𝜈13𝑑 𝜈31)𝐴∗. 𝐸1𝑑𝐸3(𝜈32 + 𝜈12𝑑 𝜈31)𝐴∗. 𝐸1𝑑𝐸30 0 0(1 − 𝜈12𝑑 𝜈21𝑑 )𝐴∗. 𝐸1𝑑𝐸2𝑑 0 0 0𝐺12𝑑 0 0𝐺13 0𝑆𝑌𝑀 𝐺23] with 𝐴∗ defined as: 𝐴∗ =(1−𝜈12𝑑 𝜈21𝑑 −𝜈23𝑑 𝜈32−𝜈31𝜈13𝑑 −2𝜈21𝑑 𝜈32𝜈13𝑑 )𝐸1𝑑𝐸2𝑑𝐸3 (3.15) In the above equation, 𝐸𝑑, 𝐺𝑑 and 𝜈𝑑 respectively represent the damaged elasticity modulus, shear modulus andPoisson’sratio. As Equation (3.15) shows, in this research CODAM formulation was extended to a 3D form to become compatible with COMPRO 3D process modelling and implemented within the LT element for modelling both the processing and the in-plane progressive damage of composite material. 3.4 Linking Process Simulation with In-Plane Damage Modelling 3.4.1 COMPRO and CODAM Inconsistencies In the modular approach simulation, multiple modules run in sequence and each module is responsible for modelling a specific phenomenon. They receive state variables as input from the central database and pass the output to it. As all these modules work with common state variables, inconsistencies between their assumptions can cause problems in modelling. 34 As outlined in Section 3.2.3, at each time-step in the CHILE approach, the incremental mechanical strains arising from the applied external load and thermochemical strains generated in the material during the process are used to calculate the stress increments. However, CODAM uses the total strain as the basis for calculating the secant stiffness of the damaged material. The summation of all incremental mechanical strains during the curing process does not have a physical meaning given that each increment of mechanical strain is associated with a different state of the material as curing advances and therefore the material stiffness evolves. Moreover, in the CHILE approach the total strain, which can be calculated based on total displacement, is not an appropriate measure for predicting damage in the material, since part of the material deformation that occurs before resin vitrification (i.e. the deformation during the viscous fluid state of the resin) does not lead to significant stress (and hence damage) in the material. Physically, the quantity that better represents the state of the material for the purposes of damage modelling is the accumulated stress, i.e. the stresses that are computed using Equation (3.2). 3.4.2 Combined Modelling Approach In the combined approach, first, the CCA framework is used to determine the material stiffness and stresses of the undamaged material as they evolve during processing. CODAM is subsequently used to quantify damage, if any, and perform the corresponding stiffness reductions and stress computations. 35 To address the incompatibility issue of the process and damage modelling approaches, mentioned above, an effective mechanical strain parameter (?̂?) is introduced to be used instead of the total strain in CODAM. This parameter is basically used to translate the mechanical stresses into their equivalent strain, because it is assumed that in the curing material the high level of stress will represent the existence of damage, not high level of total strain. On the other hand in CODAM is a strain-based damage model. Therefore the effective mechanical strain, which is used to updated the damage parameters in the combined approach, is defined as ?̂? = (𝑪0)−1 𝝈0 (3.16) where 𝝈0 is the current stress calculated from the processing model and 𝑪0 is the current stiffness of the undamaged material. It should be noted that ?̂? is a fictitious quantity and does not represent the real state of strain in the material. It is merely the strain equivalence of the current state of stress in the material. In other words, ?̂? is the mechanical strain that would have resulted in a state of stress 𝝈0, in a virtual, linear elastic material with constant stiffness 𝑪0. The damage model is then formulated on the basis of the above-defined effective mechanical strain thus retaining the advantages of continuum damage mechanics formulations in the strain space. The proposed model formulation under a specific condition, in which the material is elastic and the elastic properties of the material are constant, is identical to the CODAM model. In such case the effective mechanical strain (?̂?) is equal to the total mechanical strain. To show that CODAM is a special case of the proposed formulation, the effective mechanical strain is rewritten based on the assumption of constant material stiffness. ?̂? as a function of incremental undamaged mechanical stress is: 36 ?̂? = (𝑪0)−1𝝈0 = (𝑪0)−1∑Δ𝝈𝑖𝟎 (3.17) Δ𝛔i0 is the incremental undamaged stress vector in the 𝑖𝑡ℎ time increment. In the case that the material stiffness is constant during the whole analysis the equation can be rewritten as: ?̂? = ∑(𝑪0)−1𝛥𝝈𝑖𝟎 =∑Δ𝝐𝑖𝑚𝑒𝑐ℎ = 𝝐𝑚𝑒𝑐ℎ (3.18) where 𝝐𝑖𝑚𝑒𝑐ℎ is the incremental mechanical strain vector in the 𝑖𝑡ℎ time increment. 3.4.3 Algorithm At the beginning of each time increment, the nodal temperatures are passed to the thermochemical module from the central database. Then the thermochemical module calculated the temperature, Degree of Cure (DOC) and free stress vector (𝛥𝝐𝑓𝑟𝑒𝑒) on the integration points and passes them to the stress/deformation module. The stress/deformation module obtains the incremental displacements (𝛥𝒖) of an element from the central database. Then it calculates the incremental strains using the standard strain-displacement relationship as follows: 𝛥𝝐 = 𝑩 𝛥𝒖 (3.19) where 𝛥𝝐 includes the contributions from both mechanical and free strains (thermal expansion and resin cure shrinkage). To compute the incremental mechanical strains, the free thermal and cure shrinkage strains, 𝛥𝝐𝑓𝑟𝑒𝑒, obtained from the thermochemical module, need to be subtracted: 𝜟𝝐𝑚𝑒𝑐ℎ = 𝜟𝝐 − 𝜟𝝐𝑓𝑟𝑒𝑒 (3.20) Then the incremental effective stress is calculated: 37 𝜟𝝈𝟎 = 𝑪𝟎 𝜟𝝐𝑚𝑒𝑐ℎ (3.21) According to the CHILE model, the material is assumed to be linear elastic in each increment of time and 𝑪𝟎 is the undamaged material stiffness at the current increment. The total effective stress 𝝈0 is then calculated by accumulating the incremental stresses up to the current time step and used to compute the effective mechanical strain vector (?̂?) from Equation (3.16). The total effective stress, the undamaged material stiffness and the effective mechanical strain are passed to the damage module. In the damage module, the effective mechanical strain, ?̂? is then used to update the fibre and matrix damage parameters based on the CODAM formulation. The damage parameters are then used to calculate the appropriate stiffness reduction factors, which in turn act on the current undamaged stiffness matrix, computed from the stress-deformation sub-module to obtain the damaged stiffness matrix (𝑪𝑑). The stresses are then calculated as follows: 𝝈 = 𝑪𝒅?̂? (3.22) Finally, the stress vector and damaged stiffness matrix (𝝈 and 𝑪𝒅) at all integration points of the element are used to compute the element stiffness matrix and internal nodal force vector that are then passed to the central database. The flowchart of this analysis is presented in presented in Figure 3.1. 3.4.3.1 Pseudo-Code The following describes the various computational steps involved in a given time step. 38 Loop through in-plane integration points (𝑖 = 1 to 2 and 𝑗 = 1 to 2). Obtain in-plane local coordinates of the integration point (𝜉𝑖, 𝜂𝑗). Loop through thickness integration points in the layers (𝑘 = 1 𝑡𝑜 𝑛). Obtain through-thickness local coordinate (𝜉𝑘). Call subroutine TEMPERATURE to calculate the temperature at this integration point. Call CCA subroutine to calculate the undamaged stiffness matrix, 𝑪𝟎, and incremental free strain vector 𝜟𝝐𝑓𝑟𝑒𝑒. Call subroutine B_MATRIX to calculate the B matrix. Calculate the incremental strain vector 𝜟𝝐 = 𝑩𝜟𝒖. Calculate the incremental mechanical strain 𝜟𝝐𝒎𝒆𝒄𝒉 = 𝜟𝝐 − 𝜟𝝐𝒇𝒓𝒆𝒆. Calculate the incremental effective stress 𝜟𝝈𝟎 = 𝑪𝟎 𝜟𝝐𝒎𝒆𝒄𝒉 Update the total effective stress 𝝈𝒕+𝜟𝒕𝟎 = 𝝈𝒕𝟎 + 𝜟𝝈𝟎 Calculate the total effective strain ?̂? = (𝑪0)−1 𝝈0 Call subroutine DAMAGE to update the damage parameters and stiffness matrix 𝑪𝒅 based on ?̂?. Calculate the stress 𝝈 = 𝑪𝒅?̂? Calculate the element secant stiffness matrix and the internal nodal force vector and pass to the main FE solver for the next increment. 39 3.6 Tables Table 3.1 Inputs and outputs of CCA module Input Output Thermochemical module Boundary temperature Degree of Cure (DOC), Part temperature history Stress-deformation module Part temperature history, DOC, Loading history Fibre and resin elastic properties, displacement field, free strains, mechanical strain, stress, deformations 40 3.8 Figures Figure 3.1 The flowchart of the combined approach analysis 41 Chapter 4: Verification Different features of the implemented process and damage modelling are examined in this chapter. First, the capability of the LT elements in predicting deformations and stresses of elastic materials are tested through the modelling of homogeneous and composite materials. Second, the implemented process modelling framework within the LT element is tested through the analysis of curing process of composite laminates with different geometries. Finally the implemented CODAM model is verified through a single element and the Over-height Compact Tension (OCT) tests modelling. All the results are verified through a comparison with the results of the identical problem using commercial FEM codes, namely ABAQUS, COMPRO and LS-DYNA. 4.1 Modelling of the Elastic Plate To illustrate the performance of the LT element, implemented in ABAQUS as a user element, in detecting through-thickness properties of a plate, two examples are modelled and presented in this section. The first example is used to test the LT element results in modelling of homogeneous plates. In this example a homogeneous cantilever plate with isotropic material was modelled under a vertical load at its free end. The second example is a composite laminated cantilever plate with the same dimensions, load and boundary conditions as the first example. All the simulations in this section are carried out as static problems using the User Element (UEL) implemented in ABAQUS implicit. 42 The LT element model results are compared with the results of the analysis of the same problem using the ABAQUS 20-noded built-in 3D elements (C3D20). In the latter case a fine mesh with many multiple elements through the thickness of the plate (i.e. in the z-direction) was used. To prevent shear locking problems in bending, the aspect ratio of the 3D elements were chosen to be less than 4 in all simulations with 3D ABAQUS built-in elements. 4.1.1 Isotropic Homogeneous Material This example is designed to test the LT element performance in predicting through-thickness stresses in homogeneous and isotropic plates with isotropic materials. The plate is made of steel and its in-plane dimensions are 50mm by 10mm with a thickness of 1mm. Theplate’smaterialelastic properties are presented in Table 4.1.The plate dimensions chosen to be a relatively thin plate. The plate is fully constrained at one end and is under a uniform pressure loading of 1 N/mm2 on its free edge (Figure 4.1). The problem was modelled with two types of mesh, a course and a fine mesh. The coarse mesh contains 5 (5×1) and fine mesh contains 125 (25×5) elements. The element has 8 in-plane nodes, 5 numerical layers and quadratic shape functions through the thickness of each layer. As discussed above, the problem was also modelled in ABAQUS with a very fine mesh (0.1mm × 0.1mm × 0.1mm) using its built-in higher order continuum elements (C3D20) with a total number of 500,000 elements. 43 4.1.1.1 Results Figure 4.2 shows the deflection profile of the plate along its length (x-direction) obtained using LT elements compared with the 20-noded continuum element results. The displacement is calculated at the center of the bottom plane of the plate along its length (𝑦 = 5𝑚𝑚, 𝑧 = 0𝑚𝑚). The figure shows that the predicted response obtained from the LT model even with a coarse mesh are very close to those obtained from the full 3D model. Regarding in-plane stresses, Figure 4.3 shows the longitudinal stress (𝜎𝑥) distribution through the thickness predicted by the LT model and 3D simulation. The stresses are calculated at the plate center point (x = 25mm , y = 5mm). The figure shows that the LT model has results for in-plane longitudinal stresses agree very well with the results obtained from 3D model. Two different integration schemes are used for in-plane numerical integration of LT elements; the full integration scheme with 9 in-plane integration points (3 × 3) and reduced integration scheme with 4 in-plane integration points (2 × 2). Figure 4.4 (a) shows the coarse mesh results for 𝜏𝑥𝑧 through the thickness of the element at the center point of the plate for all the integration points when using a full integration scheme. As seen from the figure, the results at individual Gauss points’ results are far from the 3D analysis results. However the average stress distribution through the thickness of the element is fairly close to the 3D model results. The following formula shows the stress averaging procedure using Gauss integration technique: 44 𝜎(𝑧) =∫ 𝜎(𝑥, 𝑦, 𝑧) 𝑑𝐴Ω∫ 𝑑𝐴Ω =∑ 𝑤𝑖. 𝜎𝑖(𝑧)𝑛𝑖=1∑ 𝑤𝑖𝑛𝑖=1 (4.1) in which 𝜎(𝑧) is the average stress, 𝜎𝑖(𝑧) is the stress at the 𝑖𝑡ℎ Gauss point and 𝑤𝑖 is the corresponding integration weight in the in- plane numerical integration, and 𝛺 is the in-plane area of the element. Figure 4.4 (b) shows the predicted 𝜏𝑥𝑧 values through the thickness at all integration points as well as their average values at the centre of the element when using a reduced integration scheme. The figure shows that the values at the Gauss points in a reduced integration scheme are closer to the exact results when compares to the full integration scheme. The average of all the integration point results is even closer to the results obtained from ABAQUS 3D element modelling than the equivalent prediction from the full integration scheme, shown in Figure 4.4 (a). Based on the results of the full and reduced integration scheme, the reduced integration scheme was selected for modelling with the goal of reducing the computational cost and avoiding locking problems. All the results are reported at the center of the in-plane element as an average of the Gauss point values. 45 Figure 4.5 shows the 𝜏𝑥𝑧 distribution obtained from both coarse and fine mesh of the LT elements as well as ABAQUS 3D elements. The figure shows that the results using the LT elements converge to the 3D modelling results by refining the mesh. From this example it can be concluded that the LT element is capable of predicting both the in-plane and out-of-plane stresses in isotropic plates using only a single element through the thickness of the plate. 4.1.2 Composite Plate To examine the capability of the LT element in capturing out-of-plane properties stresses through in composite laminates, a laminated composite plate with the same size and dimensions as the plate in Section 4.1.1 was modelled under identical boundary and loading conditions. The laminate is made of Hexcel AS4/8552 CFRP composite material with a layup of [0/90/0] (see Figure 4.6). The elastic properties of this material are presented in Table 4.2. The problem was modelled using 125 (25×5) LT elements with 3 numerical layers through the thickness. The numerical integration was performed using a reduced Gaussian integration scheme in the x-y plane (4 integration points) and the Simpson’s integration scheme with 3 integration points through the thickness of each numerical layer. The problem was also modelled in ABAQUS using a fine mesh of higher order continuum elements (20-noded element). The element size in the ABAQUS 3D model was 0.1mm × 0.1mm. 46 4.1.2.1 Results Figure 4.7 shows the predicted deflection (as vertical displacement) at the bottom of the plate along the x-axis obtained from the LT element analysis compared with the ABAQUS 20-noded continuum element modelling results. The deflection is calculated at the center of the bottom surface of the plate along the length of the plate (𝑧 = 0𝑚𝑚, 𝑦 = 5𝑚𝑚). It can be seen that the deflection calculated using the LT element agrees very well with the ABAQUS 3D model results. Figure 4.8 shows the predicted in-plane axial stress distribution through the thickness of the laminate. The stresses are evaluated at the plate center (x=25mm, y=5mm). The good agreement between the LT and 3D modelling results in Figure 4.8 confirms the fact that the LT element is capable of predicting the bending stresses accurately which was also shown by Arafath (Arafath, 2007). To examine the capability of the LT element in capturing the transverse stresses in laminated plates, the transverse shear stress and strain are calculated through the thickness of the laminate at its in-plane centroid (i.e. 𝑥 = 25𝑚𝑚, 𝑦 = 5𝑚𝑚). While the Classical Laminated Plate Theory (CLPT) is not capable of predicting the transverse stresses through the thickness of laminates, Figure 4.9 shows that the transverse stresses calculated using the LT element agrees well with the 3D element modelling results. Also, Figure 4.10 shows that unlike CLPT, the LT element is capable of capturing the discontinuous transverse shear strain distribution through the thickness. 47 4.1.3 Laminated Plate under Tension: Free-edge Effect In composite laminates, the mismatch between the behaviour of adjacent layers can introduce a complex state of stresses, which can include out-of-plane inter-laminar stresses in laminates even when they are loaded in-plane. The transverse inter-laminar stresses can reach a significant level near the free edges, making it prone to delamination in laminates. Therefore, it is important to capture these through-thickness stresses in numerical modelling of composite laminated plates. To examine the performance of LT elements in predicting free edge stresses during the in-service loading of composites, two laminated composite plates with two different layups were modelled under tensile loading. Both laminates are made of AS4/8552 CFRP with a span of 50mm, width of 10mm and a thickness of 1mm. Figure 4.11 schematically shows the plate boundary conditions and loading that consists of a uniform axial elongation of 0.25mm applied at one end (5% strain along the x-axis). In the LT element analysis, the plate is modelled with a single element through the thickness with 3 numerical layers (one numerical layer for each material layer). The problem is also modelled using ABAQUS 3D element with 9 elements (20 noded continuum element) through the thickness of the laminate. 4.1.3.1 Results The normal transverse inter-laminar stress 𝜎𝑧 is obtained at the interface of the 0º and 90º plies near the free edge in the middle of the plate (𝑥 = 25𝑚𝑚). The transverse stresses develop in a relatively small region close to the free edge. Therefore the mesh should be fine enough near the 48 boundary to be capable of detecting these stresses. In the LT element analysis of this example a uniform mesh is used. The LT elements width in the y-direction is 0.1𝑚𝑚 and their length is 1.67𝑚𝑚 (3000 elements), whereas in ABAQUS 3D element modelling the elements are 0.11𝑚𝑚 × 0.11𝑚𝑚 × 0.11𝑚𝑚 (364500 elements). Figure 4.12-(a) shows the transverse normal stress at the interface of 0º and 90º layers near the free edge of the laminate with [0/90/0] layup. The high peel stresses at the free edge are generated due to the mismatch of Poisson’s ratio in the 0º and 90º layers, which results in tensile 𝜎𝑦 stresses in the 0º plies and compressive 𝜎𝑦 stresses in the 90º plies. The requirement for these 𝜎𝑦 stresses to vanish at the free edge introduces high out-of-plane stresses at the interface between 0º and 90º layers close to the free edge. Figure 4.12-(b) shows the out of-plane stress for [90/0/90] layup. It can be seen that the transverse 𝜎𝑧 stresses that were tensile in the case of [0/90/0] layup, become compressive by changing the stacking sequence to [90/0/90]. The good agreement between the LT element and 3D element results shows that the layerwise element is capable of predicting the transverse free edge stresses accurately with only one element through the thickness. 4.2 Process Modelling 4.2.1 Flat Geometry As discussed in Section 4.1.3, capturing the transverse normal stresses 𝜎𝑧 at the free edge of laminated plates, which can potentially lead to delamination, is important in modelling composite laminates. These stresses are caused either due to the process-induced deformations or 49 in-service loading of the material. In Section 4.1.3 the free edge stress distributions were investigated for a solid elastic laminate under a tensile load. In this section we study the free edge stresses during the curing process and subsequently when the laminate is subjected to a tensile load. In other words, the effect of process-induced residual stresses is considered on the in-service behaviour of the laminate by modelling the composite from the beginning of its curing process to its subsequent in-service loading all using the same finite element mesh and modelling framework. 4.2.1.1 Problem Description A cross-ply, CFRP composite laminate with a span of 102mm, a width of 9.75mm and a thickness of 1.14mm was modelled seamlessly from the beginning of processing to the end of its tensile loading condition using LT element which was implemented as a User Element (UEL) in ABAQUS implicit. The plate is made of Hexcel AS4/8552 CFRP material with [0/90/0] layup. One element with 3 numerical layers (one layer for each ply) was used through the thickness. The modelling was carried out in one finite element mesh, through 3 stages. In the first stage the composite plate was virtually cured on a rigid frictionless tool (Figure 4.13) in an autoclave using two different cure cycles as shown in Figure 4.14. In this stage a uniform cure cycle temperature was applied to both the part and the tool. Since the laminate is thin, the temperature gradient through the thickness of laminate was neglected and a uniform temperature was assumed for the whole part and tool. Then in the second stage the part was removed from the tool and finally, in the third stage, a tensile load was applied to one of the short edges of the plate while the other edge was clamped. 50 The LT element used two different integration schemes for in-plane and through-thickness numerical integration. A reduced Gauss integration scheme, with 4 integration points (2 × 2) was used for in-plane integration, and for the through-thickness integration, the Simpson’s integration scheme was used with 3 integration points through the thickness of each layer of the element. The in-plane and out-of-plane stresses were obtained through the thickness and at the layers’ interfaces near the free edge. Also, the evolution of the stresses was investigated by calculating the stresses at different times during the curing process. All the stresses were calculated at the centre of the element at different integration layers through the thickness by averaging over the four in-plane integration points. To verify of the LT element results, the same problem was modelled using COMPRO and the stresses were compared with the LT element results. The modelling in COMPRO was carried out using 3D higher order elements (20-noded ABAQUS internal elements) with three elements through the thickness of the laminate. The in-plane mesh for both models was identical. 4.2.1.2 Results 4.2.1.3 Deformation The deformation profile of the laminate along the length of the plate, which was manufactured using the Cure Cycle 1 (shown in Figure 4.14), is calculated and compared in both models 51 (COMPRO and LT elements) in Figure 4.15. The deformation is measured at the centre of the bottom surface of the plate along its length. It can be seen that the part has bent slightly after tool-removal, however, the deformations are very small because of the material symmetry. The figure also shows that the deformations calculated using COMPRO are identical to the results obtained using LT elements. 4.2.1.4 Transverse Normal Stress Comparison for Different Cure Cycles The transverse normal stress (𝜎𝑧) is calculated at the interface between of 0º and 90º plies along a path perpendicular and close to the free edge at the middle of the plate (𝑥 = 51𝑚𝑚). Figure 4.16 shows 𝜎𝑧 near the free edge after applying a tensile strain of 10% to the laminate for two different cure cycles. Also shown are the results for the case without cure modelling. The figure shows that the simulation, which includes process modelling and carries over the process-induced residual stresses to the in-service loading analysis, results in a higher transverse normal stress at the free edge. Indeed a significant amount of residual stress is neglected in the simulation without process modelling. The figure also shows that the Cure Cycle 2 introduces higher transverse stresses in the laminate at the free edge. 4.2.1.5 Free Edge Stresses Before and After Tool-Removal Laminated composite plates usually deform after tool-removal because of residual stresses developed in the material during the cure cycle. These deformations can either increase or relax the residual stresses. In the example shown in this section, the in-plane and out-of-plane stresses near the free edge are monitored before and after tool-removal. Figure 4.17 shows the bending average 𝜎𝑦 stresses in the bottom (0º ply) and middle layer (90º ply) for Cure Cycle 1. Since 52 CFRP laminas have higher shrinkage strain and CTE in the transverse directions than in the fibre directions, the top and bottom layers (0º plies) tend to contract in the y-direction more than the 90º ply due to the cure shrinkage and cool down process. Consequently, 𝜎𝑦 in the 90º ply is compressive while 𝜎𝑦 in the 0º ply is tensile. The positive value of 𝜎𝑦 in the 0º ply and its negative value in 90º ply confirm this. It can be seen that the 𝜎𝑦 stresses diminish as they approaches the free edge, which was expected because there are no traction on the free edge of the laminate. Figure 4.19 shows the inter-laminar shear stress at the interface of the bottom and middle layers of the laminate. As expected, the shear stress approaches zero at the free edge since there are no tractions on the free edge of the laminate. In addition, it can be seen that the shear stress distribution trend becomes a little smoother after tool-removal since the tool part interaction and the associated residual stresses in the laminate are relaxed. The transverse normal stress (𝜎𝑧) near the free edge for Cure Cycle 1 is calculated at the interface between 0º and 90º plies and presented in Figure 4.20. It can be seen that there is a high normal stress at the free edge, which asymptotically reduces as it moves away from the free edge. The figure also shows that the peel stress (𝜎𝑧) on the free edge was increased after the tool-removal. To obtain a better understanding about the reason of increase in σz at the interface of layers 1 (0º) and 2 (90º), 𝜎𝑧 is obtained through the thickness of the plate at the free edge and another point which is farther from the free edge (at 1mm distance from the free edge). As Figure 4.21 53 shows, the distribution of σz on free edge shows that there is a high compressive normal stress on the bottom surface of the plate at the free edge, which is due to the tool-part interaction. This compressive stress holds the inter-laminar stress lower before the tool-removal. Therefore, after the tool-removal, as the stress on the bottom surface of the plate diminishes the distribution of the 𝜎𝑧 changes, and 𝜎𝑧 is increased on the interface of layers 1 and 2 as it is shown in Figure 4.21. The figure also shows that the transverse normal stress goes to zero when it gets farther to the free edge. The transverse and in-plane stresses calculated using the LT element and 3D element show a very good agreement as shown in Figures 4.20, 4.21 and 4.22. These demonstrating the capability of the LT element in capturing the process induced stresses, both in-plane and out-of-plane. 4.2.1.6 Residual stress development during the curing process To control the process induced residual stresses in laminates; one should understand the source of residual stress evolution in the material. In this section the time history of residual stress at the free edge of a composite laminate is investigated. As the spatial trends of the stresses distribution near the free edge are the same at different stages of curing (it scales by different factors at different stages of curing), to obtain a better understanding of its development during the cure, they are measured and represented at the free edge, at different times during curing. 54 Figure 4.22 shows the development of the average in-plane stress, 𝜎𝑦 (calculated using Equation (4.1)), and maximum transverse normal stress, 𝜎𝑧, during the Cure Cycle 1 which are calculated at four specific times during the simulation (at 55, 175, 210 minutes and after tool-removal). The figure that shows the results of residual stresses does not change significantly during the first temperature ramp as the resin is not cured and the material is not stiff enough to develop residual stresses. During the temperature hold from (50 minutes to 175 minutes) the residual stress increases by a very small value, which is due to the cure shrinkage effect in the resin. During the cool down stage when the resin has already cured the residual stresses rise significantly. Based on this figure, the main portion of the final residual stress in the material is due to the cool-down process at the end of curing. The reason for this residual stress built-up is the fact that the material is gelled, vitrified and has developed enough stiffness that can induce stress in the material due to the mechanical strains introduced because of the CTE mismatch between the layers. The process modelling using the LT element has been verified for flat composite parts by comparing the residual stress development during the cure with COMPRO simulation results that employ ABAQUS 3D elements. 55 4.2.2 Curved Geometry In this section the curing of a curved composite part is modelled using the layerwise element in ABAQUS implicit and compared with the results for the same problem modelled using ABAQUS 3D elements. The composite part is made of Hexcel AS4/8552 CFRP with a thickness of 10mm cured under the cure cycle presented in Figure 4.24. The degree of cure calculated using COMPRO and layerwise element is also presented in this figure. Both unidirectional and cross-ply laminates are modelled. The material is cured on a rigid tool and the bottom surface of the part is fully bonded to the tool. Both layerwise and 3D elements are used for process simulation assuming the conditions, i.e. the part is fully constrained in the z-direction. The schematic of the part with dimensions and coordinate system are shown in Figure 4.23. The 3D (COMPRO) analysis is performed with 4 elements through the laminate thickness (radial direction), 30 elements in the tangential direction, and one element along the z-axis, while in the LT element analysis the same number of elements are used in tangential and z-direction, with a single element through the thickness which has four numerical layers. The 3D analysis is performed with 120 elements, while 30 elements were used in the LT element model 4.2.2.1 Results The warpage of the composite part is calculated after the tool-removal step using layerwise element and COMPRO 3D analysis. Figure 4.25 shows the deformation of the composite part with lay-up of [0º/90º]s after tool-removal. To make the deformation clear, displacements are 56 scaled by a factor of 10. It can be seen that the part springs in after tool-removal due to the different cure shrinkage and thermal strains in different directions. This phenomenon in the curing process of curved laminated composites has been shown using analytical modelling and experiments by Fernlund and Albert (Fernlund, 2005 and Albert et. al., 2002). The deformation diagram is represented as two different components of tangential (𝑢𝜃) and radial (𝑢𝑟) displacement as a function of 𝜃 in the cylindrical coordinate system (𝜃 is shown in Figure 4.23). Transverse and radial displacements at the bottom of the two laminate (𝑟 =50 𝑚𝑚) lay-ups are plotted in Figure 4.26 (for unidirectional 0º part) and in Figure 4.27 (for cross-ply [0º/90º]s part). The bending (circumferential) residual stress (𝜎𝜃) is calculated through the thickness for the cross-ply laminate at the 𝜃 = 45° section before and after tool-removal. 𝜎𝜃 distribution through the thickness, plotted in Figure 4.28, shows that the 0° plies are in compression while 90º plies are in tension. It also shows that 𝜎𝜃 is almost zero at the bottom layer of laminate before tool-removal, since it is confined by the tool, then after tool-removal the compressive stress increases in this layer. To verify the capability of LT element in predicting the out-of-plane process induced stresses, transverse normal stress (𝜎𝑟) distribution at 𝜃 = 45° is plotted in Figure 4.29 before and after tool removal. As the figure illustrates, the highest value of 𝜎𝑟 is on the 0º and 90° plies on top of the laminate. The positive sign of the transverse normal stress shows that the 0º and 90° plies have tensile normal at the interface between them. These types of stresses may lead to 57 delamination; however in this example the value of 𝜎𝑟 at the end of curing is not high enough to cause delamination. It can also be seen that the normal stress acting between the tool and the bottom of the laminate is reduced to zero after tool-removal. The inter-laminar tension at the interface between 0º and 90º plies at 𝜃 = 45° section is monitored from the beginning of the process all the way to the tool-removal stage. Evolution of 𝜎𝑟 during the curing process is plotted in Figure 4.30. The figure illustrates that the main part of the final residual inter-laminar stress is formed in the last step of curing during the cool-down phase since the resin is cured and stiff enough at this stage of manufacturing to result in appreciable residual stress in the material. As seen from the comparison of layerwise element and 3D continuum elements results in this section, the capability of layerwise element in capturing process induced deformations and stresses is verified against the COMPRO 3D modelling with multiple elements through the thickness of composite part. 4.3 Progressive Damage Modelling As discussed in Section 3.3, CODAM model is implemented in the layerwise element and linked to the process modelling. In this section the performance of CODAM model in the layerwise element without process modelling is investigated. Therefore, a numerical study of damage development and propagation in CFRP composite laminates is performed. In particular, the OCT tests were simulated in ABAQUS implicit. The OCT test was modelled using the implemented User Element (UEL) introduced in Section 3.4. 58 The OCT specimen geometry introduced by Kongshavn and Poursartip (Kongshavn and Poursartip, 1999) was selected to produce stable crack growth in compact tension specimens. Kongshavn demonstrated that crack growth in the OCT specimen with a sharp notch tip is stable under displacement control and is large enough so that the boundaries do not greatly affect the damage zone size or shape. The results of the OCT simulation using the LT element are compared with the results of modelling the same problem using CODAM in LS-DYNA. CODAM damage parameters are set to be identical in both ABAQUS and LS-DYNA simulations. Calibration of damage initiation strain is performed by Forghani (Forghani et al. 2011) and the value for fibre and matrix damage initiation strains are assumed to be 1.1% and 0.8% respectively. However the saturation strain for fibre and matrix damage can be estimated based on the fracture energy of the laminate measured in the experiments. The OCT virtual specimen geometry and dimensions are shown in Figure 4.31. The laminate thickness is 1mm and the size of the elements in the damage area is 1𝑚𝑚 × 1𝑚𝑚. Figure 4.32 shows the rectangular are (51𝑚𝑚 × 30𝑚𝑚) the damage area, in which the FE mesh is uniform with smaller elements with the objective of predicting damage propagation more accurately. Both unidirectional and quasi-isotropic lay-ups are used in the simulations carried out in this chapter and the material is Hexcel AS4/8552 CFRP. 59 Due to the fact that solution convergence of the progressive damage modelling is ABAQUS implicit is difficult, the minimum allowable time increment is set to be smaller than its default value in ABAQUS. The default minimum allowable time increment in ABAQUS is 1.0 ×10−5𝑠, whereas it is set to be Δ𝑡𝑚𝑖𝑛 = 1.0 × 10−15𝑠 in this study. The number of allowable cutbacks (subincrements) for an increment in ABAQUS solver is set to 15 from its default value which is 5. Cutback is the process of decreasing time increment size automatically by ABAQUS solver with the objective of achieving. When ABAQUS solver fails to solve the equilibrium equation in an increment it decreases the increment size to one fourth and re-starts the iterations to solve the problem. 4.3.1 Unidirectional laminate For the first example to check the capability of the UEL in progressive modelling, an OCT test of a unidirectional 90º laminate (fibres oriented parallel to the notch direction) is modelled using the UEL. Fibre and matrix damage saturation strain is assumed to be 130% in this example. The saturation strain for unidirectional laminate modelling is not calibrated with the experiments; it is simply set to an arbitrarily large value to make the specimen failure more gradual in order to achieve convergence and to observe the deformation details of the model. In the first attempt of OCT simulation the ABAQUS implicit solver failed to converge in the first few steps after damage initiation. As Figure 4.33 shows the area of the damaged elements (matrix damage) occurred in a large area around the crack tip which may be the reason for solution divergence. Thus the model was modified to have only one row of elements along the crack tip and the crack was forced to grow in a horizontal line along the crack tip. This 60 assumption is valid for a symmetric problem so that it could be used for a 90° unidirectional laminate modelling. The simulation was repeated with a band of elements along the notch. In other words, only one row of elements along the notch was damaged, while other elements were behaving as elastic material. Although the analysis goes further than previous case, the solution failed to converge after a few number of elements were damaged (Figure 4.34). To investigate the convergence problem, the same OCT test was modelled using Simple Damage Model (SDM) implemented by (McClennan, 2004) as a User Material in ABAQUS implicit for 2D element. The modelling was carried out with both linear (4-noded) and quadratic (8-noded) elements. Figure 4.35 shows the result of the OCT test simulation for both linear and higher order (quadratic) elements. The figure illustrates that in the modelling with linear elements the solution converged to the end of the complete failure of the specimen while the solution diverged after couple of elements damaged in the analysis with higher order elements. As the UEL is a higher order element the same reason may have caused divergence for both CODAM in the UEL and SDM in higher order elements modelling. The ABAQUS’s solution procedure which is used in this problem is explained as follows with the objectives of solving the equilibrium problem. ABAQUS’s default implicit solver, the direct sparse solver, is used in this study. The direct sparse solver uses Newton’stechnique for solving the nonlinear equation. The following equation shows the process of solving the equations using Newton’smethod: 61 𝑲𝒖− 𝑭 = 𝑹 (4.2) where 𝑲 is the stiffness matrix, 𝑭 is the external load vector and 𝑹 is the residual (of equilibrium) load vector. For the exact solution of Equation (4.2) the residual force vector should be zero. Therefore the solver attempts to find the solution by minimizing the norm of the residual force vector in multiple iterations. When any of the residual force entries is greater than the specified minimum error after each attempt, the solver reduces the increment size and re-solves the equation. In case the increment size becomes less than the minimum allowable or number of attempts becomes greater that the allowable number of attempts, the solver aborts the iteration. In such a case, the node with the highest residual force is usually the node that makes the solution unstable and causes convergence failure. For the OCT simulation the nodal residual forces were examined for the UEL simulation and the highest residual force was for the middle node of an element on the crack band undergoes damage. Figure 4.36 shows the behaviour of the middle node for the crack band elements which are damaged (the displacements are scaled by a factor of 20). The figure shows that the middle nodes of the crack band elements have an unexpected behaviour; some nodes have positive displacement in y-direction and some have negative displacement. In the higher order elements there are two integration points located above and below of the middle node. During the crack growth in the material, due to symmetry both integration points damaged fully at the same time and the material stiffness becomes zero on both sides of the middle node so that the node becomes unstable and consequently the solution fails to converge. 62 To solve the aforementioned convergence issue, based on the symmetric geometry of the problem, a symmetric boundary condition was applied to the problem by fixing the middle nodes of the crack band elements in the y-direction. Figure 4.37 shows the OCT specimen at the end of the simulation after applying symmetric boundary condition on the middle nodes (i.e. setting the y-direction displacement to zero). It can be seen that the crack grows along a horizontal line along the crack tip. The applied load is plotted versus the crack mouth opening displacement (CMOD) in Figure 4.38. There were no convergence problems in this simulation as it continued successfully to the complete failure of the specimen. 4.3.2 Quasi-isotropic Laminate In this section the OCT test of a quasi-isotropic laminate with the lay-up of [90/45/0/-45]s is modelled using CODAM material model implemented in the UEL in ABAQUS implicit and the results are compared with the simulation of the same problem using CODAM in LS-DYNA. Unlike the previous problem, the fibre and matrix damage saturation strains were calibrated with the experiments of Zobeiry (Zobeiry, 2010). 4.3.2.1 Saturation Strain Calculation As mentioned before the damage saturation strain is calibrated based on the fracture energy obtained from the physical OCT test performed by (Zobeiry,2010) with the same lay-up, material, and dimensions. The fracture energy calculated for this experiment was 𝐺𝑓 = 110𝐾𝐽𝑚2. 63 Equation (4.3) shows the relation between fracture energy (𝐺𝑓) and energy density of the element (𝛾𝑠). As the element height in this example is ℎ𝑒 = 1𝑚𝑚 the element energy density should be 𝛾𝑠 = 110 𝑀𝑃𝑎 (𝑀𝐽𝑚3). 𝐺𝑓 = ℎ𝑒 . 𝛾𝑠 (4.3) To find the corresponding damage saturation strain four one-element unidirectional strain tests were carried out for the [90/45/0/-45]s lay-up with different saturation strains (15%, 20%, 25% and 30%). The equivalent stress of the element is plotted versus equivalent strain in Figure 4.39. The one-element test results are compared with the result of modelling the same one element test in LS-DYNA explicit solver, using the built-in material in which CODAM model has been implemented (MAT219). The figure shows that the behaviour of one element in these two models is the same. In addition, the area under the stress-strain curve represents the energy density of the element and shows that the higher the saturation strain the larger the energy density. The energy density is plotted versus damage saturation strain in Figure 4.40. It shows that the saturation strain of 20% corresponds the energy density of 110 𝑀𝑃𝑎. 4.3.2.2 Simulation Results Figure 4.41 shows a comparison of the results of the OCT simulation using LS_DYNA explicit and LT element within ABAQUS implicit. It shows that the results from the LT element analysis in ABAQUS implicit agrees well with the results obtained from LS-DYNA explicit analysis from the beginning of the loading until the load reaches to the peak value, and the implicit solver 64 fails to converge. This figure also shows that the results obtained from the ABAQUS implicit are smoother than the explicit analysis results, while LS-DYNA explicit solver is able to predict the softening behaviour of the specimen better than the ABAQUS implicit solver. Another issue with the ABAQUS implicit analysis is that it was not able to simulate the specimen behaviour to ultimate failure (i.e. the softening regime in the load-displacement curve). The solution has been failed to converge after the peak load due to the convergence difficulties of the implicit solver. Since solution convergence in damage modelling improves by having a smoother or more gradual softening behaviour, the problem was modelled with higher damage saturation strain. Figure 4.31 shows the load-CMOD diagram of the two OCT test simulations (damage saturation strain values of 20% and 250%), and it confirms that the simulation with a higher damage saturation strain value (𝜖𝑠 = 250%) can be successfully solved to the end of simulation, where the material reaches ultimate failure. It confirms that the OCT test with smoother softening behaviour has a better chance of convergence compared to the other simulations that correspond to a more sudden and steeper softening response. 65 4.4 Tables Table 4.1 Mechanical properties of Steel Parameter Value Unit Elasticity Modulus (𝐸) 200 𝐺𝑃𝑎 Poisson Ratio (𝜈) 0.3 Shear Modulus (𝐺) 79.3 𝐺𝑃𝑎 Table 4.2 Mechanical properties of fully cured AS4/8552 FRP Parameter Value Unit 𝐸1 122.6 𝐺𝑃𝑎 𝐸2 9.72 𝐺𝑃𝑎 𝜈12 0.268 𝜈21 0.021 𝜈23 0.47 𝐺12 5.21 𝐺𝑃𝑎 𝐺23 3.38 𝐺𝑃𝑎 66 4.6 Figures Figure 4.1 Schematic of the flat homogenous, isotropic plate with material properties of steel under end shear load Figure 4.2 Deflection of the bottom plane of the plate shown in Figure 4.1 calculated using the LT element with two different meshes compared with the analysis results using the ABAQUS 3D continuum elements result L= 50mm t = 1 mmP = 1 N/mm2xyz-2.5-2-1.5-1-0.500 10 20 30 40 50Deflection (mm) x (mm)ABQUS 3D ElementLayerwise - Fine meshLayerwise - Coarse mesh67 Figure 4.3 Through-thickness distributionofσx obtained at the centre of the plate (x=25 mm, y=5mm) using the LT element with two different meshes compared with the analysis results using the ABAQUS 3D continuum elements result 00.20.40.60.81-200 -100 0 100 200Z (mm) σx(MPa)ABQUS 3D ElementLayerwise - Coarse meshLayerwise - Fine mesh68 (a) (b) Figure 4.4 Through-thickness distribution of the transverse shear stress τxz obtained at all in-plane integration points of the LT elements with (a) full and (b) reduced integration scheme. Also, in-plane integration points numbering for a typical LT element is shown. The average values compared with the results of ABAQUS 3D continuum elements modelling for the plate-bending problem shown in Figure 4.1. The average values agree with the 3D modelling of the problem. 00.20.40.60.81-6 -4 -2 0 2 4Z (mm) τxz(MPa)Integration points1 and 384,67,952ABAQUS 3D ElementsAverage00.20.40.60.81-2 -1 -1 0 1Z (mm) τxz(MPa)Integration points1 and 3ABAQUS 3D ElementsAverageIntegration points2 and 4123 4xy69 Figure 4.5 Transverse shear stress τxz distribution through the thickness of the plate shown in Figure 4.1 calculated at the centre point using both coarse and fine mesh of LT elements with reduced integration scheme compared to the corresponding ABAQUS 3D element modelling result. Figure 4.6 Schematic of a laminated plate with [0/90/0] layup under an end shear load 00.20.40.60.81-1.5 -1 -0.5 0Z (mm) τxz(MPa)ABQUS 3D ElementLayerwise - Coarse meshLayerwise - Fine meshL = 50 mmt = 1 mmP = 1 N/mm20 90 xyz70 Figure 4.7 Deflection of the bottom surface of the laminated plate shown in Figure 4.6 using LT elements and compared with ABAQUS 3D continuum elements Figure 4.8 Through-thickness distribution of the axial stress σx measured at the centre of the plate (x=25mm, y=5mm) -5-4-3-2-100 10 20 30 40 50Deflection (mm) x (mm)ABAQUS 3D ElementLayerwise element0.00.10.20.30.40.50.60.70.80.91.0-200 -100 0 100 200z (mm) σx(MPa)ABAQUS 3D ElementLayerwise element71 Figure 4.9 Transverse shear stress distribution through the thickness of the laminated plate calculated at the centre point of the plate using reduced integrated LT element and compared with ABAQUS 3D element modelling results Figure 4.10 Transverse shear strain distribution through the thickness of the laminated plate calculated at the centre point of the plate using reduced integrated LT element and compared with ABAQUS 3D element modelling results 0.00.20.40.60.81.0-2 -2 -1 -1 0z (mm) xz(MPa)ABAQUS 3D ElementLayerwise element0.00.20.40.6.81.0-6.E-04 -4.E-04 -2.E-04 0.E+00z (mm) ϵxz(mm/mm)ABAQUS 3D ElementLayerwise element72 (a) (b) Figure 4.11 Schematic of a laminated plate with (a) [0/90/0] and (b)[90/0/90] layup under axial load. The free edge peel stress σz is calculated near the free edge using LT element and ABAQUS 3D continuum element. L = 50 mmt = 1 mm0 90 xyzu = 0.25 mm L = 50 mmt = 1 mmu = 0.25 mm 0 xyz90 Free EdgeFree Edge73 (a) (b) Figure 4.12 Transverse normal inter-laminar stress distribution at an interface between 0° and 90° layers near the free edge calculated using the LT element and compared with ABAQUS 3D continuum element modelling for (a)[0/90/0] and (b)[90/0/90] layups -40-200204060801001201400 1 2 3σz(MPa) Y (mm)ABAQUS 3D ElementLayerwise element-30-25-20-15-10-505100 1 2 3σz(MPa) Y (mm)ABAQUS 3D elementLayerwise element74 Figure 4.13 Schematic of a cross-ply composite plate being cured on a rigid invar tool Figure 4.14 Two different cure cycles used for curing of the laminate shown in Figure 4.13 102 mmToolPart0 directionx (u)z (w)y (v)040801201602000 200 400 600Temperature(°C)Time (min)Cure Cycle 1Cure Cycle 275 Figure 4.15 Deformation profile after tool-removal for the laminate shown in Figure 4.13 undergoing Cure Cycle number 1 shown in Figure 4.14 Figure 4.16 Transverse normal stress calculated at the interface between layer 1 (0°) and layer 2 (90°) after tool-removal and subsequent uniaxial loading for two different cure cycles. Also the results for the case without cure modelling are shown. -0.09-0.07-0.05-0.03-0.010.010 20 40 60 80 100Deformation (mm)x (mm)Layerwise elementCOMPRO-40-2002406081000 0.2 0.4 0.6 0.8 1σz(MPa)Distance from the edge (mm)Cure Cycle 1Cure Cycle 2Without Cure Modelling76 Figure 4.17 Average stress on layer 1 (bottom 0° ply) and layer 2 (90° ply) distribution along the path near the free edge before and after tool-removal, compared with COMPRO 3D modelling Figure 4.18 Average σy stress evaluated the bottom layer after tool-removal calculated using COMPRO 3D element and LT element modelling for the plate shown in Figure 4.13 under Cure Cycle 1 -120-80-40040800 0.2 0.4 0.6 0.8 1σy(MPa)Distance from the edge (mm)Layer 1 - Before tool removalLayer 1 - After tool removalLayer 2 - Before tool removalLayer 2 - After tool removal400 0.2 0.4 0.6 0.8 1σy(MPa)Distance from the edge (mm)Layerwise elementCOMPRO77 Figure 4.19 Inter-laminar transverse shear stress τyz at the interface between layers 1 (0°) and 2 (90°) y-axis (with the origin at the free edge) before and after tool-removal, compared with COMPRO 3D results -20-16-12-8-400 0.2 0.4 0.6 0.8 1τyz(MPa)Distance from the edge (mm)Before tool removal - COMPROAfter tool removal - COMPRO Before tool removal - Layerwise After tool removal - Layerwise78 Figure 4.20 Transverse normal stress (σz ) at the interface between layers 1 (0°) and 2 (90°) y-axis (with the origin at the free edge) before and after tool-removal, compared with COMPRO 3D results -10010203040500 0.2 0.4 0.6 0.8 1σz(MPa)Distance from the edge (mm)Before tool removal - COMPRO After tool remova-COMPRO Before tool removal - Layerwise After tool removal - Layerwise79 Figure 4.21 Transverse shear stress distribution through the thickness of laminate on the free edge and at a point far from the free edge before and after tool-removal. σz is increased on the free edge and also it is reduced after tool-removal at the bottom of the plate. 00.380.761.14-40 -20 0 20 40 60Z(mm)σz(MPa)Before tool removal - on the free edge After tool removal - on the free edgeBefore tool removal - far from the edge After tool removal - far from the edge80 Figure 4.22 Residual stresses evolution during the cure cycle, compared with COMPRO 3D modelling Figure 4.23 Schematic of the curved part with layup of [0/90]s made from AS4/8552 CFRP cured on a rigid tool. -30-20-10010203040040801201602000 50 100 150 200 250σ(Mpa)T (°C)Time (min)Temperature z y - layer 1 y - layer 2σσσTool removalzrθ10 mmR=50 mmθ81 Figure 4.24 The one hold cure cycle is used for the curved part curing and the resin degree of cure during the processing simulation in COMPRO and Layerwise element Figure 4.25 Curved part deformation after tool-removal in r-θplane.Tomakethedeformation visible in the figure, displacements are scaled by a factor of 10 0408012016020000.20.40.60.810 100 200 300Temperature (°C)Degree of CureTime (min)Layerwise elementCOMPROTemperature82 (a) (b) Figure 4.26 (a)Radial and (b)tangential displacement of the curved part at r = 50mm calculated from the LT element modelling compared with COMPRO 3D modelling for a unidirectional part ([0] lay-up) -0.4-0.3-0.2-0.10.00.10.20 20 40 60 80 100Ur (mm)θLayerwise elementCOMPRO-0.06-0.04-0.020.000.020.040.060 20 40 60 80 100Uθ(mm)θLayerwise elementCOMPRO83 (a) (b) Figure 4.27 (a)Radial and (b)tangential displacement of the curved part at r = 50mm calculated from layerwise element modelling compared with COMPRO 3D modelling for a cross-ply part ([0/90]s lay-up) -0.6-0.4-0.20.00.20.40 20 40 60 80Ur (mm)θLayerwise elementCOMPRO0.000.050.100.150 20 40 60 80Uθ(mm)θLayerwise elementCOMPRO84 Figure 4.28 In-plane stress distribution through the thickness of the curved [0°,90°]s part at θ=45 before and after tool-removal calculated using layerwise element modelling compared with COMPRO 3D modelling Figure 4.29 Transverse normal (or radial) stress distribution through the thickness of the curved [0°,90°]s part at θ=45 at different curing stages calculated from layerwise element modelling compared with COMPRO 3D modelling 0.02.04.06.08.010.0-200 -150 -100 -50 0 50 100Thickness coordinate (mm) σθ(MPa)Before tool-removalAfter tool-removalLayerwiseCOMPRO0°0°90°90°0.02.04.06.08.010.0-4 -2 0 2 4 6Thickness coordinate (mm) σr(MPa)After tool-removalBefore tool-removalLayerwiseCOMPRO0°0°90°90°85 Figure 4.30 Inter-laminar transverse normal stress variation at the interface between the top 0º and 90º plies at θ=45 section of the curved part during the cure cycle obtained from modelling with layerwise element and compared with COMPRO 3D modelling -10123456040801201602000 100 200 300 400σr(MPa)T (°C)Time (min)Temperature - COMPRO - Layerwise elementσrTool removalσr86 Figure 4.31 Schematic figure of the OCT test modelled using layerwise element with damage modelling capability in ABAQUS implicit. Figure 4.32 A rectangular area (51mm×30mm) around the notch tip, which has a uniform fine mesh with element size of 1mm×1mm W=80 mma = 33 mmd=19.1 mmc=38.7 mm100 mm100 mm80.65 mm0º fibresxy30 mm51 mm87 Figure 4.33 The FE model of an OCT test on a 90° unidirectional lamina made of AS4/8552 CFRP. The modelling with higher order layerwise element failed to converge. The red circle shows the area with damaged elements. Figure 4.34 Simulation of OCT where the damage zone is limited to one row of elements 88 (a) (b) Figure 4.35 OCT test modelled using CODAM in ABAQUS implicit wit 2D (a) linear (4-noded) and (b) quadratic (8-noded) elements. The linear element modelling successfully converged to the end of complete specimen failure however the modelling with the higher order elements failed to converge after damage occurred in the first element 89 Figure 4.36 The middle node of damaged elements show an unexpected behaviour. The highest residual force (which makes the results diverge) occurs on the middle node of the element which is about to get damaged because the stiffness at both integration points above and below the middle node goes to zero. Figure 4.37 The symmetry of the problem is applied to the model by fixing the middle node of the crack-band elements. The solution is converges successfully until complete failure of the unidirectional (90º) lamina. 90 Figure 4.38 The load applied load versus CMOD in the OCT test simulation of a unidirectional (90º) composite plate in ABAQUS implicit using the higher order layerwise element. Figure 4.39 Stress vs. strain from unidirectional strain simulation of one-element to find the fracture energy for different values of fibre and matrix saturation strains. 00.20.40.60.80 1 2 3Load (KN)CMOD (mm)02004006008000 20 40 60Stress (Mpa)Strain (%)15% - LS-Dyna20% - LS-Dyna25% - LS-Dyna30% - LS-Dyna15% - Layerwise20% - Layerwise25% - Layerwise30% - Layerwise91 Figure 4.40 Energy density vs. damage saturation strain plotted to find the fibre and matrix damage saturation strains which corresponds to the energy density of γs=110 MPa Figure 4.41 The OCT test simulation of [90/45/0/-45]s laminate made from Hexcel AS4/8552 CFRP using layerwise element in ABAQUS implicit comparing with simulation of the same laminate using an explicit solver in LS-DYNA. 6010014018010 20 30 40Energy Density(MPa)Saturation Strain (%)00.511.522.530 0.5 1 1.5 2 2.5Load (KN)CMOD (mm)ABAQUS (LT element)LS-DYNA92 Figure 4.42 The OCT test simulation of [90/45/0/-45]s laminate made from AS4/8552 using layerwise element in ABAQUS implicit with to different damage saturation strain parameter for CODAM. It shows the behaviour of the model with higher saturation strain is smoother and consequently its chance of convergence is higher than the one with lower damage saturation strain. 01234560 2 4 6 8Load (KN)CMOD (mm)Saturation Strain = 250%Saturation Strain = 20%93 Chapter 5: Numerical Examples 5.1 Effect of Residual Stresses on Damage Initiation of Laminated Plates In this section a simple numerical example is presented to demonstrate the effect of process-induced residual stresses on promoting the initiation of in-plane damage in a cross-ply composite laminate subjected to a uniaxial tensile load. To show the effect of considering residual stresses in the in-service analysis of material, a composite laminated plate was modelled subjected to two different cure cycles. The results of these two analyses are compared with each other as well as with the results of modelling the laminate without considering any process-induced residual stresses in the material. The problem was modelled using a layerwise element UEL implemented in ABAQUS implicit. 5.1.1 Problem Description A cross-ply, CFRP composite laminate with a span of 102mm, a width of 9.75mm and a thickness of 1.14mm is considered for this study. The plate is made of Hexcel AS4/8552 material with a [0/902]S lay-up. First the composite plate was virtually cured on a rigid invar tool (Figure 5.1 - (a)) in an autoclave using two different cure cycles as shown in Figure 5.2. Then the cured plates were removed from the tool and subjected to an applied uniform in-plane displacement as shown in Figure 5.1 - (b). Due to symmetry, only half of the plate was modelled. 94 5.1.2 Results At the end of the cure cycle, the process-induced stresses in the transverse (90 degree) ply in the longitudinal direction (i.e. loading direction) of the laminate cured with a 1-hold and 2-hold cure cycles (shown in Figure 5.2) are 38 MPa and 42 MPa, respectively. With the damage initiation strain for fibre and matrix assumed to be 1.1% and 0.8%, respectively, no damage was initiated during the curing process. To better illustrate the effect of matrix damage on the stiffness of the laminate and locate the point where damage initiates, the derivative of the applied force with respect to displacement as a function of the applied displacement is plotted in Figure 5.3. The figure compares the damage initiation of an as-designed part with that of the as-manufactured part that includes the effect of process-induced residual stresses due to the two different cure cycles. As shown in the figure, the damage initiation starts much earlier in the as-manufactured laminate that accounts for the process-induced stresses than that in the as-designed laminate. The figure also shows that the damage initiation may significantly vary based on the cure cycle used for manufacturing. 5.2 Curved-Part Delamination This example explores the effect of curing process on the formation of inter-laminar stresses in composite laminates. An L-shaped part is cured and then subjected to a loading until delamination occurs in the part. The idea is to compare the state of inter-laminar stresses leading to onset of delamination in presence and absence of process-induced stresses. 95 5.2.1 Problem Description Curved laminated composite parts are encountered in many composite structures such as angle bracket, co-cured web or frame, and one of the final failure modes in such structures is delamination in the curved part. Inter-laminar stresses which result in delamination can be introduced during the curing process and/or in-service loading of the material. There have been many studies on the curing of thick curved composite laminates by focusing on different aspects of cure modelling. Also many studies have been carried out on the delamination of curved parts such as L-Shaped laminates. However the focus of this research is on linking these two simulations within a common numerical framework. To demonstrate this, the problem is simplified in some of the details with respect to cure and delamination modelling. The simulation of the curing stage and delamination initiation of a curved part are conducted within a common framework. To investigate the inter-laminar tensile stress and delamination of curved parts, the same curved part studied in Section 4.2.2 was modelled. The laminate layup is [0/90]s and its dimensions are shown schematically in Figure 5.4. The problem was modelled in two stages. During the first stage the CFRP material was cured on a tool (either convex or concave). Then, after tool-removal, a uniform displacement was applied to both ends of the part and continued until delamination initiated in one of the 3 layer interfaces shown in Figure 5.5. The criterion for initiation of delamination in the part was exceedance of the inter-laminar tensile strength. Secondly, when the tensile stress at these interfaces reach to the damage initiation stress of the resin, it is assumed that delamination initiates at that interface. The damage initiation stress value 96 in the resin is the corresponding stress when the damage initiation strain is reached. This damage initiation strain followed the value used by Forghani (𝜖𝑖 = 0.8%) (Forghani 2011). The simulation results of the model are compared with the modelling of the same curved laminate under displacement applied to its both ends without any process modelling. 5.2.2 Results To view the stress distribution in the part under the applied displacement before delamination, a displacement of 1mm is applied on both ends of the part as shown in Figure 5.5 and the out-of-plane stress (𝜎𝑟) is plotted in Figure 5.6. The transverse normal stress distribution in the part shows that the highest inter-laminar tensile stress occurs at Interface 1 and at the centre of the part(𝜃 = 45°). The transverse normal stress distribution along the circumferential of the part is also plotted at 3 layer interfaces in Figure 5.7. The maximum tensile stress occurs at Interface 1 and consequently delamination is likely to initiate at that location. The curing process of the laminate was simulated using the cure cycle shown in Figure 5.8 on a rigid tool and the part was considered to be fully bonded to the tool during the complete curing process. Schematics of the two different tools, used for the curing simulation, are shown in Figure 5.9. After the curing simulation and tool-removal, a uniform displacement was applied to both ends of the part the magnitude of which continued until damage initiated in the Interface 1 at the centre point. As Figure 5.10 shows, the inter-laminar tensile stress is lower for a given applied 97 displacement when the process modelling is included in the simulation. Based on the results of Section 4.2.2, the interlaminar stress after the curing process is negative which means that the inter-laminar stresses are in compression. While in the modelling without process simulation this stress is assumed to be zero before displacements are applied to the part. It can be seen that by linking process modelling to damage modelling and considering process-induced residual stresses in the structural analysis, the delamination capacity of the laminate during its in-service life is different from that which ignores the process-induced stresses. In this case the laminate was found to have a slightly higher capacity. Comparing the results of simulations using convex and concave tools shows that there is no significant difference between these two types of tooling geometries. The main reason of having same results in both conditions is that in the material model used for AS4/8552 in COMPRO the material modulus development occurs in a very short time period, during which there is no significant deformation due to the thermal strains and cure shrinkage strain. Also according to the results of Section 4.2.2, most of the residual stress is built-up during the cool down stage, where the material is a fully elastic solid. Therefore the residual stress caused due to the cool down process is only a function of the final boundary condition of the part (after tool-removal). In other words, since most of the residual stress developed when the material modulus was constant, having different boundary conditions does not introduce different final residual stress in the material. 98 5.3 Figures (a) (b) Figure 5.1 Schematic of a [0/90]s laminate made from AS4/8552 CFRP : (a) cured on a rigid invar tool in autoclave using two different cure cycles and (b) the laminate is subsequently removed from the tool and subjected to an external uniform in-plane displacement. 102 mmToolPart0º fibrex (u)z (w)y (v)uu99 Figure 5.2 Two different cure cycles used for curing of the laminate in Figure 5.1 Figure 5.3 Laminate stiffness vs. applied displacement for two different cure cycles and without manufacturing simulation. Matrix damage initiation occurs with lower applied displacement for the ones with cure modelling because of the residual stresses in the matrix. 040801201602000 200 400 600Temperature(°C)Time (min)Cure Cycle 1Cure Cycle 2Matrix damage initiation pointFibre failure point100 Figure 5.4 Schematic figure of the curved-part modelled for detecting delamination Figure 5.5 A displacement applied to both ends of the curved-part until delamination occurs in one of the plies interfaces zrθ10 mmR=50 mmθuuInterface 1Interface 2Interface 3101 Figure 5.6 Transverse normal stress distribution in the curved part when a displacement of 1mm is applied to both ends without process modelling. The stress has not reached the critical value for delamination σcr=79.4MPa (which corresponds 0.8% strain in matrix). The highest stress occurs at Interface 1 and at section θ=45º. Figure 5.7 Transverse normal stress on 3 potential delamination interfaces when a displacement of 1mm is applied to both ends of the curved part without process modelling. -400-300-200-1000100200300σr (Mpa)-40-20020400 20 40 60 80σr (MPa)θInterface 1Interface 2Interface 3102 Figure 5.8 The temperature cycle applied to the cuved part for process modelling (a) (b) Figure 5.9 The curved-part is modelled with two different (a) convex and (b) concave rigid tools. 040801201602000 100 200 300Temperature (ºC)Time (min)Rigid ToolRigid Tool103 Figure 5.10 Transverse normal stress at Interface 1 and at section θ=45ºas a function of the applied displacement. While the figure illustrates the modelling with curing simulation shows lower value of the residual stress in the material, there is no significant difference between the results for convex and concave tools. -3003060900 0.5 1 1.5 2 2.5 3σr(Mpa)Displacement (mm)w/o cure modellingConvex toolConcave tool104 Chapter 6: Conclusion This thesis presented a modelling framework for seamless simulation of composite laminates from the beginning of the curing process through the in-service behaviour and failure of the structure, implemented with the Layerwise Theory (LT) elements. The conclusions based on the presented work in this thesis are as follows: Curing process and in-service loading of the material in the elastic region using the LT elements with one element through the thickness, results in very accurate prediction of deformation and residual stresses comparable to accuracy of the predictions made using multiple ordinary 3D elements through the thickness. Unlike Classical Laminated Plate Theory (CLPT) and 3D FEM modelling with single element through the thickness of the laminate, LT is capable of capturing transverse shear stresses and strains distribution through the thickness for both flat and curved geometries. It was shown that existence of cure-induced residual stresses can significantly affect the initiation of damage in composite laminates, and ignoring residual stresses may lead to over-prediction of the strength and unsafe design. Due to extreme non-linear nature of damage propagation in laminated composites, implicit FE solver faces difficulties with convergence of the nonlinear solution of equilibrium equation. The convergence of the solution becomes more difficult when the material response has a steeper softening response during the damage propagation (i.e. less fracture energy). Moreover, the analysis results shows that the higher the number of layers in the LT element, the lower the chance of convergence in the damage simulation. 105 Bibliography Albert, Carolyne, and Göran Fernlund. "Spring-in and Warpage of Angled Composite Laminates." Composites Science and Technology 62.14 (2002): 1895-912. Arafath, A. A., Reza Vaziri, and Anoush Poursartip. "Modelling Process-induced Deformations in Composite Structures Using Higher Order Elements." International Conference of Composite Materials (ICCM). 2007. Arafath, Abdul Rahim Ahamed. "Efficient Numerical Techniques for Predicting Process-Induced Stresses and Deformations in Composite Structures." PhD Thesis, Deparment of Civil Engineering, the University of British Columbia, Vancouver, BC, Canada. (2007) Arciniega, Roman A., and JN Reddy. "Consistent Third-Order Shell Theory with Application to Composite Cylindrical Cylinders." AIAA Journal 43.9 (2005): 2024-38. Barbero, EJ, and JN Reddy. "Modeling of Delamination in Composite Laminates using a Layer-Wise Plate Theory." International Journal of Solids and Structures 28.3 (1991): 373-88. Beuth, Jack L., and SH Narayan. "Residual Stress-Driven Delamination in Deposited Multi-Layers." International Journal of Solids and Structures 33.1 (1996): 65-78. Bogetti, Travis A., and John W. Gillespie. "Process-Induced Stress and Deformation in Thick-Section Thermoset Composite Laminates." Journal of Composite Materials 26.5 (1992): 626-60. Bogetti, Travis A., and John W. Gillespie. "Two-Dimensional Cure Simulation of Thick Thermosetting Composites." Journal of Composite Materials 25.3 (1991): 239-73. Chen, PC, and RL Ramkumar. "RAMPC-an Integrated Three-Dimensional Design Tool for Processing Composites". IN: Materials-Pathway to the future; Proceedings of the Thirty-third International SAMPE Symposium and Exhibition, Anaheim, CA, Mar. 7-10, 1988 (A88-42326 17-23). Covina, CA, Society for the Advancement of Material and Process Engineering, 1988, p. 1697-1708. 1988. 1697-1708. Cook, Robert D. Finite Element Modeling for Stress Analysis. Wiley, 1994. 106 Crossman, FW, and ASD Wang. "The Dependence of Transverse Cracking and Delamination on Ply Thickness in graphite/epoxy Laminates." Damage in composite materials, ASTM STP 775 (1982): 118-39. Ersoy, Nuri, et al. "Modelling of the Spring-in Phenomenon in Curved Parts made of a Thermosetting Composite." Composites Part A: Applied Science and Manufacturing 41.3 (2010): 410-8. Fernlund, G., et al. "Finite Element Based Prediction of Process-Induced Deformation of Autoclaved Composite Structures using 2D Process Analysis and 3D Structural Analysis." Composite Structures 62.2 (2003): 223-34. Fernlund, Göran. "Spring-in of Angled Sandwich Panels." Composites Science and Technology 65.2 (2005): 317-23. Forghani, A., et al. "A Non-Local Approach to Simulation of Damage in Laminated Composites." Damage in Composites: DEStech Publication Inc. (2012): 279. Forghani, Alireza. “A Non-Local Approach to Simulation of Damage in Composite Structures” PhD Thesis, Department of Civil Engineering, the University of British Columbia, Vancouver, BC, Canada. (2011) Garg, Amar C. "Delamination—a Damage Mode in Composite Structures." Engineering Fracture Mechanics 29.5 (1988): 557-84. Gong, X. Lu, et al. "Residual Stress Distribution and its Influence on the Mechanical Behaviour of Composite Laminates". International Conference of Composite Materials (ICCM). 1999. Hubert, Pascal, and Anoush Poursartip. "A Review of Flow and Compaction Modelling Relevant to Thermoset Matrix Laminate Processing." Journal of Reinforced Plastics and Composites 17.4 (1998): 286-318. Jeronimidis, G., and AT Parkyn. "Residual Stresses in Carbon Fibre-Thermoplastic Matrix Laminates." Journal of Composite Materials 22.5 (1988): 401-15. 107 Johnston, Andrew A. “An integrated model of the development of process-induced deformation in autoclave processing of composite structures” PhD Thesis, Department of Materials Engineering, the University of British Columbia, Vancouver, BC, Canada, . (1997) Johnston, Andrew, Reza Vaziri, and Anoush Poursartip. "A Plane Strain Model for Process-Induced Deformation of Laminated Composite Structures." Journal of Composite Materials 35.16 (2001): 1435-69. Johnston, Pascal Hubert Andrew, Reza Vaziri, and Anoush Poursartip. "A Two-Dimensional Finite Element Processing Model for FRP Composite Components". Proceedings of the Tenth International Conference on Composite Materials: Processing and manufacturing. Woodhead Publishing , 1995. 149. Kim, RY, and HT Hahn. "Effect of Curing Stresses on the First Ply-Failure in Composite Laminates." Journal of Composite Materials 13.1 (1979): 2-16. Krajcinovic, Dusan, and Kaushik Mallick. "Micromechanics of the Progress Induced Damage Evolution in Thermosets." Journal of the Mechanics and Physics of Solids 43.7 (1995): 1059-86. Li, Chun, Michael R. Wisnom, and L. Graeme Stringer. "Simulation of Process-Induced Residual Stresses in Thick Filament Wound Tubes." Simulation and Modeling (2002): 271-82. Loos, Alfred C., and George S. Springer. "Curing of Epoxy Matrix Composites." Journal of Composite Materials 17.2 (1983): 135-69. Martin, Roderick H., and Wade C. Jackson. "Damage Prediction in Cross-Plied Curved Composite Laminates." Composite materials: Fatigue and fracture. 4 (1993): 105-26. McClennan, Scott Andrew. "Crack Growth and Damage Modeling of Fibre Reinforced Polymer Composites." MAc Thesis, Departments of Materials Engineering, the University of British Columbia, Vancouver, BC, Canada, (2004) Reddy, JN. "An Evaluation of Equivalent-Single-Layer and Layerwise Theories of Composite Laminates." Composite Structures 25.1 (1993): 21-35. 108 Reddy, JN. "A General Non-Linear Third-Order Theory of Plates with Moderate Thickness." International Journal of Non-Linear Mechanics 25.6 (1990): 677-86. Reddy, Junuthula Narasimha. Mechanics of Laminated Composite Plates and Shells: Theory and Analysis. CRC, 2003. Reddy, Junuthula Narasimha. Mechanics of Laminated Composite Plates: Theory and Analysis. CRC press, 1997. Reddy, JN. "A Simple Higher-Order Theory for Laminated Composite Plates." Journal of applied mechanics 51.4 (1984): 745-52. Reissner, E., and Yl Stavsky. "Bending and Stretching of Certain Types of Heterogeneous Aeolotropic Elastic Plates." Journal of Applied Mechanics 28.3 (1961): 402-8. Robbins, DH, and JN Reddy. "Modelling of Thick Composites using a Layerwise Laminate Theory." International Journal for Numerical Methods in Engineering 36.4 (2005): 655-77. Schellekens, JCJ, and R. De Borst. "A Non-Linear Finite Element Approach for the Analysis of Mode-I Free Edge Delamination in Composites." International Journal of Solids and Structures 30.9 (1993): 1239-53. White, Scott R., and Yeong K. Kim. "Process-Induced Residual Stress Analysis of AS4/3501-6 Composite Material." Mechanics of Composite Materials and Structures an International Journal 5.2 (1998): 153-86. White, SR, and HT Hahn. "Process Modeling of Composite Materials: Residual Stress Development during Cure. Part II. Experimental Validation." Journal of Composite Materials 26.16 (1992): 2423-53. Whitney, James Martin, and Arthur W. Leissa. "Analysis of Heterogeneous Anisotropic Plates." Journal of Applied Mechanics 36.2 (1969): 261-6. Whitney, JM, and NJ Pagano. "Shear Deformation in Heterogeneous Anisotropic Plates." Journal of Applied Mechanics 37.4 (1970): 1031-6. 109 Williams, Kevin V., Reza Vaziri, and Anoush Poursartip. "A Physically Based Continuum Damage Mechanics Model for Thin Laminated Composite Structures." International Journal of Solids and Structures 40.9 (2003): 2267-300. Wimmer, G., et al. "Computational and Experimental Investigation of Delamination in L-Shaped Laminated Composite Components." Engineering Fracture Mechanics 76.18 (2009): 2810-20. Wisnom, MR, et al. "Mechanisms Generating Residual Stresses and Distortion during Manufacture of polymer–matrix Composite Structures." Composites Part A: Applied Science and Manufacturing 37.4 (2006): 522-9. Yan, YJ, et al. "Development in Vibration-Based Structural Damage Detection Technique." Mechanical Systems and Signal Processing 21.5 (2007): 2198-211. Zhu, Qi, et al. "Three-Dimensional Viscoelastic Simulation of Woven Composite Substrates for Multilayer Circuit Boards." Composites Science and Technology 63.13 (2003): 1971-83. Zobeiry, Navid. "Extracting the Strain-Softening Response of Composites using Full-Field Displacement Measurement." PhD Thesis, Department of Civil Engineering, the University of British Columbia, Vancouver, BC, Canada. (2010)
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Integration of process simulation with damage modelling...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Integration of process simulation with damage modelling of composite laminates using layerwise elements Bakhtiarizadeh, Hamidreza 2015
pdf
Page Metadata
Item Metadata
Title | Integration of process simulation with damage modelling of composite laminates using layerwise elements |
Creator |
Bakhtiarizadeh, Hamidreza |
Publisher | University of British Columbia |
Date Issued | 2015 |
Description | During the last decades, several simulation models have been introduced in numerical modelling of composite structures. Research works in this area are focused on either manufacturing process or in-service behaviour analysis of composite materials. Although, these two phenomena are not separable and the effects of curing process on in-service behaviour of composite material is inevitable in practice, process-induced residual stresses are rarely considered in the in-service analysis of materials. The main objective of this research is to combine the process simulation and in-service analysis of composite materials within a common finite element framework. This would make it possible to model a complete life cycle of composite material from the beginning of the curing process all the way to its failure. In this approach the process-induced residual stresses are carried over to the in-service analysis and failure of the material. The model is implemented in an existing framework for process modelling developed at UBC Composites Group. Minimizing computational cost and accuracy are two important objectives in modelling of composite structures and usually there is a trade-off between these two objectives. In this research with the objective of minimizing computational cost and having the capability of capturing through-thickness stresses accurately, the Layerwise element is selected for the finite element modelling framework. The performance of the integrated process/damage simulation framework is tested through modelling of flat curved composite laminates that undergo various processing (cure) cycles and are subsequently subjected to mechanical loads. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2015-08-28 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivs 2.5 Canada |
DOI | 10.14288/1.0165785 |
URI | http://hdl.handle.net/2429/54689 |
Degree |
Master of Applied Science - MASc |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2015-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/2.5/ca/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2015_november_bakhtiarizadeh_hamidreza.pdf [ 3.16MB ]
- Metadata
- JSON: 24-1.0165785.json
- JSON-LD: 24-1.0165785-ld.json
- RDF/XML (Pretty): 24-1.0165785-rdf.xml
- RDF/JSON: 24-1.0165785-rdf.json
- Turtle: 24-1.0165785-turtle.txt
- N-Triples: 24-1.0165785-rdf-ntriples.txt
- Original Record: 24-1.0165785-source.json
- Full Text
- 24-1.0165785-fulltext.txt
- Citation
- 24-1.0165785.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0165785/manifest