awhiw v . . . \ 30k... 2. ... c in _.,_ . .2. 1%.? 4113. . Viudfl; .. an; L .1 Jrnmh‘cflz h . :r....;.l. ; .1. 4a 633. This is to certify that the i i dissertation entitled I A PENALTY-BASED INTERFACE TECHNOLOGY FOR ' CONNECTING INDEPENDENTLY MODELED I SUBSTRUCTURES AND FOR SIMULATING GROWTH OF DELAMINATION IN COMPOSITE STRUCTURES presented by ‘ my fl Antonio Pantano has been accepted towards fulfillment of the requirements for PhD . Mechanics > degree in [ WKW ‘ Major professor Date TUNE“ [@2002 MS U is an Affirmaliw Action/Equal Opportunity Institution 0- 12771 " LIBRARY Michigan State University PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE MAR 0 @004 pair) a gag: NOV 13mii°i2 8 209$ 6/01 cJCIRC/DateouepGS-p. 1 5 Sl' A PENALTY-BASED INTERFACE TECHNOLOGY FOR CONNECTING INDEPENDENTLY MODELED SUBSTRUCTURES AND FOR SIMULATING GROWTH OF DELAMINATION IN COMPOSITE STRUCTURES By Antonio Pantano A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mechanical Engineering 2002 ABSTRACT A PENALTY—BASED INTERFACE TECHNOLOGY FOR CONNECTING INDEPENDENTLY MODELED SUBSTRUCTURES AND FOR SIMULATING GROWTH OF DELAMINATION IN COMPOSITE STRUCTURES By Antonio Pantano An effective and robust interface element technology able to connect independently modeled finite element subdomains is presented. This method has been developed using penalty constraints and allows coupling of finite element models whose nodes do not coincide along their common interface. Additionally, the present formulation leads to a computational approach that is very efficient and completely compatible with existing commercial software. A significant effort has been directed toward identifying those model characteristics (element geometric properties, material properties and loads) that most strongly affect the required penalty parameter, and subsequently to developing simple “formulae” for automatically calculating the proper penalty parameter for each interface constraint. This task is especially critical in composite materials and structures, where adjacent sub-regions may be composed of significantly different materials or laminates. The present interface element has been implemented in the commercial finite element code ABAQUS as a User Element Subroutine (UEL), making it easy to test the approach for a wide range of problems. .Qx‘h huh I ”3"" A54. 1 .gn I‘ll II: (If? Once the reliability and the effectiveness of the interface element were established, new capabilities were implemented in the F E code in order to simulate delamination growth in composite laminates. Thanks to its special features, the interface element approach has several advantages over the conventional FE ones. An analysis of the literature on delamination techniques shows that generally delamination growth is simulated in a discretized form by releasing nodes of the FEM. The presented interface element allows this limitation to be overcome. It is possible to release portions of the interface surface whose length is smaller than that of the finite elements. In addition, for each portion the value of the penalty parameter can be changed at will, allowing the damage model to be applied to a desired fraction of the interface between the two meshes. The approach implemented in the interface element is one of the most commonly adopted. It mixes features of the strength of materials approach and fracture mechanics. The strength-based part predicts the onset of damage when the interfacial stress reaches the interlaminar tensile strength. Meanwhile, the relation between the work of separation and the critical value of the energy release rate provides a reliable softening model, required to describe how the stiffness is gradually reduced to zero after the onset of damage. Results for double cantilever beam DCB, end-loaded split (ELS) and Fixed- Ratio Mixed Mode (FRMM) specimens are presented. These results are compared to measured data to assess the ability of the present damage model to simulate delamination growth for Mode I, Mode II and Mixed Mode I+II. Copyfightby Antonio Pantano 2002 To my Mother and Father lIIOPC II 1 ‘: AWL.) 3d\hxfi acadsn‘ 'V, In mam... ACKNOWLEDGMENTS I will always be grateful to my beloved parents for years of unconditional support. I hope they will be proud of me, now and in the future. I would like to reserve a special gratitude to my sweet wife Giulia. Without her love this result would have been harder to achieve. I definitely recognize the immense fortune I had to meet a person like my research advisor Prof. Ronald Averill. His uncommon humanity equals only his exceptional academic knowledge. I also greatly appreciate Prof. Gary J. Burgess, Dahsin Liu and Hungyu Tsai for making efforts to serve my thesis committee. This work was sponsored by NASA Langley Research Center under grant NAG- 1-2213. The authors are grateful for the helpful comments of Dr. Jonathan Ransom, the grant technical monitor. Partial support was also provided by the State of Michigan Research Excellence Fund. Antonio Pantano June 2002 East Lansing, Michigan vi IBIO LBIC CHAP CHA CH TABLE OF CONTENTS LIST OF TABLES ....................................................................................... XII LIST OF FIGURES ................................................................................... XIII CHAPTER 1 INTRODUCTION TO THE INTERFACE ELEMENT ....... 1 1.1 Introduction and Literature Review ....................................... l 1.2 Organization of the Thesis ..................................................... 5 CHAPTER 2 GENERAL DESCRIPTION OF THE INTERFACE ELEMENT ............................................................................. 6 2.1 Hybrid Interface Method ........................................................ 7 2.2 Penalty Hybrid Interface Method ........................................... 9 2.3 Determination of the Penalty Parameters ............................ 1 1 2.4 Automatic Round-Off Error Control ................................... 14 2.5 Implementation as Abaqus User Element Subroutine ......... 17 CHAPTER 3 TIMOSHENKO BEAM ....................................................... 18 3.1 Extensional Load — One Element ........................................ 18 3.1.1 Lagrange Multiplier Method .................................................. 19 3.1.2 Penalty Method ...................................................................... 21 3.1.3 Comparison between the Two Methods ................................ 22 3.1.4 Relation between Penalty Parameter and Beam Properties... 23 3.2 Extensional Load — Two Elements ...................................... 24 3.2.1 Lagrange Multiplier Method .................................................. 25 3.2.2 Penalty Method ...................................................................... 26 3.2.3 Comparison between the Two Methods ................................ 28 vii CHAD 3.3 3.4 3.5 3.6 3.7 3.8 CHAPTER 4 4.1 4.2 4.3 3.2.4 Relation between Penalty Parameter and Beam Properties 29 Bending Loads — One Element — Full Integration ............... 30 3.3.1 Lagrange Multipliers Method ................................................ 31 3.3.2 Penalty Method ...................................................................... 35 3.3.3 Comparison between the Two Methods ................................ 38 3.3.4 Relation between Penalty Parameter and Beam Properties 39 Bending Loads — One Element — Reduced Integration ....... 43 3.4.1 Lagrange Multipliers Method ................................................ 43 3.4.2 Penalty Method ...................................................................... 46 3.4.3 Comparison between the Two Methods ................................ 49 3.4.4 Relation between Penalty Parameter and Beam Properties 51 Bending Loads — Two Element — Reduced integration ....... 54 3.5.1 Lagrange Multiplier Method .................................................. 55 3.5.2 Penalty Method ...................................................................... 60 Full versus Reduced Integration .......................................... 64 Interface Element ................................................................. 67 3.7.1 Stiffness Matrix ..................................................................... 67 3.7.2 Best Values of the Penalty Parameters .................................. 69 Numerical Results ................................................................ 74 3.8.1 One Element — Clamped Beam .............................................. 74 3.8.2 Two Elements — Clamped Beam ........................................... 75 3.8.3 Two Elements — Two Materials — Clamped Beam ................ 76 3.8.4 Three Elements — Simply Supported Beam ........................... 77 PLANE STRESS QUADRILATERAL ELEMENT ........... 79 Finite Element Model ........................................................... 79 Cubic Spline Interpolation Functions .................................. 82 4.2.1 General Form of the Cubic Spline ......................................... 83 4.2.2 Natural Cubic Spline over Three Points ................................ 87 2D Beam with a Vertical Interface — Single Variable ......... 9O viii C HA 4.4 4.5 4.6 4.7 CHAPTER 5 5. 1 5.2 5.3 5.4 5.5 Clamped Rectangular Plane Stress Element ........................ 98 4.4.1 Axial Load ........................................................................... 100 4.4.2 Transversal Load ................................................................. 1 01 Penalty Parameter and Element Properties ........................ 103 4.5.1 Axial Extension ................................................................... 103 4.5.2 Vertical Flexure ................................................................... 104 Building an interface Element Using Abaqus UEL ........... 105 4.6.1 Automatic Choice of the Optimal Penalty Parameter .......... 106 Numerical Results .............................................................. 108 4.7.1 Patch Test ............................................................................. 108 4.7.2 Isotropic 2D Cantilever Beam with a Vertical Interface ..... 112 4.7.2.1 Axial Load ..................................................................... l 13 4.7.2.2 Bending Load ................................................................ 1 14 4.7.3 Tension-Loaded Plate with a Central Circular Hole ............ 1 16 4.7.4 Composite 2D Cantilever Beam with a Vertical Interface .. 128 4.7.5 Clamped-Clamped Asymmetric Beam ................................ 132 PLATES ............................................................................. 135 Finite Element Model ......................................................... 135 Clamped Rectangular Plate Element ................................. 140 5.2.1 Transversal Load ................................................................. 1 41 5.2.2 Bending Moment ................................................................. 142 Penalty Parameter and Element Properties ........................ 143 5.3. 1 Transversal Load ................................................................. 143 5.3.2 Bending Moment ................................................................. 145 Building an interface Element for Plates ........................... 146 5.4.1 Automatic Choice of an Optimal Penalty Parameter ........... 147 Numerical Results .............................................................. 150 5.5.1 Plate Bending ....................................................................... 151 5.5.2 Simply Supported Plate under Sinusoidal Load .................. 152 ix CHAPIEI CHAPIEI CHAPTER 6 6.1 6.2 6.3 6.4 CHAPTER 7 7.1 7.2 7.3 3D LINEAR BRICK ELEMENT ...................................... 157 Automatic Choice of the Penalty Parameter ...................... 157 Cubic Spline Interpolation Functions ................................ 159 Building an Interface Element for 8-Node Brick ............... 160 Numerical Results .............................................................. 160 6.4.1 3D Cantilever Beam with 3 Square Section ........................ 160 6.4.1.1 Axial Load ..................................................................... 161 6.4.1.2 Bending Load ................................................................ 164 6.4.2 3D Cantilever Beam with a Quadrilateral Section .............. 168 6.4.3 3D Shafts Assembly ............................................................ 176 6.4.3.1 Bending Load ................................................................ 177 6.4.3.2 Torsion .......................................................................... 185 INTRODUCTION TO THE MODELING OF DELAMINATION GROWTH .......................................... 193 Causes of Delamination ..................................................... 193 7.1.1 Free-edge stresses ................................................................ 193 7.1.2 Matrix cracks ....................................................................... 194 7.1.3 Impact .................................................................................. 196 7.1.4 Residual thermal stresses and moisture ............................... 197 Experimental Studies on Delamination Growth ................ 197 7.2.1 Definition of the strain energy release rate G ...................... 197 7.2.2 Test specimens for evaluating G .......................................... 200 7.2.2.1 Test specimens for Mode l ............................................ 200 7.2.2.2 Test specimen for Mode 11 ............................................ 202 7.2.2.3 Test specimen for Mixed Mode l+lI ............................. 202 Delamination Criteria ......................................................... 203 7.3.1 Fracture Mechanics Based ................................................... 204 7.3.1 .1 Numerical evaluation of the strain energy release rate. 204 7.3.1.2 Fracture mechanics based delamination criteria ........... 21 1 7.3.2 Interlaminar strength and mixed approaches ....................... 213 CHAPIEF CHAPI APPEX REFER CHAPTER 8 A NOVEL INTERFACE TECHNOLOGY FOR MODELING DELAMINATION ...................................... 2 19 8.1 Damage Model ................................................................... 220 8.2 Mixed Mode Approach ...................................................... 224 8.3 Energetic Approach for a Correct Distribution of Stress... 231 8.4 Friction Model .................................................................... 237 8.4.1 Evaluation of the force at the interface ................................ 237 8.4.2 Implementation of a Friction Model .................................... 239 8.4.3 Numerical Test for the Friction Model ................................ 240 8.5 Numerical Results .............................................................. 242 8.5.1 Double Cantilever Beam Test #1 ......................................... 242 8.5.2 Double Cantilever Beam Test #2 ......................................... 246 8.5.3 End-Loaded Split (ELS) Beam Test .................................... 247 8.5.4 End-Loaded Split (ELS) Beam Test - With Friction .......... 249 8.5.5 Fixed-Ratio Mixed Mode (FRMM) Test ............................. 251 CHAPTER 9 CONCLUSIONS ................................................................ 254 APPENDIX A USER MANUAL ............................................................... 257 REFERENCES ........................................................................................... 268 xi Table 1. Table 2. Table 3. Table 4. Table 5. Table 6. Table 7. Table 8. Table 9. Table 10. Table 11. Table 12. Table 13. LIST OF TABLES Penalty parameters associated with each of the constraints enforced ........... 68 Penalty parameters relations for a Timoshenko beam element. .................... 73 Single beam clamped - Outcomes ................................................................. 75 Two beams clamped - Outcomes ................................................................... 76 Two beams — two materials - clamped - Outcomes ....................................... 77 Three elements simply supported beam - Outcomes ..................................... 78 Tip axial displacement numerical results for beam extension. .................... l 13 2D Cantilever Beam under axial load - Mesh 10x8-10x4 .......................... 1 14 Tip deflection numerical results for beam flexure ....................................... 1 15 Tip deflection in a 2D beam under bending load — mesh 10x8-10x4 .......... 1 15 Expression for computing the penalty parameters in a plate element. ........ 150 Nondimensionalized results for the center deflection of the plate .............. 152 Tip deflection numerical results for beam flexure ....................................... 164 xii Figure 1. Figure 2. Figure 3. Figure 4. Figure 5. Figure 6. Figure 7. Figure 8. Figure 9. Figure 10. Figure 11. Figure 12. Figure 13. Figure 14. Figure 15. Figure 16. Figure 17. Figure 18. Figure 19. LIST OF FIGURES 2D (a) and 3D (b) Interface element configurations ........................................ 7 Beam under an axial load applied at the tip — One Element .......................... 19 Beam under an axial load applied at the tip — Two Elements ....................... 24 Beam under bending loads — One Element .................................................... 31 Beam under bending loads —- Two Elements ................................................. 55 Interface element for the Timoshenko beam element .................................... 68 Beam under axial loads — Two Elements ...................................................... 69 Beam under bending loads — One Element .................................................... 74 Two beam elements joined through an interface element ............................. 76 Two beam elements joined through an interface element ............................. 77 Three elements simply supported beam ........................................................ 78 Geometrical dimensions of the rectangular element. .................................... 80 Representation of the cubic splines over the first interval ............................. 89 Representation of the cubic splines over the second interval ........................ 89 2D Cantilever beam under axial load applied at the tip ................................. 90 Single element beam connected to a fixed penalty frame ............................. 98 Single element beam under axial load ...................................................... 100 Single element beam under a transversal load ............................................. 101 Patch of incompatible 2D meshes connected by one interface element. ..... 109 xiii Figure 30. Figure 31. figure I: ""2 (L (I :1 L Figure 20. Figure 21. Figure 22. Figure 23. Figure 24. Figure 25. Figure 26. Figure 27. Figure 28. Figure 29. Figure 30. Figure 31. Figure 32. Figure 33. Figure 34. Figure 35. Figure 36. Figure 37. Figure 38. Patch of incompatible 2D meshes - Displacements distribution in the direction 1 (U1). ........................................................................................... 110 Patch of incompatible 2D meshes - Displacements distribution in the direction 2 (U2). ........................................................................................... 110 Patch of incompatible 2D meshes - Displacements distribution in the direction 1 (a) and in direction 2 (b) under uniform shear deformation. ..... 1 11 2D Cantilever beam with vertical interface ................................................. l 12 Stability of the solution for different values of y- axial load. ..................... 1 14 Stability of the solution for different values of y - bending load ................ 1 15 Tension-loaded plate with a central circular hole (not to scale). ................. 116 (a) Interface element and (b) conventional FE mesh. .................................. l 18 (a) Interface element and (b) conventional FE mesh - zoom of the deformed configuration. ............................................................................................... 1 19 (a) Interface element and (b) conventional FE solution - horizontal displacement distribution (U1). .................................................................... 120 (a) Interface element and (b) conventional FE solution - vertical displacement distribution (U2). .......................................................................................... 121 (a) Interface element and (b) conventional FE solution - stress distribution in the direction 1 (on). .................................................................................... 122 (a) Interface element and (b) conventional FE solution - stress distribution in the direction 2 (022). .................................................................................... 123 Stress Concentration Factor Along the axis 2 .............................................. 124 Interface element FE mesh #2 - zoom of the deformed configuration. ....... 125 (a) Interface element and (b) conventional FE solution - horizontal displacement distribution (U1). .................................................................... 126 (a) Interface element and (b) conventional FE solution - vertical displacement distribution (U2). .......................................................................................... 127 2D Composite cantilever beam with vertical interface ................................ 128 Beam extension along the thickness at the free end under axial load. ........ 130 xiv Figure 39. Figure 40. Figure 41. Figure 42. Figure 43. Figure 44. Figure 45. Figure 46. Figure 47. Figure 48. Figure 49. Figure 50. Figure 51. Figure 52. Figure 53. Figure 54. Figure 55. Figure 56. Figure 57. Figure 58. Figure 59. Figure 60. Beam vertical deflection along the thickness at the free end under axial load. ..................................................................................................................... 130 Beam extension along the thickness at the free end under transversal load. 131 Beam vertical deflection along the thickness at the free end under transversal load. ............................................................................................................. 131 Clamped-clamped two-piece asymmetric beam. The data in parentheses represents the x and y coordinates of the points indicated .......................... 132 Deformed configuration obtained through conventional Abaqus FE solution. ..................................................................................................................... 133 Deformed configuration obtained using interface elements approach. ....... 133 Axial displacements along the thickness at the interface ............................ 134 Transversal displacements along the thickness at the interface ................... 134 Geometrical dimensions of the plate element .............................................. 138 Single plate element connected to a fixed penalty frame — Top view ......... 141 Single plate element connected to a fixed penalty frame — Lateral view 141 Single plate element under transversal load applied at the free side ........... 142 Single plate element under bending moment applied at the free side ......... 143 Original and deformed configurations of the plate ...................................... 151 Penalty interface FE model of a flat plate. .................................................. 153 Original and deformed mesh as seen from the bottom of the plate. ............ 153 Distribution of traversal displacements w - Penalty interface frame FEM. 154 Distribution of traversal displacements w - Conventional FEM. ................ 154 Variation of the transverse displacement w along one of the two sides of symmetry of the model ................................................................................ 156 Three-dimensional solid 8-node linear brick element. ................................ 157 Definition of a local coordinate system ....................................................... 158 3D Cantilever beam with a square section .................................................. 161 XV Figure 61. I. Figure 03. I Iiguebi ( FigureM. ( Figurcbi. i Figuretib. Figureb’. figurcDS. Iigurer. Figure”i_i, Figure71, Figure 73, Figure 76. Figure 77. “we ’8. F1Litre 79; Figure 61. Figure 62. Figure 63. Figure 64. Figure 65. Figure 66. Figure 67. Figure 68. Figure 69. Figure 70. Figure 71. Figure 72. Figure 73. Figure 74. Figure 75. Figure 76. Figure 77. Figure 78. Figure 79. Figure 80. Figure 81. Figure 82. Figure 83. Figure 84. Cantilever beam with a square section — Axial displacement U3 — 3D view ..................................................................................................................... 162 Cantilever beam with a square section — Axial displacement U3 — Side view ..................................................................................................................... 162 Cantilever beam with a square section - Stress S33 .................................... 163 Cantilever beam with a square section — Lateral displacement U1 ............. 163 Cantilever beam with a square section -— Lateral displacement U2 ............. 164 ABAQUS reference model — Displacement U2 .......................................... 165 Interface element model — Displacement U2 — (Classical deflection = 4.0) 165 ABAQUS reference model — Displacement U3 .......................................... 166 Interface element model — Displacement U3 ............................................... 166 ABAQUS reference model - Stress S33 — (Classical = 6.0E+4) ................ 167 Interface element model — Stress S33 — (classical = 6.0E+4) ..................... 167 3D Cantilever beam with a quadrilateral section ......................................... 168 Quadrilateral section of the beam. ............................................................... 169 Original and deformed mesh ....................................................................... 170 ABAQUS reference model — Displacement U2 .......................................... 171 Interface element model — Displacement U2 ............................................... 171 ABAQUS reference model — Displacement U1 .......................................... 172 Interface element model — Displacement U1 ............................................... 172 ABAQUS reference model — Displacement U3 .......................................... 173 Interface element model — Displacement U3 ............................................... 173 ABAQUS reference model — Lateral view 1 -— Stress S33. ......................... 174 Interface element model — Lateral view 1 — Stress S33 ............................... 174 ABAQUS reference model — Lateral view 2 — Stress S33 .......................... 175 Interface element model — Lateral view 2 — Stress S33 ............................... 175 xvi Figure 35. Figure 86. Figure 87. Figure 88. Figure 89 Figure 9‘5“ Figure 91 Figure 9: Figure 93 Figure 0. Figure 9 Figure 9 Figure 9 Figure ‘5‘ Figure F Ifigure ‘ Figure Figure 85. Shaft assembly. ............................................................................................ 176 Figure 86. Original and deformed mesh under bending load ........................................ 178 Figure 87. ABAQUS reference model — Displacement U1 .......................................... 179 Figure 88. Interface element model — Displacement U1 ............................................... 179 Figure 89. ABAQUS reference model — Displacement U2 .......................................... 180 Figure 90. Interface element model — Displacement U2 ............................................... 180 Figure 91. ABAQUS reference model — Displacement U3 .......................................... 181 Figure 92. Interface element model — Displacement U3 ............................................... 181 Figure 93. ABAQUS reference model - Stress Sll ..................................................... 182 Figure 94. Interface element model — Stress Sll .......................................................... 182 Figure 95. ABAQUS reference model — Stress S22 ..................................................... 183 Figure 96. Interface element model — Stress 822 .......................................................... 183 Figure 97. ABAQUS reference model — Stress S33 ..................................................... 184 Figure 98. Interface element model — Stress S33 .......................................................... 184 Figure 99. Deformed mesh under torsion load .............................................................. 185 Figure 100. ABAQUS reference model — Displacement U1 .......................................... 186 Figure 101. Interface element model — Displacement U1 ............................................... 186 Figure 102. ABAQUS reference model — Displacement U2 .......................................... 187 Figure 103. Interface element model — Displacement U2 ............................................... 187 Figure 104. ABAQUS reference model — Displacement U3 .......................................... 188 Figure 105. Interface element model — Displacement U3 ............................................... 188 Figure 106. ABAQUS reference model — Stress 812 ..................................................... 189 Figure 107. Interface element model — Stress 812 .......................................................... 189 Figure 108. ABAQUS reference model — Stress Sll ..................................................... 190 Figure 109. Interface element model — Stress Sll .......................................................... 190 xvii Figure l 10. .- Figure ll 1. Figure Ill. Figure l 13 Figure l 14 Figure l 15 Figure llt Figure 11 Figure I: Figure 1: Figure I Figure I Figure l Figure 110. ABAQUS reference model — Stress $22 ..................................................... 191 Figure 111. Interface element model — Stress 822 .......................................................... 191 Figure 112. ABAQUS reference model — Stress S33 ..................................................... 192 Figure 113. Interface element model — Stress S33 .......................................................... 192 Figure 114. Double cantilever beam (DBC) specimen ................................................... 201 Figure 115. End-loaded split (ELS) specimen ................................................................ 202 Figure 116. Fixed-ratio mixed-mode(FRMM) specimen ................................................ 203 Figure 117. 2D FEM before crack opening (a) and after crack opening (b) ................... 208 Figure 118. Interfacial constitutive model ...................................................................... 214 Figure 119. Bilinear interfacial constitutive damage model. .......................................... 220 Figure 120. Division of the interface element in intervals .............................................. 223 Figure 121. Graph fix/5,20 versus 82/520. ........................................................................... 226 Figure 122. Updated interfacial constitutive model for mode I delamination ................ 227 Figure 123. Updated interfacial constitutive model for mode II delamination ............... 227 Figure 124. Graph Gn/G”c versus Gl/G.c ......................................................................... 229 Figure 125. Final interfacial constitutive model for mode I delamination ..................... 230 Figure 126. Final interfacial constitutive modeI for mode II delamination .................... 231 Figure 127. Fracture mechanics model. .......................................................................... 232 Figure 128. Normal stress distribution along the entire interface. .................................. 232 Figure 129. Boundary conditions and the geometry of the test for the friction model 240 Figure 130. Test Friction — Graph Force-Displacements. ............................................... 241 Figure 131. Geometry and boundary conditions for the DCB test specimen. ................ 242 Figure 132. Original and deformed (5mm opening) models of the DCB test specimen using a 120x2 mesh. .................................................................................... 243 Figure 133. Force vs. displacement results of experimental and numerical DCB test 244 xviii Figure 1." Figure In: Figure 13! Figure 13' Figure 14. Figure 14. Figure 14- Figure 134. Force vs. displacement results of DCB test. Convergence of the solution with intervals number. Mesh 300x8. ................................................................... 245 Figure 135. Convergence of the solution with the number of increments. Mesh 300x8. 245 Figure 136. Geometry and boundary conditions for the DCB test #2 ............................ 246 Figure 137. Force versus displacement results of experimental and numerical DCB test #2 ..................................................................................................................... 247 Figure 138. Geometry and boundary conditions for the ELS test specimen. ................. 248 Figure 139. Original and deformed (30mm tip displacement) models of the ELS test specimen for a 300x8 mesh. ........................................................................ 248 Figure 140. Force vs. displacement results of experimental and analytical ELS test for varying number of loading increments. ....................................................... 249 Figure 14l.Original and deformed (30mm tip displacement) models of the ELS test specimen for a 300x8 mesh. ........................................................................ 250 Figure 142. Force vs. displacement results of experimental and numerical ELS test in presence of friction. ..................................................................................... 251 Figure 143. Geometry and boundary conditions for the F RMM test .............................. 252 Figure 144. Original and deformed (20mm tip displacement in absence of delamination) models of the FRMM test specimen. ........................................................... 252 Figure 145. Force vs. displacement results of experimental and numerical F RMM test.253 xix 1.1 l CHAPTER 1 INTRODUCTION TO THE INTERFACE ELEMENT 1.1 Introduction and Literature Review The ways in which analysis and design are performed have changed extensively during the past decade. Automated design algorithms are now plentiful and increasingly robust, and no longer are individual components of a structure designed in a vacuum. Facilitated in part by product data management (PDM) and product lifecycle management (PLM) software and intemet-based data sharing, integrated design activities are now possible. Further, the speed and economy of modern computers have enabled engineers to perform many large scale analyses that were unheard of a decade ago. It is now common to analyze and simulate the response of entire aircraft, spacecraft, automobiles, ships, and other structural assemblies to a variety of complex and combined types of loading. The CPU time required to perform such analyses is very small compared to the time required for engineers to create the mathematical and computer models, and the latter effort is usually the most expensive component of a large scale analysis or simulation. With model sharing and large scale analysis activities on the rise, it is becoming 1110-; ll 16 he i i 12‘ Ir" '. ' 117.0 h evident that improved technology for building computer models is needed. One issue that arises often is the need to perform a unified analysis of a structural assembly using sub- structural models created independently. These sub-structural models are frequently created by different engineers using different software and in different geographical locations, with little or no communication between the teams of engineers creating the models. As a result, these models are likely to be incompatible at their interfaces, making it very difficult to combine them for a unified analysis of the entire assembly. Finite element interface technology has been developed during the past decade to facilitate the joining of independently modeled substructures. Unconventional approaches have been employed to connect special elements based on analytical solutions to finite element models [1-2]. In [1] a near-field solution for a dynamically propagating crack at the interface of two dissimilar anisotropic elastic materials was successfully implemented into a hybrid-displacement finite element formulation. In [2] Lagrange multiplier terms are used to couple a special element, based on a stress analytical solution, to standard finite elements. In order to take advantage of parallel computing. Farhat and colleagues [3-4] developed a domain decomposition approach for partitioning the spatial domain into a set of disconnected subdomains, each assigned to an individual processor. Lagrange multipliers were introduced to enforce compatibility at the interface nodes. In [5] non- conforming “mortar” elements are employed to connect incompatible subdomains using a conjugate gradient iterative technique in a domain decomposition scheme designed for parallel computers. The finite element interface technology developed at NASA LaRC [6-10] and V HIST. d Cur“. 2*“ U .5! .l {JV-1.," elsewhere [11] allows the connection of independently modeled substructures with incompatible discretization along the common boundary. This approach has matured to a point that it is now very effective. However, because the interface technology utilizes Lagrange multipliers to enforce the interface constraint conditions, the resulting system of equations is not positive-definite. Hence, state-of-the-art sparse solver technology cannot be utilized with the interface technology. “Recently, an alternative formulation for the finite element interface technology based on Lagrange multipliers has been developed [12]. The alternative approach recasts the interface element constraint equations in the form of multi-point constraints. This change allows an easier implementation of the formulation in a standard finite element code and alleviates the issues related to the resulting non-positive definite system of equations. The method seems to provide reliable results, but the formulation of the interface method is still quite complicated. A possible remedy for these shortcomings is to modify the current hybrid variational formulation of the interface element by enforcing the interface constraints via a penalty method as opposed to the current Lagrange multiplier approach. The primary consequences of this modification will be (i) a simple formulation that is easily implemented in commercial finite elements codes, (ii) a positive-definite and banded stiffness matrix and (iii) a reduced number of DOFs, since the Lagrange multiplier degrees of freedom (DOFs) will not be present. Thus, the proposed approach should greatly enhance the computational efficiency of the interface element technology. From an accuracy point of view, the penalty method enforces the constraints only appro'ximately, depending on the value of the penalty parameter chosen, while the C 10C im e PFC) pa. :3‘ (1! .,. rf r k 'E Lagrange multiplier approach enforces the constraints exactly. The penalty method interface approach was recently attempted using a single global value of the penalty parameter to enforce all constraints [13]. This study demonstrated the validity and the effectiveness of the penalty approach in an interface element. However, there is need for specific guidelines regarding the selection of an appropriate value of the penalty parameter, especially when the substructures to be connected have different material and/or section stiffnesses. There is a large body of literature related to the determination of an optimal value of the penalty parameter (see, e.g. [14-26]). However, there is no existing criterion for choosing the penalty parameter in the framework of the interface element under investigation. The determination of such a criterion is one of the goals of the effort presented herein. An effective and robust interface element technology has been developed using the penalty method. A significant part of the work has been directed toward identifying those model characteristics that most strongly affect the required penalty parameter, and subsequently to developing simple “formulae” for automatically calculating the proper penalty parameter for each interface constraint. The new approach has been validated through a wide variety of one, two-dimensional and three-dimensional problems that have been investigated extensively from both an analytical and a computational point of view. Finally, the penalty based interface element technology has been implemented into an existing commercial code. [.2 0’! A 1 Applicant Chapter 3 dCWIUpCi meshes 0 Ch. for mode capabiii‘. composi Preseme 1.2 Organization of the Thesis A general description of the interface element is presented in Chapter 2. Applications of the hybrid and penalty interface methods to ID beams are provided in Chapter 3. A version of the interface element for 2D plane elasticity and plate elements is developed in Chapter 4 and Chapter 5, respectively. The interface element for connecting meshes of 3D solid 8-node brick elements is described in Chapter 6. Chapter 7 introduces to the development of the novel interface element technology for modeling crack growth, by reviewing the literature on the subject. The additional capabilities implemented in the interface element in order to simulate delamination in composite laminates and the numerical tests performed to validate the approach are presented in Chapter 8. it'd): CORE: inter} then and ifidi ' ,_ O O- . 4 CHAPTER 2 GENERAL DESCRIPTION OF THE INTERFACE ELEMENT Consider two independently modeled subdomains Q, and Q, as shown in Figure 1(a) and 1(b), respectively, for a 2D and for a 3D geometry. The two substructures are connected to each other using an interface element acting like “glue” at the common interface. The interface element is discretized with a set of nodes that are independent of the nodes at the interface in subdomains Q, and Q, . Both in the hybrid interface method and in the penalty method, the coupling terms associated to the interface element are arranged in the form of a “stiffness” matrix and assembled with the other finite element stiffness matrices as usual. n. q. Q; i» . O . S O . V Flgur 2.1 r multip} le‘l» FAIL"; (b) Figure 1. 2D (a) and 3D (b) Interface element configurations. 2. 1 Hybrid Interface Method The hybrid interface method [6-12] introduces two vectors of Lagrange multipliers ’11 and ,1, in the total potential energy (TPE) of the system to satisfy the displacement continuity conditions. Thus the TPE of the system assumes the form: 7t=7rnl+7rnz+IA(V—u,)ds+I/i,,(V—u2)ds (1) S S The nodal displacements of the sub-domain Qj are identified by qj and qj. The 7 superscript 0 identifies the degrees of freedom (DOFs) that are not on the interfaces, while i denotes DOFs that are on the interfaces. The displacement field u j of the sub- domain O1 is expressed in terms of the unknown nodal displacements q; as: u, = N j 61’. (2) where N, can be the matrices of linear Lagrange interpolation functions. The displacement field V is approximated on the entire interface surface in terms of unknown nodal displacements q, as: V=Ta (D where T is a matrix of cubic spline interpolation functions [27-38]. The first variation of 7: is taken with respect to all the DOFs and the vectors of Lagrange multipliers xii and 22: 5 ' r = 0 4 ”IQi)~CIis‘ISsq§»‘12941»4’2 H On the interface part of the subdomains the following equations result: u_l=V, xil=tj, A1+xiQ=0 forj=l,2 (5) where t, is the traction on the interface. These equations show that: 0 Displacement continuity is enforced, o itj are interface tractions, o the total traction is zero on the interface. The first variation of ayields the system of equations: 'K," K,“ 0 0 0 M, 0 F ’q,’ " if," K," K,"” 0 0 o 0 0 q;’ f,” 0 0 K 3' K3” 0 0 M: q; f! 0 0 K," K,” 0 0 0 C 9 the parameter ,6 is reduced according to: ,B""“' = g 3 "U W The interface element stiffness matrix is recalculated using the new value fl = ,8 . This approach reduces the risk that round-off errors could adversely affect the solution. Thus, the initial value of ,6 can be increased, in order to get a higher degree of 16 3d 1 b} that eler The Zillion“ accuracy, knowing that it will be automatically reduced if rounding errors don’t allow that precision to be realized. 2.5 Implementation as Abaqus User Element Subroutine To test the behavior of the penalty method based “interface” element and the accuracy of the obtained relations for an automatic choice of the penalty parameters, the element has been implemented in the commercial finite element code ABAQUS as a User Element Subroutine (UEL) [39]. The UELsubroutine receives all the necessary information about geometry and material properties of the two connected meshes of finite elements from the input file. Then, the stiffness matrix of the interface element is built as previously defined: —! p GIN _Gll.\ 0 L 0 —G," G; The appropriate values of penalty parameters for each constraint are assigned automatically inside the subroutine. 17 H} .\le PW FM 1 \ l CUFF. ‘A CHAPTER 3 TIMOSHENKO BEAM In this chapter the Hybrid Interface Method and the Penalty Frame Method are applied to connect Timoshenko beam elements. The first goal is to define the stiffness matrix that could be associated with an interface element able to accomplish the link. The Hybrid Interface Method defines the upper limit to the accuracy of the Penalty Frame Method and lets us explore the relation between the Lagrange multiplier and the penalty parameter for the Timoshenko beam elements. Very simple problems are studied in order to pursue our objectives. FE analyses involving axial and bending loads are considered separately. Problems related to flexural DOFs are modeled both with conventionally formulated FEM of the Timoshenko beam and with a Reduced Integrated Element. 3. I Extensional Load — One Element The first problem to be studied is one of the simplest available: the extension of a beam under a uniform load P applied at the tip. This case is analyzed using a single linear beam element connected to one fixed point V by a displacement continuity constraint, which is imposed through a Lagrange multiplier or a Penalty parameter. The configuration of the geometry and the mesh is plotted in Figure 2. 18 Fio 3.1 9‘01: Thu L V 1 E, A 2 O t a v ll. [12 Figure 2. Beam under an axial load applied at the tip — One Element 3.1.1 Lagrange Multiplier Method The hybrid interface method introduces a Lagrange multiplier A, in the total potential energy (TPE) of the system to satisfy the displacement continuity condition. Thus the TPE of the system in this study assumes the following form. 2r=7r.+A-(v-u.)—r.-v (28) where 7:, is the TPE of the bar; r, and v are, respectively, the reaction force and the displacement at point V. Expanding and simplifying: 7r=lI-;—(o'x-8,)dV—;fi-u,+fl,-(v—u,)—r,-v (29) 7r=%-JEA(%) dx—ifi-u,+l,-(v—u,)—r,-v (30) where E is the Young’s modulus. A is the cross sectional area and fl are the applied nodal forces. If the displacements are approximated linearly by: 19 Setting to it Fri i ' it mug e. Thus. the re Since ti, uziuij =u,(l—%)+u,£—xL-) (31) _/‘=1 7: takes the following form. EA " dN dN 7r=-2—-0[l ”El[;u’— dx ——]dx— P u,+/i, (v— u,)— r,~v (32) Setting to zero the first variation of 7: , taken with respect all the DOFs, produces the following equations. an N, _= _. 4:0 33 6,, “iiwiiufi: ——’]dx- < ) " dN gleA- [[dN2][ZuJ——’—]dx—P=O (34) auz 0 dx .1 dx an 3:41—VFO (35) 275:,,_u,:0 (36) 6A Thus, the resulting FE model is: F EA EA 1 — —— 0 —1 a L L Fur F0 EA EA u P —— —- 0 0 l 2 = L L vi W (37) O 0 0 1 M) ,0 _ —l 0 l 0_ r. = A = -P (38) u, = v = 0 (39) PL “2 — E (40) I116 with 006. 3.1.2 Pr permit} p: APpmim; lite tirst t2 Pi’lirnercr, The ”suit The solution provided by the hybrid interface method coincides with the exact theoretical one. 3.1.2 Penalty Method In the penalty method the displacement continuity constraint is imposed through a penalty parameter 7, . Therefore the TPE of the system takes the form: 7r=7r,+%7,-(v—u,)2—r,-v (41) l 2 l 2 7r=IJ—2-(oi,'£,)dV—;f,-u,+57,~(v—u,) —r,-v (42) a—l-IIEA(fl)2dx—if-u +l -(v—u)2—r-v (43) 2 0 dx 1:!“ r 271 l r Approximating the displacements linearly: EA "[ 7r=——- dN dN, l 2 2 Eur-XMEujgx—de-P-ufigrr'(v-u.) -r.-v (44) 0 ’ l The first variation of 7r in this case is taken with respect all the DOFs, but not the penalty parameter. §E=EA-:J(%)[};ujddijjdx-n‘(v—u.)=0 (45) SizEA-§[dgz)(;uj%jclx-P=O (46) a3:)r—=7,c(v—u)—r,=0 (47) The resulting FE model is: 21 50h 3.1.3 Milieu rEA + y EA 7 ‘ _ l _— — l L EA EA u' 0 — —— —— 0 u2 = P (48) L L v r ’71 0 71 I Solving the system of equations: It = -P (49) u. = f— (50) 7t u.=[i+_L_]p (51) “ 7, EA 3.1.3 Comparison between the Two Methods When the penalty parameter approaches infinity, the penalty method solution tends to the one obtained using the Lagrange multiplier method. Limit u, = Limit—I: = 0 (52) rr—m n—m 7, Limit u, = Limit -1— + i P = (LJP (53) 71"”) ‘- 71—>00 7, EA EA Moreover, the constraint that has been enforced is: C, (u,.u,.v)=(v—u,) (54) If we define u,,,u,7,vr as the solutions derived from the penalty formulation, it is possible to verify the well-known relation between the Lagrange multiplier and the penalty parameter. 22 2, :7, -C, (u,,,u,,.v,) :7, -(v, —u,,)= 7, [om-1:] = —P (55) 3.1.4 Relation between Penalty Parameter and Beam Properties In this section the obtained solutions are investigated in order to find a relation between material and geometrical properties of the beam and the penalty parameter. As underlined previously, the exact displacement of the tip of the beam matches that from the Lagrange multiplier finite element formulation. .W L u, = — P 56 - EA ( ) The Penalty parameter solution u§’"”"”"' differs from the exact one by the presence of an additional term (P/ 7,): ufcnul’y : [i + L] P (57) 7, EA The ratio uf"”"”"’/ ug'w” is evaluated in order to identify the relationship between the two solutions: 1 L penalty ~—-—— + _ P u2 _ 7, EA (’5?) eruc't : 1 + (58) “2' F 1;] P 71 I EA The penalty parameter 7, is now substituted by: EA 7: = ,3 ' (T) (59) It follows that the ratio between the solutions becomes independent of material and geometrical properties of the beam. 23 ICU 111C 3.2 C00 Figun penalty 1 =1+— 60 fl ( ) As underlined in the previous chapter, the accuracy of the solution comes to depend directly on the value assigned to the parameter ,6. For example if ,8 is equal to 1000 the Penalty solution will differ from that obtained by use of the Lagrange multipliers by 0.1%. The degree of precision of the solution cannot be indefinitely increased, since round off amplification error would rise. 3.2 Extensional Load — Two Elements The previous analysis has been modified in order to analyze the imposition of a continuity constraint between two beams elements. Two linear beams elements are connected using a Lagrange multiplier or a Penalty parameter continuity constraint. The configuration of the geometry and the mesh is depicted in Figure 3. 1 E, An, L 2 V 3 E2, A2, L 4 ll] “2 v “3 “4 Figure 3. Beam under an axial load applied at the tip —- Two Elements 24 3.2.1 The 101 Where , ll 1 'F—._. n- , l .r:-- ‘1. I. l Dept-nee SUhstitut 5611an [h 3.2.1 Lagrange Multiplier Method The total potential energy of the system 7: is: It =7r, +7r2 +2, -(v—u,)+2, ~(v—u3) Where 7r, , 7r, are, respectively, the TPE of the first and of the second bar. ”=LI%(0'x'gx)dV+l,i-%(Ux'gx)dV—';f'°u’ +31 .(v—u,)+,1,,.(v—u3) dx (61) (62) fl=%'JEA[Z—x‘u]2+dxéijz/lztd—u)dx‘gfr'm+11'(V—“2)+’12‘(V‘”3)(63) ”—4-;—JEA(-:1ix—uldH-ijHA(duldx—Zfl.ul+fl1.(v_u2)+1’2.(v—u3) (64) dx Displacements are approximated linearly. u: uN =u(l—— +u — ,Zi " -' ‘ L 2 L Substituting: E A " dN dN E. A, 2" dN, dN at, Arm“ . (xi-Am rim-— 0 L l _I _P.u4+,l, ~(v-u3)+42'(V-u3) Setting the first variation with respect all the DOF to zero we get: 67! " dN dN, _ z E A . _' dx= 0 all, I l I( dx )[Zu “1 dx 1] 6a " dN dN ———:EA - 2 . I (1x =0 au. ' ' ll dx )[Zu dx) 2‘ 672' " dN dN —=EA- ——'— . 1 dx :0 au, 2 2 (ii dx )lzu" ax) A” (65) (66) (67) (68) (69) 0" an "fl/ dN —=E7A‘)‘ 2 .l —P=O 70 au, Kamila/ax] () a, a, A: < > QE=V—u,=0, git—=v—u,=0 (72) ax, ‘ 612 ‘ The resulting FE model is: F M —§fl'— O O O 0 O L L [ ‘ , W 0 f .M M 0 0 0 _1 0 ' L L “2 0 0 O EQA2 _ EA, 0 0 _1 u3 O L L < u4 r:( P } (73) 0 o — 53A? 152/12 0 0 0 v 0 L L 21 0 O 0 O 0 0 1 1 A, ,0 O —i O 0 l 0 0 i ’ ’ _ O O —1 0 1 0 0 4 Solving this system of equations: f. = —P (74) A, = P (75) A, = —P (76) u2=u3=v=—L—P (77) EIAI L u4 = — P + L P (78) EIAI EZAZ 3.2.2 Penalty Method In the penalty method the displacement continuity constraint is imposed through 26 two penalty parameters 7, and 7,. Thus the TPE of the system in this study changes its form. 71:77] 4.7-[2 +%7I -(v—u2)2 +%72 -(V—u3)2 where 72', , II, are, respectively, the TPE of the first and of the second bar. 4 flzfig-(o'x.g_,)dV+‘;[-;—(o'r.g,)dV—;fl.u, +%y, ~(v—u3)~+%}’2-(V—u3)2 I. 7 7_ (v— 14,) The linear Lagrange interpolation functions are used. 3 x x u: uN =u[l———)+u,[—j ; .l .l l L - L Substituting: _£_A_.-.’ EN. :13: _E_4_ LI: -2 I[,udxl][zu"dx]dx+2 I’qux 0 I I, J —P-u4 +%y, ~(v—u2)2 +-;—72 -(v—u3)2 The first variations of I: assume the subsequent forms: 53—: =E,,-A :K-d—N dx ][:u,%]dx=o 67: ” dN dN EzEIA'.J( dx2)[zuj dx —]dx— 7, (v— 212): 0 j 67! 1‘ dNl dN, ‘673:E2A2°J(Ex_)[zuj dx ]dx—72-(v—u3)=0 J 27 fl%.!E Am Jizdxj m%.23‘EA A,[::]dx meu+— y,(v— u,)+ l 5 2“]? (79) (80) (81) (82) (83) (84) (85) (86) I. dN 2’1:ng2. [dN2)[Zu ’de-P=0 (87) 6114 0 dx .1 1 dx 67: —=71‘(V’u2)+72'(V-u3)=0 (88) av Thus, the FE mode] is determined. — EIAI __ EIAI 0 O 0 _ L L , ‘ , ‘ E,A, E,A, 0 f. " L L + 7| 0 0 “7| J ”2 0 0 0 EzA2 7 _ EzA, fl], u3 )=< 0 } (89) L 7 L u4 P 0 0 — E°A° E24 0 (V, L0 L L _ 0 -7. m 0 7. + 72_ The solution to this system of equations takes the following form: f1=_P (90) u2=—L—P (91) 1A1 .=[L+_L__]p (92) 7| EIAI u3=[—L—+i+—l—]P (93) ElAl 71 72 u4=[—L—+i]P+[ L +—1—]P (94) EIAI 71 E2A2 72 3.2.3 Comparison between the Two Methods The convergence to the Lagrange solution when penalty parameters approach infinite is verified. 28 ——- The} Limitv=Limit —1—+—L— P: —-—L— P (95) ("m ("M 71 EIAI ElAl Limit u3 = Limit[——L—+-]—+-l—]P= {—L—jP (96) 5/1—2’2—W’ .VI-72 ‘WJ ElAl 7, 72 EIAI Limitu4=Limit i—+—1— P+ L +—1- P: L + L P (97) ""3"” ("7?”) EA 7| EzAz 72 EIAI E2142 The constraints we have enforced are: G, (u,,u2,u3,u4,V)=(V—”2) (98) G2 (u,,u2.u3,u4,V) =(v—u3) Then, as in the previous paragraph, if we define u,,,u,,,u3,,u4,,v, as the solutions derived from the Penalty formulation, it is possible to verify the well-known relation between the Lagrange multipliers and the penalty parameters: 1 L L = -G u ,u,,u ,u..,v = -v_—u. = - —+— P——P =P 99 ’11 71 1(ly -7 3y 4/ y) 71(/ 2f) 71 [[71 EH41] EIA, ] ( ) 2’2. :72 '02 (ulrauZy’u37’u4r’v7)=}/3 .(VV _u37): (100) :72.[[_L+——L )P-[i+'—L-+‘L]P]:_P y, ElA, 71 EIAI 72 3.2.4 Relation between Penalty Parameter and Beam Properties Like in the single beam case, the exact theoretical displacement of the tip of the beam matches that from the Lagrange multiplier finite element formulation. uj"""’=—L P+ L P (101) E,A, E2A, The Penalty parameter solution uf"”"”" differs from the exact one due to the presence of 29 two additional terms (P/ 7,) and (P/ 7, ). uf"""”-"= -—L——+i P+ L +i P (102) E|A| 7| EzAz 72 U .W I to u4 [7c nu/n The ratio of u is evaluated in order to underline the dependencies between the two solutions. wally :[ L— + 1]P+[2L;4_2+ - ljp 71-7-1] ”4L ' EIAI 7| EA 77 1 7| 77 = + ’ 103 ugl’ml [LLA _, L ZJP ( ) El Al +E2/4 Now, if the penalty parameters 7I and 7, are substituted by: E|A| Ll-IBLLILW] (104) 72 —fl-( L j (105) the ratio between the solutions becomes independent of materiaI and geometrical properties of the beams. Pt naln u4 l.’ .rat‘l 4 =1+L 106 ,3 ( ) 3.3 Bending Loads — One Element — Full Integration The conventionally formulated linear finite element model of the Timoshenko beam is generally not used in commercial finite element codes due to its shear locking behavior. Instead, a Consistent Interpolation Element (CIE) or Reduced Integration Element (RIE) is adopted. For this reason, only the 1 element case is studied with a 30 Fl; [01 a} U335, ihe f conventional FE formulation. The ptu'pose is to compare the results with those generated by the Reduced Integration Element. Now, we will focus our attention on the deflection of a beam under typical bending loads: transversely distributed load F, concentrated moment M and concentrated force Q applied at the tip. 7% F Q t y, 7r . w ft. x r ., r, [V ,5. ,V l' T' .. n v rv v ‘|. ._ . I ‘ ‘ 7‘ ‘ ' . , .‘ 1 | |<|:,. w.“ . ; ‘ ' . I | , 7 .. 1 . , "" ‘ , . . :i 1“ ,‘ix‘ ] ’»"ifh"’jl‘ “1,1‘ ‘ l EI‘ ' L. .} 7|.l v | :"'l NH: ‘7! |.. l l. l 7n r 'I A L ‘ 1 I' I All 11‘ I Al' ii '1 V 1 E, I, k, G,A 2 o w, w, ¢ .. W2 Y’v 9’: L '16 Figure 4. Beam under bending loads — One Element This case is analyzed using a single linear Timoshenko beam element connected to one fixed point V by a displacement continuity constraint; this is imposed through a Lagrange multiplier or a Penalty parameter. The configuration of the geometry and the mesh is depicted in Figure 4. 3.3.1 Lagrange Multipliers Method The hybrid interface method introduces two Lagrange multipliers A, and A? in the total potential energy (TPE) of the system to satisfy the continuity conditions for transverse displacement and rotation. Thus the TPE of the system in this study assumes the following form. 31 7t =7r, +11 -(w‘, —w,)+/l2 -(‘I’v —‘l’,)—rl -w‘, —ml *1", (107) Where 7r, is the TPE of the bar; r, and ml are, respectively, the vertical reaction force and the reaction moment at V. 1 as! V, 5w 2 " 2 2 7T] =—2'I[E1(a )2 +kGA(LIJ+—a—-) :ldx—aflF-W]dx—;Ql.wl —;M’ -\P, (108) 0 x X The first variations of 7rl take the following form. a —j[51[9"_’]§‘3‘1+MG/1[ iw.)(ay+9‘i‘£]]dx— 0 6x 6x ax ax (109) RF. 6w)dx 2.1—Q w— 2M ‘I’ O The Lagrange linear interpolation functions are used. 2 x x w: wN =w l——- +w, — 110 Z " " i L) (L) ( ) 2 LP=ZLPJNI=LP,(1—i)+w,[i) (111) 1:1 i L b L Thus, 7:, can be written: 67!, = I. [{EILJJ 0 QMZM ‘21:] kGA[ +L_I{kGA[;‘PJNJ 7.;ij %](26w %J}dx— dillzwlh [(F-Z6w,-N,]dx— 0 (112) 0 —in .WI —iMl .‘Pr l=I i=1 Taking the first variation of 7r with respect all the DOFs and setting them to zero: J ——=j{kGA[Z~P,N‘W 1 dx +2 wjddx ‘Z‘xl—‘Ndx—[fingdx—tzo (113) 32 672' 1’ dN dN “W __= 131 111 __7__2 +kGA ‘IJNN + —"—N dx- =0 114 611’. Ii [22:2,] [ZMZWWN A? () L dN I. a”=II"GAlZ~P,~F’.—Zz-+ij ’szlidhlw-Nflmwo “157 0 .l 6w, _, dx dx 0 = E] ~11 J 3 +kGA ZWNN3+ZW 'N, dx—M=O(116) 6W2 6.7 :1: 1 4" dx .1 j I] ,/ J dx ;”=L-q=o (117) -§Z—l=wv—Wl=0 (119) §E=wv 41’, =0 (120) Then, the resulting FE model is: kAG _kAG _kAG _mo 0 o —1 0 FL L 2 L 2 _ _ _(_] MG 151 kAGL kAG kAGL E] W. 2 —— —+— —— 0 0 0 —| 2 L 3 2 6 L ‘1’. 0 “kg 549 ELG 1g 0 0 0 0 W2 _Q_fl L 2 L 2 111, 2 (121) _kAG kAGL_fl kAo _E_I_+kAGL 0 o 0 0 W. M 2 6 L 2 L 3 q," ,1 0 0 o o o o 1 0 o 0 0 o 0 1 A“ m' A. o —1 o 0 0 1 0 0 0 ‘ “ 0 _ 1 _ 0 —| o o 0 | o 0‘ Taking into consideration that V is fixed (W, = 0 and ‘1’, = O), the system can be solved. For the purposes of this study it is convenient to identify contribution to the solution 33 coming 1 1117111} 111 Instead. 1 changes I coming from the three types of loads. If only the concentrated moment M at the tip is present, the solution is: W _ 6ML2 3 12EI+kAGL2 _ 12ML2 2 12151 + kAGL2 2:0 L=—M b (122) (123) (124) (125) (126) (127) Instead, when only the concentrated transverse load Q at the tip is applied, the solution changes to: _ (12EIL + 4kAGL3 )Q _ _ kAG(|2EI + kAGLz) : 6ng 2 (1251 +kAGL2) 31=Q 32=-LQ Finally, a transversely distributed load F yields the following outcome. 34 (128) (129) (130) (131) (132) (133) (134) (135) In 61311 11111 3.3.2 P1 In lhIOU‘gh {u becomes: Selljp . (I (65le + 2kAGL“ ) F W7 :— '1 (136) 12EIkAG + (kAGL)‘ 3 2 = “F , (137) 12EI+kAGL“ 1, = LF (138) 1, =—L—“2F— (139) In this case the solutions provided by the hybrid interface method are not the exact theoretical ones. 3.3.2 Penalty Method In the penalty method the displacement continuity constraints are imposed through two penalty parameters 7l and 72. Thus the TPE of the system in this study becomes: (140) V I 2 1 2 ”=”1+§71‘(Wv—W|) +§7z~(‘l’v—‘Pl) —r,-wv—ml-‘IJ with 7:, equal to: 6x x it, =%I[EI(§L£) +kGA(‘I’ +gfi] de—IIF'WIdx’2QfiW1’2My‘V. (141) The displacements are approximated linearly: w=Zw1Nj =wl (l—£]+w2[—x-) (142) F. - L L ~17 =ij1v, =91, (1—1)+w,(1) (143) 1:1 L L Setting the first variation of a with respect to all the DOFs to zero we get: 35 (Y"’lt't‘1z|'.- ,l J a” " dN dN dN " —= kGA ~111v——' ___,__. dx— F-N dr— aw, (Ii [2 1 "13+ij dx ax” K ‘) (144) 1. -w(. - > = o ./ 672 dN dN dN _= E] \P J——'— +kGA ‘PNN+ ——’N air— I. 67: L dN2 dN‘ dNZ . 11101271721 1.11 dx 0 1, dN dN 5” =I 5, 21,714sz 3.10/1 ZW/N/N,+ij—’N, dx—M=0(147) 0 _/ dx (ix _/ J dx 671 —= -w‘,-—w —r=0 148 aw" 7| ( 1) 1 ( ) 672' 'aTPTZ72'(LP1--LP|)_”7|:0 (149) Then, the resulting FE model is: WEE. + __kA_G_ _"AG ""142 _ 0 _ — - L 7' 2 L 2 7' _ _ {fl} _ kAG £31 + kAGL + kAG kAGL _ £1 0 W1 2 2 L 3 72 2 6 L 7 ‘1’. 0 1040 MC: kAG kAG w1 FL —— — — 0 0 “ = — — — L 2 L 2 ‘1’, Q I 2 I (150) _ kAG kAGL __ fl 1040 g1+ kAGL o 0 w‘ M 2 6 L 2 L 3 )1,” ,1 —7I O O 0 7l 0 m _ 0 —7: 0 0 0 7:, ’ l ‘ Knowing that V is fixed (w‘, = 0 and ‘1", = O ). the system is solved. Concentrated moment M: wl = O (151) 36 l1. l11, = (152) (IZEIL + kAGL‘ + 61.27, ) M ”I, : — 153 " (12EI+kAGL2)7, ( ) (12EI+kAGL2+12L7,)M 154 _ (12E1+kAGL2)7, ( ) Concentrated transverse load Q w, = —9 (155) 71 711, : _£ (156) 72 2 2 2 3 [12E1L7,7, + (kAGL) -(L 71 + 7,) + 4kAG -(3EIL 7, + 31517, + L 7,7, )]Q w, = _. 157 “ kAG(12EI+kAGL2)7,7, ( ) 12EIL + kAGL3 + 6L2 ~11,=( 2 7,)Q (158) (12E1+kAGL )7, Transversely distributed load F: w, = —E (159) 71 ‘1’, = FL‘ (160) 272 w2 = [12E1L7,7, +(kAGL)2 -(L27, +27, )+4kAG.(3EIL27, +6El7, + 137,7, )]LQ (161) — 2kAG(12EI+kAGL2)7, 7, (12EIL+ kAGL‘ +6L27,)LQ l11, = 2 ‘ (162) 2(1251 +kAGL )7, 37 3.3.3 Comparison between the Two Methods When the penalty parameters approach infinity the Penalty method solutions tend to the Lagrange ones. Only the outcome regarding node 2 will be considered for brevity. Concentrated moment M: 1251L+kAGL3 +6L2 M _ 2 Limit w, = Limit —( 72) = 6ML , (163) we» - rz-m (1251 + kAGLZ )7, 12E] + kAGL- 1251+1cAGL2 +12L M 2 Limit‘l’, = Limit 72) = ”ML , (164) 73+.» ~ 72% (1213] + kAGLz )72 12E] + kAGL‘ Concentrated transverse load Q 1251L + kAGL‘ + 6L2 , 2 Limit‘l’, =Limit ( , HQ = 6“) , (165) mm 7),... (1251 + kAGL‘ )7, (1251 + kAGL) T ransversely distributed load F: Limit w2 = rlnyz—MD 1251L7,7, +(kAGL)2 ~(L27, + 7,)+416<1G.(351L27l +3517, + L37,7,)]Q Limit —* 166 7.7.72 kAG(12E1+kAGL2) 7,7, ( ) (1251L+4kAGL-‘)P ’— kAG(12EI+kAGL2) 12515 + kAGL3 + 6L2 L 3 Limit‘I’, =Limit ( 2 7,) Q = 3“) 2 (167) 7.72% 7mm 1 2(1251 + kAGL )7, (1251 + kAGL ) As before, the relations between the Lagrange multipliers and the penalty parameters are verified. The constraints we have enforced are: G, =(w,, —w,) (168) G, =(L11, 411,) (169) 38 (1111111111 ____L (7.7111111! IILIIIVI ‘1’ I' \— 33'4 R1 A) material a ——_ ‘ Concentrated moment M: A,=7,-G,=7,-(w,,—w,)=7,'[O—O]=O (170) M ’72:72°Gz:72'(va—kpl):72'l:0_—]:_M (171) 2 Concentrated transverse load Q: 41=7rG|=71'(Wv-“’|)=7|'[0—[—%]]=Q (172) l 62:72'02:72'(qjv-qj|):72'[0—[LQH:_LQ (173) Transverse/v distributed load F: . L A=7.~o.=r.-(w..—w.)=r.IO—[——7—II=FL (174) 15L2 15L2 . =1"3:7' ._\p :3' — =—— 1, 7- 0- r- (*P. .) 1- [0 [272]] 2 (175) 3.3.4 Relation between Penalty Parameter and Beam Properties As before, the obtained solutions are examined in order to find a relation between material and geometrical properties of the beam and the penalty parameter. The exact theoretical solution for displacement of the tip of the beam does not match that from the Lagrange multiplier finite element formulation. Concentrated Moment M: The ratio w,”"""”-" to w,“’“"””"" is evaluated in order to underline the dependencies between the two solutions. 39 3.0“. II Thus. Paramuu Then (I: (".1 ,—. (1251L + LAGL3 + 61.27, ) M w,”""""" (1251 + kAGL2 )7, (1251 + kAGL2 ) [hagrunge : 2 = 1 + (176) Wz‘ _. _§_/!£_ m 6L - 72 1251 + kAGLZ Now, if the penalty parameter 7, is substituted by: (1251 + kAGL2 ) n=fl~ (nn 6L The previous ratio becomes independent of material and geometrical properties of the beam. ) pend/(1' 1 —2———=1+— 178 ,6 ( ) L'xut't 2 Thus, the accuracy of the solution depends directly on the value assigned to the parameter 6. The ratio ‘I’§"’"”"" over L11,"'*"“"*"’ is: (1251 + kAGL2 +12L7,)M pend/'1' 1251 + kAGL2 1251 + kAGLz “’7 ,,= l , I“ =1+< I (179) \qugnmgt IZML 12L . 72 1251 + kAGL2 Setting: (12EI+kAGL2) 180 72 — fl ,2, ( I We get the same relation found previously. ' =1+i (181) Concentrated Transversal Load 0: 40 X1111. if ”1613111) The Mill) 5911,“, h ,0“, I .ugrungc The ratio 117,”"""/’-" over w, is: penalty 2 lugrwtgt' _ ”’2 2 [12EIL7,7, +(kAGL) -(L37, + 7,)+4kAG ~(351L37, + 3517, + 57,7, )]Q ”FT—7 T i “7 Qéhfiifloflig“ — “ (1251L+41AGL3)Q kAG(12EI+kAGL2) kAGLz(12EI+kAGL2) kAG(12EI+kAGL2) + (12E1L+4kAGL3)7, +(12EIL+4kAGL3)7, Now, if the penalty parameters 7, and 7, are substituted by: 2kAG(12E1+kAGL2)| 71 = 16 ' 3 (1251L +4kAGL) 2164a2 (1251 +kAGL2) 72 = fl. 3 (1251L+4kAGL ) The ratio between the solutions becomes: “1' ,nenultt~ 1 1 The ratio ‘I’f'mm' over WWW" is: 1251L + kAGL3 + 6L27, Q ( ) ,,,,...,1. (1251 + maL2 )7, 1 (1251 + kAGLZ) : = + \IJ;"’X’"”X“ W 6 L2 Q _ 61472 (1251 + kAGLz) Setting: (1251 + kAGLz ) 72 = fl ' 6L It follows: 41 (182) (183) (184) (185) (186) (187) 71111111 51711111 Ill? 121110 The 1211i Scltln_ {I} fund/(y 1m, =1+—]— (188) ‘1’: ‘ fl T ransverselv Distributed Load F: The ratio w,”"""”—" over w,’"’”"”“" is: “fixnalty wéugrwtge : [1251L7, 7, + (kAGL)2 (L37, + 27, )+ 4kAG ~(3E1L2 7, + 6El7, + L37, 7, )] FL “NWT—7M 2kAG(1251+kAGL2) 7,7, _ 7 (189) _ (651L + 2kAGL3 ) FL _ kAG(12EI + kAGLZ) kAG(2451+2kAGL'-’) kAGL2(12EI+kAGL2) + (12E1L+4kAGL3)7, + (12EIL+4kAGL3)7, Now, if the penalty parameters 7, and 7, are substituted by: 2kAG(24EI + 2kAGLz) ' (1251L + 4kAGL3) n=fl (mm 21cAGL2 (1251 +kAGL2) ' (1251L + 4kAGL3) 72 =fl (191) The ratio between the solutions turns into: 779W" 1 1 1 =1+—-+—=1+— (192) W2" 216 216 fl The ratio ‘I’f'm’l" over W,”‘“"”"”" is: (1251L + kAGL3 + 6L27, ) FL WW. (1251 + tcAGL2 )7, _ 1 (1251 + kAGLZ) 193 ' = = + (pguxrangt' 3L2 F 6 L72 ( ) I (1251 + kAGLi ) Setting: 42 11 117111.11 3.4 B 3.4.1 1 total 1701: displaccn Iltllmtint: (1251 +kAGL2) , = . 194 7. 16 6L ( ) It follows: penalty lillt‘xutl : 1 + i (195) ‘1’; fl 3.4 Bending Loads — One Element — Reduced Integration 3.4.] Lagrange Multipliers Method The hybrid interface method introduces two Lagrange multiplier A, and A, in the total potential energy (TPE) of system to satisfy the continuity conditions for transverse displacement and rotation. Thus the TPE of the system in this study assumes the following form. 71:71, +/l, ~(w,,—w,)+/i,-(‘I’,,—‘I’,)—r, -w,,—m, -‘I1 (196) V 1" 6912 aw2 " 2 2 7r=-— EI— +kGA LP+-— dx— F-wdx- ~w— M-‘P 197 .,01[(,x) ( ,le 011 1 go: .1) The first variations of 71 is: I. I. 67:, = EI(§P—Jfl+kGA[W+iw—XW+@J dx— ((F-aw)dx— 6x 6x 6x 6x 0 0 :2 'W. -EZ:M. "P, 1=l 1:1 (198) The Lagrange linear interpolation functions are used to approximate L1’ in: L 2 1 (5(a) (d. (199) 2 0 6x 43 Instead, in the expression: %I[kGA (11 + 9‘1) (at): (200) the term [‘1’ + 52) will be substituted by ( ‘I’,+‘P, w,—w,) + . x L This operation will produce the same finite element model as the Reduced Integration procedure. Thus, 67:, is: on, =[({51[Z~11, dgx1](za11%]}ax+ 0 J I +I{kGA(lP‘ +172 + W? ’ “IX-[6711. +1681, +-1—6w2 —l6w.)}dx— (201) 0 2 L L I L 2 2 —([F-Zaw, ~N,]dx—iQ, -w, —>2:M, A11, 0 Taking the first variation of n with respect all the DOFs and setting them to zero we get: 671 w ‘I’ w ‘I’ 1' ———=kGA _L__l__._i___2. _ F-N dx— :0 202 817, [L 2 L 2I 0k ‘) 1‘ ( ) t 7 I’ dN _a_7_r_:kGA[_ll_+£i+h+—L\P2)+I El Z‘P,-—Jfl (Ix—312:0 (203) 691, 2 4 2 4 0 , dx dx I. 5” =kGA(_m+fl+k+31-_)+((F.N,)dx+Q=o (204) (31., L 2 L 2 0 “ 1 dN 5” =(-KL+-——W'+1"3—+sz+l 51 281, 1sz 611, 2 4 2 4 , dx dx 0 Hair—M = 0 (205) 67: R=41fli=0 (206) 44 =55. Then, the resulting FE model is: rI610 142 L kAG 2 MG L kAG OOOOIJ O COCO O—‘OOO 1-‘OOO 2E] 66 3.7 11 mcnm 3) 10 Cl] 3.7.1 ammb paramc‘. CIIIISCL‘ U The im Ol fish (lnl 3. 7 Interface Element The work done so far gives us the tools to achieve two different tasks: 1) to define the stiffness matrix of an “Interface” element for the linear Timoshenko beam elements; 2) to choose automatically the best values of the penalty parameters. 3.7.1 Stiffness Matrix We first consider the uniaxial load case. We observe the stiffness matrix of the assembled system (Section 3.2.2) composed by two beams connected through the penalty parameter. We notice that the matrix could be seen as the result of assembling three consecutive elements. — EIAI _E1Al 0 0 0 - L L f N [f3 E,A, EIA, 0 1 " L L +71 0 0 ’71 “2 0 0 0 52A, + 7 _E,A, _y lu3 )=< 0) (398) L " L 2 u4 P 0 0 —§3i 153/13 0 IV) l0. L L L. 0 —7I -72 0 71 +72_ The intermediate interface element would be constituted by three points with one degree of freedom each, and would have the following stiffness matrix and displacement vector. 71 0 ‘71 ”2 O 72 —72 u3 (399) ‘71 _72 71+72 v Changing the order of the axial displacements, it is: 3 | 3 o .F (400) I 3 3 + .3" l 35 C dement andloac nodes \)i F'gure 6. TO “The - pai'iimclcr Table 1. The same path can be followed starting from our study of the two linear beam elements connected using a Penalty parameter continuity constraint, fixed in the center and loaded with bending loads (Section 3.5.2). In this case the outcome is: F 7, 0 —7, 0 0 0 1 rwI l 0 7: 0 —7_, O 0 ‘1’, -r. 0 7. +73 0 -r. 0 < W. ( (401) 0 -72 0 7: + r. 0 -7. ‘1’. O O —73 0 y, 0 w2 _ 0 0 O ‘74 O 74 (KP: l Now, if we join the two stiffness matrices, an interface element made of three nodes with three DOFs each would be defined, as in Figure 6. 1 V 2 ¢ 0 a “1,W19 ‘1'1 “v,Wv9 Tv “2,W2a ‘1’2 Figure 6. Interface element for the T imoshenko beam element To write the new 9x9 stiffness matrix, we need first to redefine the name of the penalty parameters to be associated with each of the constraints. Table l achieves this task. Table 1. Penalty parameters associated with each of the constraints enforced Penalty parameter Constraint enforced 711 (”v ' ”1) 72, (w, — w,) 7 31 (WV ‘ ‘1") 712 (u, — u,) 722 (w, " W2) 7 32 (‘1’), " KIJ2) 68 Thus, the stiffness matrix and the displacement vector of the interface element are: '- _ f y” 0 0 ‘711 O O 0 0 0 ul 0 72, O 0 —}/3, 0 O O 0 W1 0 0 73' 0 0 ‘731 0 O 0 ‘1’] ‘711 O 0 711+712 0 0 "712 0 0 "v 0 ‘721 0 O 722 +721 0 0 ‘722 0 < WV 7 (402) 0 0 ’73] 0 0 731 +732 0 0 —732 L111- 0 0 0 ‘71: 0 0 7|, 0 0 u2 0 0 0 0 ‘722 0 0 722 0 W2 _ 0 o 0 0 0 m. 0 0 732 , l‘l’z. 3.7.2 Best Values of the Penalty Parameters Consider the problem, plotted in Figure 7, of two linear Timoshenko beams connected through an interface element. 4—0 —0 o—> a Fixed: uv= 0 1 El9ll9klaG19Al 2 V 3 E29129k29G29A2 4 ul L1 uz uv “3 Lz u4 Figure 7. Beam under axial loads — Two Elements Results obtained in sections 3.1.4 and 3.2.4 lead us to affirm that the penalty parameters 7,, and 7,, should be: r” = fl[——] (403) 69 15,21, 712 - fl[Tl (404) The ratio between the solutions of the problem becomes independent of material and geometrical properties of the beams. u penulti- penalty 1 ' , . = 4 , = 1 + — (405) am! exact u 1 “4 26 Based on the results of sections 3.4.4 and 3.5.2, for the other two DOFs, w and ‘l’ , a similar conclusion can be reached for each kind of bending load. However, in this case we cannot find the values of the penalty parameters that make true both of the following equations: H’ penalty , perm/11‘ 1 l _ 4 __ __ ,Iugrangc — lama/age - I + (406) ’1 W ,6 ' In- cnaln‘ LPPU'“ . {PP . 1 l _ 4 _ _ LPlugru/rge — LPLugra/Ige _ l + (407) 1 4 fl We can satisfy exactly only one of the two equations, and assign to the other one a value of the penalty parameter greater than required. Concentrated Moment M In order to get equation (406) to be true in the presence of a concentrated moment M, we need the following. r31 = fl - —2]21‘ (408) 7.. = .6 - 3%5- (409) Instead, to fulfill the relation (407), it is required that: 70 A tho option 35 a c- E111 731 :28“ L1 EJ7 732 :26? ‘- (410) (411) A choice needs to be taken between the two groups of two values. The most reasonable option is to take the higher value and satisfy (406), that is: 2E1 731:.B'—Zl'_l2 732 :fl 1 As a consequence, this selection determines the next equation to be true. penalty penalty ‘1’, _ ‘11, Lagrange _ Lagrange _ w, ‘1’, Concentrated Transverse Load 9 Equation (406) requires: 2 721 = :8 -13..- +__!2 4E111 kl/llGl 21.2 73 = 16 ' I ' 13: ,_g_ 4EII1 kl/llGl 2 722 2 fl 1%. + L2__ 4E,l2 [(214sz 21.3 732 = )6 ' b _ 1L + £22 41:72]2 k2A262 While equation (407) specifies: 71 215,1, (412) (413) (414) (415) (416) (417) There “hich onh‘0' 11101111. Transi Itquhe \the e 31=18-£31L (418) 732 :fl'EZIZ (419) There are no simple relations between the two expressions for 73, and [,2 to tell us which to choose. Then, it is convenient to compare their values and take the bigger one only when the material and geometrical properties of the two beams being connected are known. Transversely Distributed Load F Results for this load case are similar to those of the previous. Equation (406) requires: 2L 721 = l6 . 1:: + I L12 (420) 8E,ll 2klAlGl 22: 4 731 Z 16 L: + L? ”— ( 21) 4E]!l klAIGl 2L2 (422) 215 (423) 2.. = fl-fl (424) 72 (425) Again, no simple relations exist between the two expressions for 73, and 732 to tell us which to choose. When the material and geometrical properties of the two beams is known the bigger one can be identified and chosen. Table 2 summarizes the all relations among the penalty parameters and model characteristics. Table 2. Penalty parameters relations for a T imoshenko beam element. Axial load Moment Transverse load applied Transversely applied at at the tip Distributed Load the tip E A 2.. ”'lfl - - - 2 2 fl ' 3 ’6 . 4 LI 2 ywl - - 1’1 + [’1 1’1 + L1 4E1], klAlGl 8E1], 2k,A,G, 2 2 2 3 p z: 11 2, fl ' z: 1" 1: 2Elll [ + J [ + J 71/1 - )6 - _L] 4E,ll klAlG, 4E1]I klAlGl or fl-§'—I' or fl-—2E'I‘ 14 14 , 2,2, 4 1 . l 2 2L fl 3 '6 ' 4 2 2 7142 - - [’2 + L2 1’2 + [’2 4E2]2 szsz 8E2]2 2k2AzG2 2 2 2 3 '8 3 L2 ’6 ' 4 I? 2 1., 2., 2, 2., 2E2]2 + + n2 - ,8 - L, 4E2]2 k2A2(32 4E2]2 szzG2 ” E212 2E2]2 or n .— fl 1’2 fl 2 73 3.8 Numerical Results To test the behavior of the relations just obtained for an automatic choice of the penalty parameters, the developed Interface element has been implemented in the commercial finite element code “ABAQUS”. The User Element Subroutine (UEL) receives all the necessary information about geometry and material properties of the two connected Timoshenko beams from the input file. Then, the stiffness matrix of the interface element is built. Inside the subroutine, that is written in F ortran 77, a series of “IF” structures select the appropriate values of penalty parameters in all cases where a choice is required. 3.8.1 One Element — Clamped Beam We will test first the resulting displacements of a beam under: axial load P, transversely distributed load F, concentrated moment M and force Q applied at the tip. 1 F _ Q \ D d P M V 1 E,I, k, G,A 2 O G O llvana‘I’v u1,w1 ”’1 L “zawza‘l'z Figure 8. Beam under bending loads — One Element 74 This case is analyzed using a single linear Timoshenko beam element connected to an undeforrnable structure by an interface element. The configuration of the geometry and the mesh is plotted in Figure 8. The beam properties are: Length=0.5, thickness=0. l , width=0.1. Young’s modulus=200.E9, shear correction factor=0.85 In Table 3 the results for various combinations of load are presented. The word “Abaqus” marks the outcomes obtained without using the interface element, that in this case means to set w, = 0 and L1’, = 0. Instead, the word “UEL” indicates the use of the interface element. Table 3. Single beam clamped - Outcomes w2 Abaqus w2 UEL w,“”""““ my” ‘1’, Abaqus LIt, UEL w;"‘"""“' #1137" P 5.0000E-08 5.0050E-08 9.9900E-01 O 0 1 M -3.0000E-03 -3.0030E-03 9.9900E-01 -6.0000E-03 -6.0030E-03 9.9950E-01 Q -1.7942E-03 -1.7972E-O3 9.9833E-01 -3.0000E—03 -3.0015E-03 9.9950E-01 F -8.9712E-04 -8.9937E-04 9.9750E-01 -1.5000E-03 -1.5008E-O3 9.9947E-01 M+Q 4.7942E-03 -4.797ZE-03 9.9937E-01 -9.0000E-03 -9.0020E-03 9.9978E-01 M+Q+F -5.6914E-03 -5.6941E-03 9.9953E-O1 -1.0500E-02 -1.0502E-02 9.9981 E-01 We assumed all the loads to have the magnitude 1 and ,6 = 1000. 3.8.2 Two Elements — Clamped Beam Consider the problem wherein two linear beam elements are connected using an interface element and clamped at the tip of one beam. The configuration of the geometry and the mesh are plotted in Figure 9. 75 Figure 9. “2 "V . w w. y 9’2 Y’» “3 W U4 4 4w; 5”: 3 E,I,k,G,A 3: Two beam elements joined through an interface element The two beams are made of the same material, whose properties are: Length=0.5, thickness=0.1, width=0.1, Young’s modulus=200.E9, Y’4 shear correction factor=0.85 Table 4. Two beams clamped - Outcomes w, Abaqus w4 UEL w,“”"""“ / w!” ‘1’, Abaqus ‘1’, UEL ‘1’j”“"“"/‘Pf,””' A 5.0000E-05 5.0050E-05 9.9900E-01 O O 1.0000E+00 M -3.0000E-O3 ~3.0020E+00 9.9933E-04 -6.0000E-03 -6.0030E+00 9.9950E-04 Q -1.9486E-03 -1.9494E-03 9.9959E-01 -3.0000E-03 -3.0008E-03 9.9973E-01 F -7.8682E-O4 -7.8710E-04 9.9964E-01 -1.1250E-03 -1.1252E-O3 9.9982E—01 M+Q -4.948GE-03 4.9497503 9.9978E-01 -9.0000E-03 -9.0015E-03 9.9983E-01 M+Q+F -5.7355E-03 -5.7362E-03 9.9988E-01 -1.0125E-02 -1.0126E-02 9.9990E-01 In Table 3 the results for various combinations of load are presented all the loads to have the magnitude 1 and ,6 =1000. 3.8.3 Two Elements — Two Materials — Clamped Beam . We assumed Consider two linear beams elements connected using an interface element and clamped at the tip of one beam. The configuration of the geometry and the mesh are plotted in Figure 10. The two beams are made of two different materials, whose 76 properties are: M Length=0.5, thickness=0.l, width=0.1, Young’s Modulus=200.E9, shear correction factor=0.85 £96111; Length=0.5, thickness=0.05, width=0.05, Young’s Modulus=1.E9, shear correction factor=0.85 7 1 Q % O "I “2 “V "3 u4 w, : E12 In k1, Gr, A1 2w; Wv :7 W32 E22 12, 1‘2; 022 A2 1 w4 an L, 116 W. 116 L2 21,, Figure 10. Two beam elements joined through an interface element In Table 5 the results for a concentrated force Q of magnitude 1 applied at the tip of the beam are presented for various values of ,6. Table 5. Two beams — two materials - clamped - Outcomes w. Abaqus w. UEL wi”"""“/wi‘"’" ‘P. Abaqus ‘P. UEL 91:”“18 7942’” fl =1E3 -7.2408E-05 725295-05 9.9833E-01 240225-04 240355-04 9.9946E-01 fl=1E4 -7.2408E-05 724205-05 9.9983E-01 240225-04 240245-04 9.99925-01 fl=1E5 -7.2408E-05 -7.24105-05 999975-01 240225-04 240235-04 9.9996E-01 3.8.4 Three Elements — Simply Supported Beam Three linear beam elements are connected using two interface elements. All the 77 system is simply supported. The configuration and the mesh are plotted in Figure 11. F 1 .11 1.11!1111111111111111151111111 1111111 1111111111 4 2 I I3 nterf. Elm. l nterf. Elm. 2 Figure 11. Three elements simply supported beam The three beams are made of three different materials, whose properties are: 1310.421 Length=0.5, thickness=0.15, width=0.15, Young’s Modulus=400.E9, shear correction factor=0.85 Beam; Length=0.5, thickness=0.1, width=0.l, Young’s Modulus=200.E9, shear correction factor=0.85 2901722 Length=0.5, thickness=0.025, width=0.025, Young’s Modulus=100.E9, shear correction factor=0.85 In Table 6 the results for a distributed force F of magnitude 1 are presented for ,8 = IE 4 . Table 6. Three elements simply supported beam - Outcomes w4 Abaqus w4 UEL wfbaqm / wf,“ ‘1’, Abaqus ‘1’, UEL WM“ 1111:,“ Node 1 0.0000E+OO 0.0000E+00 1.0000E+00 -3.8741E-06 -3.8746E-06 999875-01 Node 2 -1.9364E-06 -1.9366E-06 9.9990E-01 -3.8704E-06 -3.8709E-06 9.9987E-O1 Node 3 -3.8529E-06 -3.8533E-06 9.9990E-01 -3.7954E-06 -3.7959E-06 9.9987E-01 Node 4 0.0000E+OO 0.0000E+00 1.0000E+00 1 .5405E-05 154055-05 1.0000E+00 78 CHAPTER 4 PLANE STRESS QUADRILATERAL ELEMENT 4. 1 Finite Element Model developed. The plane elasticity problems are described by two coupled partial differential (426) In this chapter an interface element for plane stress quadrilateral elements is equations expressed in terms of the two components of the displacement vector u and v. 611 6v 8 ea 8v 6[ [ -—(Cll_+CIZ—)—— C66 _+'— (3x 6x 6y ' 6y 6x a“ 3:8 (427) 6 6a 6v 6 __ C66 ’T—‘IT— __(Cl2_+ 22 6x 6y fix by ‘ 6x ” by For an isotropic material the constant terms c” are defined in the plane stress stress-strain relations to be: c - c - Y ll 22 1 _ V2 c 12 V CH (428) Y Where Y is the Young’s modulus. Appling the total potential energy principle, the weak form of the two governing 79 equations can be developed. The associated finite element model is: Where: 6N1+ 6N K”: QJ‘h[CH— 0N 6N, A’ 6A1 K,,_ “[612 a__Naay 8N (1+5N MON 141 1421 . (F) 141" [4221 11121111241 ax ax ‘4 8y 8y ‘6‘; +6666an Ki' = [K5217 665365ij 22ayay where h is the thickness of the structure. —’] dxdy —-——’—) dxdy ——’] dxdy (429) (430) (431) (432) (433) In the following, the stiffness matrix defined above is computed for a rectangular element having the in plane dimensions shown in Figure 12 and a thickness equal to h. Figure 12. A Geometrical dimensions of the rectangular element. Assuming the Poisson ratio v to be zero, and the displacements to be approximated through the linear Lagrange interpolation functions, the matrices K U are: 80 2 -2 —] l 2 I —1 —2 K 11 Y b —2 2 l -—l + Y a l 2 —2 -1 “‘ 6a —1 1 2 -2 2 6b -1 —2 2 1 " l —l —2 2 -2 —1 I 2 (h(aY +bY) h(aY _bY) h(_aY _bY) h(_aY +bY), 6b 3a 12b 3a 12b 6a 6b 6a aYfi __ (21' al’, bY h __qr [)7 _ar _b,Y, h(12b 3a) h(6b 3a) ( 6b 60) h( [2b 6a) 12b 6a 6b 6a 6b 30 12b 3a [1 __al' [2)” h _aY _bY h 01’ _bY h (11’ bY K ( 6b 6a) ( 12b 6a) (12b 3a) (6b 3a)) (434) 2 -2 —1 1 2 1 —I —2 K Y b -2 2 I —1 +Y a l 2 -2 —I 7 ._ 2 - Z *2 2 6a —1 1 2 -2 6b —1 —2 2 1 l —l ——2 2 —2 —1 1 2 (h(aY br') h(aY_bY) h(__aY _bY) h(_aY bY)) 3b 6a 6b 6a 6b 12a 3b 12a 71(0), b)’ h(aY bl’) h(_aY bY) h alf _bY) 6b 6a 3b 6a 3b [20 6b [2a h(_aY _bY h(_aY +br’) h(aY bY) 71(0), _bY) 6b 12a 3b 12a 3b 6a 6b 60 h _qr pr h _aY br 11 aY _bl’ h at bl: K ( 3b [2a) ( 6b 12a) (6b 6a) (312 6a)4 (435) f 1 l _l _1 1,11 (hY _hY _hY 711’ 1 4 4 4 4 8 8 8 8 Y 1 _/ _1 l [11’ _her _liY hY _ 7 4 4 4 4 _ 8 8 8 8 K’2"h 2 _/ 1 1 _1 ‘ _hY h)” hY _l_1Y 4 4 4 4 8 8 8 8 l _1 ,1 l _j _h1’ hi’ hi’ _hY| K K 4 4 4 414 K 8 8 8 1 (436) Finally, the stiffiiess matrix of the rectangular element can be defined. (h(aY+bY) h(aY_hY) h(_aY_blf) h(_aY+bY) th _hY _hY hi" \ 6b 30 12h 3a I2b 6a 6b 60 8 8 N 8 1 al' bl" al' 111' 'a)’ hl’ ar_br 111’ _hr 711’ hr h(12h-3a) h(6h+3a) h(_6h*6a) “—1212 6.21 8 8 8 8 h(_ar_br) h(_a)’+bY h(aY +87) h(al’ -bi’) _hl' hr ht’ _hY I2!) 60 6b 60' m) 3.2 121) 3a 8 8 8 8 h(_aY+hY) h(_a)’ _hl’) [1(01" _hl’) h(al' Jr) _hl’ hr 111' _hY K/=‘ 6b 60 121» 60 121) 3a 6b 3a 8 8 8 8 ‘ ’23" ’3." “’1 ”(32:25) "(‘21-le hI-ZZ-f2’.l"(—‘él'*f!;l _hY _hY [11’ I11’ h 01' _bY h a); 121'; h _aY 127 h _aY_bY , 8 8 8 8 (6b 6a) 131) +60) (31211212) (61) 120)? _’,’,Y, _hY hY hl’ _al’_bl") h(_a‘Y+bY) h(aY +177) h(aY _hY) 8 H 8 8 6b [20 3b I20 3!) 6a 6b 6a 1 by I17 I17 ht’ 01' (11’ 01’ [21’ 07 bt’ 01’ by 1 K 8 8 -8 _8 h(-3b+l2a)h(—6b—120)h(6h-6a) h(3b+6a)l (437) 81 4.2 Cubic Spline Interpolation Functions Spline functions [27-38] are mathematical tools able to build a curve constrained to pass smoothly through a set of data points. They exist in various orders, but cubic splines are the most widely used in engineering practice. They are expressed by third order interpolating polynomials that share the characteristic of having continuous first and second order derivatives. It follows that cubic splines do not have the “wiggle” problem associated with high-order interpolating polynomials. In their general form, the cubic spline can interpolate as many data points as needed. However, because our search for the proper penalty parameters involved solving really large system of equations symbolically, we had to limit the number of DOFs by using few pseudo nodes to represent the interface elements. A natural cubic spline interpolation over three points appeared as the best choice for having a limited number of DOFs while maintaining a smooth approximation of interface displacements. The word natural indicates that the second derivatives at the end points are assumed to be zero, any extension of the spline beyond the ends would be straight. The Abaqus user element developed using the natural cubic spline interpolation behaves very well in most cases; it is also very computationally efficient. Nevertheless, for some categories of problems this approach can’t achieve satisfactory results. In these cases, a general form of the cubic spline interpolation is adopted. Second derivatives at the end points are calculated by differentiating twice a cubic function which passes through the first four pseudo-nodes along the interface path and another cubic function that passes through the last four pseudo-nodes along the interface path. The Abaqus user element subroutine associated to general form of the cubic spline, allows the user to specify the desired number of nodes 82 composing the interface element. Both cubic spline formulations are described in the following two sections, starting with the general form. 4.2.1 General Form of the Cubic Spline The mathematical description of the cubic spline formulation can be found in several books [27-3 8]. Anyway, for an efficient implementation in our computer code we found very useful to adopt the matrix approach followed by Jonathan Ransom [40]. In this section, his methodology is reported. Given a series of points x, (i = 0,1,...,n) which are generally not evenly spaced, and the corresponding function values/(x), the cubic spline function denoted g(x) may be written as: __ 8.2-(1‘1) (x141 —x)3 _ gJ’I (x141) (x—x,)3 _ _ g(x) — 6 l: Ax, —Ax,(x,+, x):| + 6 Ax, Ax,(x x,) + (438) +f(x,)|:(_x’_:'A;—l_x_)]+f(x,+, )[Eixéll where Ax, = x”, — x, and g,” denotes differentiation twice with respect to x. This equation provides the interpolating cubics over each interval for i = 0,1,...,n and may be given in matrix form as: g = fig,“ + 1,1 (439) For each of the k values of x at which the spline function is to be evaluated, A x, S x, Sx k = 1,2,..., p, and p is the number of evaluation points. The T, and T, 1+1 ’ matrices may be written in the form: 83 (440) where 0 for x, < x, or x, > x”, “ J 17(x1+l —xk )3 ' ‘ (t, )1.) = g T—Ax,(x,+, -x,,) forx, S x, SxH, and] =1+l (442) 1I(x -2)— g k—AxI—L———Ax,(xk —x,) forx, S x, SxH, andj=i+2 l 0 for x, < x, or x, > x”, (7,) =l(x ”' x,) forx S x, Sx andj=i+l (443) - k.) AI 1 1+1 (xk—x,) forx S ka ,and;=i+2 1 Ax, and g,2:r(x0)3 g,“ x g.“ = : ' ( J) (444) g,xr(xn)) 84 (f(x0)l (f(x‘) ) (445) L‘f(x")2 Note that there are, at most, two nonzero coefficients in each row of the T, and T, matrices given above. Applying additional smoothness conditions (i.e., equating the first and second derivatives of adjacent interpolating cubics at x,) yields a set of simultaneous equations of the form; 1%]85 (IQ-1 ”Fix-x1078“ (x, M11188 (x,+1 ) = I I (446) 6[f(x,..)-f(x,)_f(x,)-f(x,1)] 1:12,- 3., (A702 (w,)(mlfil) If the x, are evenly separated with spacing Ax, then the Eq. (446) becomes: .f(x,..)—2f(x,)+f(x.-.) (4x02 (447) [llgvu (x"1)+[4]gt¥ (xl)+[1]gu (X1411 ) = 6|: Eqs. (446) and (3.4) may be written as: Ag,“ = Pf (448) The coefficients of matrices A and P are dependent upon the end conditions. Whether the equations are of the form of Eq. (446) or Eq. (447), there are n—l (x,,). The two necessary additional equations are obtained by specifying conditions on g,“ (x0) and g,” (x, ). For a natural spline, gr,,(x,,)= g__,,(x,,)=0. In general, these second derivatives are 85 calculated by differentiating (twice) a cubic function which passes through the first four pseudo-nodes along the interface path and another cubic function that passes through the last four pseudo-nodes along the interface path. Evaluating this cubic function: g(xo) = a0 + a,x + a,x2 + a,x3 (449) at the first four points gives: I800) '1 x0 x; 231 g x 1 2 3 4 ( ') )= x' x', x', or Na=g (450) 80%) 1 x2 x5 x2 Ig(x,)) _1 x, x32 x,“ Solving for the coefficients yields a = N "g or '- 1 I: ‘ (451) b From the cubic function, g‘_u(x)= 2a, +6a3x, where a; and a, are determined from Eq. (45]). Equation (451) is valid for evenly spaced as well as arbitrarily spaced points. Similar expressions are obtained for the cubic function passing through the last four points where coefficients of the inverted matrix similar to those in Eq. (451) are denoted H“ for k,l =1,...,4. With these end conditions, the matrices of Eq. (448) are given for equally-spaced points as: 1 0 0 l 4 l 1 4 1 A: '. (452) 1 4 _. O 0 1.4(n+lxn+l) and l-I)1 p2 [73 p4 - 6 12 6 E _Xx_ X; 6 12 6 P: E; “'1; Xx— (453) _ 771 F2 F3 F4 -(11+lxn+l) where P, (x) = 2n“, + 6n,,,x, and P, (x) = 2n,,, +6n4kx, 1+1 for k,l=l,...,4. For unevenly spaced points, the tridiagonal A and P matrices may readily be obtained from Eq. (446). In Eq. (439), the spline function g(x) is expressed in terms of the functional values f(x,) as well as second derivatives of the spline function, gm (x). However, it is desirable to express g(x) in terms of the function values f(x,) only. This manipulation is done by solving for g,“ (x,) in Eq. (448) yielding: g... = A"Pf (454) Substituting in Eq. (439) yields g(x)=T,A"Pf+T,f=(T,A"P+T,)f=Tf (455) 4.2.2 Natural Cubic Spline over Three Points In the literature it is possible to find two kinds of natural cubic splines. One type interpolates the given data using a single function for the entire domain, while another one defines different functions for each interval between two data points. To the first category belongs, for example, the following type of cubic splines interpolations functions: 87 1 1 1 3 1 3 1 3 T.(x)—Z—§y(x)+§|y(x)| ‘Z'ym‘ll +§|y(x)—2| 3 1 3 1 3 1 3 T2(x)='2‘—Zly(xll +§|y(x)—1| "11700—21 (456) 3 1 l 3 1 3 1 3 T,(X)=-Z+§y(X)+§ly(x)l ‘21)”(xl—ll +§|y(x)—2| x with: x =2— y() L We have tested this particular set of functions, built to interpolate three equally spaced data points over an interval of length L, in an interface element composed by three nodes. The displacement field V of the entire interface frame was approximated in terms of the nodal displacements as: V(x) = 'T.(X) 91" +T2(X) (12‘ +T.1(X) q; (457) Though it is very easy to implement into a computer code, this set of cubic splines doesn’t allow an exact numerical integration to be obtained for several intervals of integration. For this reason, in our path for programming an interface element, we decided to replace this type of cubic splines with one of the second type. For a three equally spaced points interface frame, if x is the coordinate along the interface, the interpolation functions over the first interval are: 5(x—x,) +(x—x,)3 T'(x):]_ 4h 4173 _ 3(x—x,) _ (x—x,)3 T2”) _ 211 2113 (458) T.(x)=—(x’x')+(x“’§')‘ - 4h 48- Where x, is the coordinate of the first node and h is the distance between two points. Figure 13 plots these functions over the first interval for the case where: h = l, x, = 0. 88 Figure Similar]: 0.2 0.4 0.6 0.8 1 Figure 13. Representation of the cubic splines over the first interval 1.2 1.4 1.6 1.8 7 2 Figure 14. Representation of the cubic splines over the second interval Similarly, the interpolation functions over the second interval are: _(x-x,)+3(x—x,)2 _(x-x,)3 T = 4x) 2h 4h2 4h3 (47» _ _3(x_le2 (x_x2)3 T,(x)—l 2h2 + 2h3 (460) 2 3 T‘(x)=(x-x,)+3(x—x,) _(x—x,) (46]) 2h 4h2 4h3 89 Where x, is the coordinate of the first node and h is the distance between two points. Figure 14 plots these functions over the second interval for the case where: h=l,x,=0. 4.3 2D Beam with a Vertical Interface — Single Variable We now provide a numerical example of application of the Penalty Frame Method. In order to avoid large matrices, we first show the solution to a single variable 2- D problem. 13 14 8 0 V1 6 1 15 V2 1 12 F 4 2 0 V3 9 10 L=2l3 Figure 15. 2D Cantilever beam under axial load applied at the tip. The extension of a cantilever beam under a uniform load F applied at the tip is 90 analyzed using two groups of bi-linear rectangular elements connected to each other by a three noded interface frame. The configuration of the geometry and the mesh are plotted in Figure 15. The single dependent variable involved in this FE problem is the axial displacement 21. We have shown in Chapter 2 that in the Penalty Frame Method the displacement continuity constraint is imposed through two penalty parameter y, and 7,. Thus, the TPE of the system in this study is given in the following form. 7: = ”01 + 7:92 +%7, “V — u, )2 ds + $7, “V — u,)2 ds (462) S S Taking the first variation of 7: and setting it to zero it follows: = 0 (463) 141.111.114- (IE-(1'2 67: = 67:9, + 872,, + 7, [(821, u, )ds + y, [(81/ V)ds — S S — y, “bu, V)ds — 7, I(6V u, )ds + y, [(6u, u,)ds + (464) S S 8 +y, 1(6V V)ds—y, [(821, V)ds—y, [(8V u,)ds S S S Since there are three elements on the interface of the sub-domain 1, it is necessary to split the related integral in three parts. In the same way, the integral regarding the sub- domain 2 needs to be divided in two parts. For example: 7, I(5u,u,)ds= b 2/3 4/3 2 =y, 1(6u, u,)ds+7, I(6u, u,)ds+7, “524, u,)ds = 0 2/3 4/3 (465) 213 413 =7.1(5q1N.+5q2N2)(q1/V1+q2N2)dS+711(542N1+5613N2)(€12N1+613N2)d5+ o 213 2 +71 1(693N1+5q4N2)(an1+q4N2)dS 4/3 91 I 2 721(5u, u,)ds = 731((5u, u2 )ds + 7,1(5u, u2 )ds = S O l (466) 7 b 1 : 721(6‘11N1 +6q2N2)(qlNl +q,N,)ds+y,J(6q,N, +6‘13]\/2)(‘12Ni +q3NzldS 0 Where in the integral over the sub-domain 1 the DOFs q,, q,, q,, q4 are, respectively, the axial displacements of the nodes: 2, 4, 6, 8. While, in the sub-domain 2 q,, q,, q3 correspond, respectively, to the axial displacements of the nodes: 9, 11, 13. N, and N, are the linear Lagrange interpolation functions. For the present problem, V is approximated on the entire interface frame in terms of the axial nodal displacements v,,v,, v3 of the pseudo nodes V,, V,, V3 as: V=T,v,+T,v,+T3v3 (467) Where 7,, T,, T3 are the natural cubic spline interpolation functions defined over the interval V, -— V, by: ,(y):,_s(y—y,)+(y—y.)— 4h 4h3 My-M) 0~903 T) = — 4 -(y) 2h 2113 ( 68) ___(y—y1) (y_y1)3 rxy)_ 4h + 4h3 and over the interval V, — V3 by: 3 ELW:HIy-yfl+3CV:hY_Iy-yfl 2h 4h2 4113 IKy-hY (y-y-)3 T . =1— + - 469 (y-h) My-xf (y-%Y T = + - -—————L— “y) 2h 4h2 4h3 92 The five elements meshing our model are all square. Then, assuming the Young’s Modulus to be 1, for a single-variable square element the stiffness matrix is equal to: (2 __l _l _l 3 6 3 6 _1 2 _1 _1 6 3 6 3 K9" _1 _1, 2 _1 3 6 3 6 _1 _1 _1 2 K 6 3 6 3 ll 3 boo/can: 5 9 10 I2 11 Ill 12 14 131 (470) (471) We saw previously that the global system of equation of the Penalty Frame Method assumes the following form: ' K,”“ K,“ 0 0 K,“ K," + 0," -0,'~‘ 0 0 —0,"' 0,“ +0,“ —0," 0 0 —0;" K,’ +0; _ 0 0 0 K," 0 0 0 K ,” K00 2 .2 — l qz 0 LqZ . I f,” f.’ 0 f2. f2 \ ) (472) J Where [K ,1 is the stiffness matrix for the domain on the left three elements, while [K ,] is the stiffness of the two elements on the right. The [G] matrices need to be evaluated by the following integrals. 03‘ = 711(NI'N11dS 9:“ = 8,101 2, )ds 9," 4971” G,“ =7.11(T.,"T1)ds (473) The degrees of freedom q,’ are equal to zero, since points 1, 3, 5 and 7 are clamped. Thus, in order to find axial displacements of the unconstrained nodes, only the [K,”] part of the matrix [K ,] is needed. [K,"] is associated to the DOFs of the points 2, 4, 6 and 8 and can be evaluated to be. (j -2 0 0) _1 2+2 _1 0 K7: 6 3 3 6 0 _1 2+2 __1 6 3 3 I 2 0 0 - K 6 34 (474) The matrix [K,] , related to the axial displacements of the nodes 9, 11, 13, 10, 12 and 14, is computed as: (2 _l 0 _l _l 0) 3 6 6 3 _1 2+2 _1 _1 ___2 _1 6 3 6 3 6 3 0 *2 i 0 -§ *2 K2’ 1 1 2 1 ._ _ 0 __ 0 6 3 3 6 _1 _2 _1 __1 2+2 _1 3 6 3 6 3 3 l I l 2 0 - — 0 - h 3 6 3} (475) In order to provide an example of the integration of the [G] matrixes, the [GI] matrix is calculated in the following lines. The four rows of the matrix, in order from row 1 to row 4, are: f7] x (711 +71] +T[] )a’ 197v] +43v2 13v3 — xv v v x: r — 0 2 ’ ’ 2x 2 3x 3 810 405 810 3 (476) 94 2 f3(3 x)(T1[xlv/ +T3[x]v2 +T3[X]V3)(/x + 0 2 4 4/3 —x [199 v1 1681v2 24] v3 £3 [ 3 2 ](T’[x]v' ”WM +T3mv3m'x : 6480 3240 - 6480 / 3 ’ (477) 4 x _. 3 fj[ 2 3 1mm v, + Tzlx] v.» + 73le v3) dx + 3 a 2 2 —x 7 241V] 168] v2 1199v3 2 *(Tllxl V1 + T2le V2 + T3le V3) {/x = - 4/3 6480 3240 6480 3 (478) 4 2 x—3 13v] 43v2 197143 4 2 (T/[x]v/ +T2[.rIV3 +T3[X]V3)d’x=- 810 + 405 + In our case the matrix [Cf] determines the contribution to the axial displacement of the points 2, 4, 6, 8 by the interface points vl,v2,v3. 4 I97 43 _ I3 \ 810 405 8/0 1199 I68] _ 241 as _ 6480 3240 6480 (’1 "7’ _ 241 168/ 1199 6480 240 6480 __ I3 43 197 l K 810 405 810 (480) The other [G] matrices are calculated to be: 2 I 0 0 ’ 9 0 1 4 1 0 0 1 2 (481) 197 1199 _ 241 _ 13 l 810 6480 6480 8/0 451 _ '. 43 168/ I68] 43 (’1 '7’ 405 240 240 405 _ 13 _ 24/ 1199 197 810 6480 6480 870 (482) 95 G? = 71 SS 2:72 239 840 39 280 4] 840 239 840 ,39 280 _4I 840 39 280 34 35 39 280 39 280 \k“ NNQ 9 v2 4] 840 39 280 239 840 _4/ 840 39 280 239 840 240 10 73 Then, the “stiffness“ matrix associated with the interface element is: ( 271 9 7l 9 0 _ I977] 8/0 _ 437/ 405 I37] 8I0 0 AClue = 0 0 7] 9 4 7l 9 7I 9 0 _ II997I 6480 _ I68] 7] 240 2417/ 6480 0 0 0 0 7l 9 4 7] 9 7I 9 24 I 7] 6480 _ 168/ 7! 240 _ I I 99 7I 6480 0 0 0 _ I977] _ 437/ 8/0 405 0 _ II997I __ [6817] 6480 3240 + 7] 24] 7I _ [6817/ 9 6480 3240 27I I3 7I _ 43 71 9 8]” 405 I37] 239(7I+72) 39(7I+72) 8I0 840 280 _ 437/ 39(7I+72) 34(7I+72) 405 280 35 _I977I _4/(7/+72) 39(7I+72) 8]!) 840 280 0 _ 73 72 _ 972 240 40 0 _ 72 _ 472 [0 5 772 _ 972 240 40 I37] 8]!) 24I7I 6480 _ II997I 6480 _ I977] 810 -4I(7I+72) 840 39( 7] + 72) 280 23 9( 7I +72) 840 7 72 240 72 IO _ 73 72 240 0 0 0 0 _ 73 72 240 _ 9 72 40 7 72 240 72 3 72 6 0 - 42 I0 _ 472 5 _ 72 I0 72 6 272 3 72 6 2 3 (483) (484) (485) (486) (487) 0 I) 7 72 240 9 72 40 73 72 240 0 72 6 2 + 7 3 I (488) And the global system of equations associated to the problem in this study becomes: 96 I 2 4 2’" - ’ + "I 0 0 _ ’92” - 43 7’ ’3” 0 0 0 0 0 0 \ 3 9 6 9 HI” 405 XI!) 9 3 9 6 9 6480 3240 6480 0 _ I + 7I J + 47I _ I + 7/ 24/ 7I _ I68] 71 _ “997/ 0 U 0 0 0 0 6 9 3 9 6 9 6480 3240 6480 0 0 _ ’ 4 P” 2 4 2 4" ’37] _ 43’” _ "2 7’ 0 0 0 0 0 0 6 9 3 9 (VII) 405 HIU _I977I _II997I 24I7I I37I 23‘)(7I+)2) 39(7I+72) _4I(7I+)‘2) _7392 72 772 I 0 0 XII) 648'!) 6480 NH) 8’41) 28!) 8'40 240 III 240 _437/ _ I68] 7/ _ I68’I7I _43)I 39(7I+72) 34(7/ 472) 39(7I+)2) _972 _472 _9)2 0 0 0 405 3240 3240 405 280 35 280 40 5 40 K _ I3 7I 24/ 7I _ II‘l‘l7I _ I977I _ 4Il7I +72) 39(7I +72) 239(7I +72) 772 _ 72 _ 73 72 0 0 - XIII 644V) 64W) XIII H40 280 N40 240 II) 240 0 0 0 0 - 2342 _ 9’2 2’2 2 + 72 _ ’ 4 ’2 0 _ I _ 2 0 240 40 240 3 3 6 6 6 3 0 (I (I 0 _2‘2 _42"? _2’2 _I+ 72 4 +272 _I 72 _I _2 _I 10 5 I0 6 6 3 3 6 6 3 6 3 0 0 0 0 292 _ 9’2 — 234") 0 - ’ 4 72 2 4 72 0 — ’ J 240 40 240 6 6 3 3 3 6 0 0 0 0 0 0 0 — ’ — I 0 2 _ I 0 6 3 3 6 0 0 0 0 0 0 0 _l _ 2 -I -1 2 2 _l 3 6 3 6 3 3 6 I) 0 (I (I (I (I (l (I -— I ~ I (I - I 2 t 3 6 6 3 I A force vector able to reproduce the ax1al load 15 applied as: . q q q ,mw=M000000000,, } 4 2 - (490) Finally, if we choose the penalty parameters to be 7, = 72 = 1000 and q = 1, the system can be solved to obtain: v ._ " 2’ ‘ 6000 “99 ”9 13-3000 (491) u u u — Here can be noticed that the error in the displacements due to the use of the penalty frame method is on the order of magnitude of —1-. 7 97 4.4 Clamped Rectangular Plane Stress Element The previous analysis, even if very useful to explain the method, does not produce any information on the possible relations among material and geometrical properties of the element and the optimal penalty parameter. That’s why the next simple problem is solved through a symbolic calculation. r/wzwx/a 679 ‘N V9 V3 6////”44 11:8 Figure 16. Single element beam connected to a fixed penalty frame Consider a rectangular plane stress element clamped to a fixed penalty frame along one side and loaded on the other. The geometry of the problem is shown in Figure 16. We saw in Section 4.l that the stiffness matrix of the rectangular element is: (h‘ul"bl') h(uY Vhl') h(7ul'_hl') h(_a)'.hYJ hY _hiY _Ill" by \ 6h 30 III? 30 III) flu 6h 61: 8 8' )s' 8 ‘h(n)' _hY) h(uY +h7) ['(44IYth) h(_aY _hl'] IIY _hl’ _h): hY l [2b 30 (h 34) 6b 61) I26 60 8 8 8 I3 h(‘uX_hY) h(qu+hY) h(ul'+hl') h(uY _h)) _h)’ hY hY ~IIY L”) 60 6h 6u (III 311 I2h 3a 15' X X N Ih_uY+hY)h(_aY~bY) h(u}’ 7hr) Margy) _hY hY h)’ _hY . 6b 6d 12h 6a I2h 30 6h 30 N 8 H 8 “’2‘ h hY 7hr _hr ”(II‘IIl h(0Y_b)’) h(-a)’_bY)h(_11)"hYJ R 8 21’ R 3b 6a 6b 60 6h I2” 36 I261 4:: -:"i 41:42” 41:41) 44;: 43) 4 a a .u .u -hl’ _h)’ hY hY h(_a)'_bl')h(7aY+hY) h(al'*hl’) hral‘_bY) 8 K )9‘ N 6h [2:] 3h I20 3b flu 6h 60 hY hY hY h); _aY M’ _ul’_hY u)’_h)’ a)’ h)’ l 8 8 '8 -8 Il 3b‘I2UI Il 66 I241) I(6b 6(1) h[_a‘b+6u)I (492) 98 In order to separate the DOFs on the interface from those outside the interface, the stiffness matrix is rearranged. To make this change clear, the first row and the last column of the new matrix show the associated DOFs. ( fl ull I0 VI] 0410 ul.’ v10 v12 \ h(al' +hl’) h(_ul' +bl’) h)' h)’ h(a)' _bl') h(_ul’ _b)’) _h)’ _hl' ‘0 6b 3a 6b 6a 8 8 lb 31 I2h 6a 8 8 h(_ul' 1J77") h(ul’ *hl') _hl' _hl’ h(_ul' _hl') h(al' _b)’) h)' h)’ all 6b 6a 6b 3a 3 8 I2!) 6u I2h 3a 8 8 II _h)’ h(ul" +hl" h(__al'+hl') hl’ _hlf h‘al'_bl') h(_u)’_h)’) ”9‘ 8 3 M (So 36 [20 8 8 (h 60 db [’0 h)’ _hl' h(_ul'*hl’) h(u)' +bl') hl' _h)’ h(_ul'_hl’) h(a)'_hl') v" Kc]: N H ‘ 3h [Ba 36 60 N 8 6h [’0 6b 6a h('ul' _hl') h(_u)'_hl') h)’ h)' [1("2. +hl') h(_uY*hl') _hl' _hl' all) 12h 3a [Db 6a 8 8 6b 3a 6b 60 R 8 ;h(_ul'_bl') h(ul' _hlfi'.) _h)’ _hl' h(_al'+hY) h(al'+hl') hl' h)’ “’2‘ I2h 6a I2I> 3a 8 8 6h (Ia 6h 30 R 8 _hl' h)’ h(ul'_hl') h(_u)'_hl’) _hl' hl' h(a)'*hl') h(_a)'+hl') ”0 8 8 6h 60‘ 6h 1’ 8 8 317 6a 3b 1‘ _hl’ hl' h(_ul'_h)') h(ul'__bl') _h)’ h)’ h(_ul’+h)') h(u)'*bl') ”2 l 8 8 66 I20 66 6a 8 8 31> I20 36 6 I (493) The [G] matrices related to the problem in this study can be easily computed to be: 6 1 2 (494) 73 54 —7 1s _ 240 240 240 G2 ‘7’ 4 12 -1 30 30 30 (495) 73 2 240 15 s! _ 9 2 2 ‘7’ 40 5 _ 7 _ 1 240 30 (496) 47 51 _ 41 168 280 1680 435 _ 5I I7 _ 3 l (’2 ‘7’ 280 35 70 _ 41, _ 3 1 1680 70 210 (497) Thus, the global system of equations associated with this analysis is reported below. As before, the first row and the last column of the new matrix show the DOFs (for clarity). 99 cl 0) ”at H u d l d I I ‘(01 II, 9" .' II I. I .'_.' I"..' .‘ II ‘0 I I) I I, I ”e! 0" "I In .4 a! I! ”I In .II I II I d! I I I 00 II. .I‘ II I. I u H. D! II I. I _I! I I! I .. OI III III Io u I} I IlfI IOI I! I I, I 99 ell d. 'nv we . .‘O I' 9.4 u I a Q 'n d I :a- no I, I! ..a’ o o ”7:. A! u N u I a 1:. ”.1.“ 6" .(.¢' 07' on u H a. I H 4:. e a ‘l .1 “IO.” .10] ”r 0.4 .1 n 4:. o u 1. l t I, II II I n 0 at n u 0! - .1 a a u 0! I! cl I! I! 41”., l oi - . l a II he I a] I! a! I) I! I'u' I)" .‘u c I ' o I ‘2 cl. cl) 0 I o 'l I I O o.‘ c e a .0 Ol 0! or. u u ' I}. I0. I I I?) .‘o' II] I! I! ‘I' I: II In I I _I' .I.' I], .l'.’ _ I! 0 I II I. I If: I II (:0 II I. 0! I? I! I! .I II I. I I * .!0' I), II I! ‘) II ‘0 I I I: .‘zi':l’ .‘_ol .II' '. a II I). I! I! II (I I! I- a I o c I n 120’ In on, ‘ (498) Taking advantage of the fact that the interface points v,,v,,v, are fixed, the global sfimleaematrix reduces to: 4.l-r ‘9; too! ”(it .gt‘oy ‘u n 1 u u o .l.9!,!_'l,91'. fielding! ‘u u' c ‘n n’ I y a] a a 2!. i! K: c a .1333!) 1.1.3131! . Mn 3.1 ‘1» u) .I-.o_L-9.u Jet. 93.) I In “I In. 3.1 _u' n a a -9' U -l a a 4.4.1 Axial Load 6: n 1.: u) . g 4 - a a ‘In J" n n l u on ' "” ‘i"" ‘ ' 7 a n In a. Iguana: “34.94332! 9.! ‘n u’ 1 ‘n u.’ a a ‘11," L11 “21.9),2}. U ‘n n.’ 6 ‘n u’ a a u 2.1 “:3th a a In in) -41 -41. .1,-J_ , 9.1.! a a l u: “I .“.’. -91‘ .‘-£1_LIW' _‘I In 6.7 I u "0’1 ' a .I_Q'-u‘ : ,..(.g! __UI _I_'_ l n I)" n u' a (__ll .II\ _II .II ‘ ‘ no a.’ a a (047 _n) g! u ‘1» n’ a a -n‘ .111 .bjl .I_gi _ on c ‘00 u’ l at 1:" _n (at Hi In an .-- -..-...- - . . ._ - a l u ”a, ‘u u’ .lral‘ul _QM _e! u u) a a .l_[._6_[) I] n In 1.2 a a L" ,Id ,2.” ,1-“ N I a In 00’ l n 1:.l I! (or If‘ (0' Il\ . - I -~' 0“ I a I n (20’ ‘n a." We are interested in solving this problem for two simple load cases: axial and transversal load. According to Figure 17, we consider the axial extension first. Figure 17. TVs iv; bV: ll 12 L————J L=a Single element beam under axial load 100 A force \ Inx‘erting can be ex containinl dllIICUll t. A nodes 9 z direction fOllO\Vlng 4.4.2 T I Figure I A force vector of the following kind is applied: force = (o,o,o,o,§,1:—,o,o} (500) Inverting [K] and multiplying by the force vector, the displacements of the four nodes can be evaluated in symbolic form. However, the results are very complex expressions containing the geometrical, material and penalty parameters. As a consequence, it is difficult to arrive at any conclusions based on those outcomes. A reasonable simplification consists in imposing the vertical displacements of nodes 9 and 11 to be zero. In this case those nodes are directly fixed in the vertical direction without the interposition of the interface frame. Under this assumption the following solution is obtained: u u =—£— 9’ II byl (501) uIO’ “Hz—9.1+: (502) th by, 4.4.2 Transversal Load In this section we consider the transversal load, as plotted in Figure 18. a... 11 12 i ’V1 *V2 I) Q ” 3V3 9 10 v Z3- 2.; 7” - id _ ‘ L—a . N \' Figure 18. Single element beam under a transversal load 101 A force vector of the following kind is applied: force = {0, 0,0,0, 0, 05—3} (503) Inverting [K] and multiplying by the force vector, the displacements of the four nodes can be evaluated in symbolic form. The same simplification applied to the axial case is assumed. This time the horizontal displacements of nodes 9 and 1 l are set to zero. Under this assumption the following solution is obtained: V9.» “1:87, _ 602Q uIO— 2 2 a hY+2b hY '2 ath+2b2hY 40(202 +b2)Q Q v10:vl7— ‘_b(a2 +2b2)hY +32 A numerical test has been performed in order to verify that the Abaqus quadrilateral plane stress element CPS4 behaves according to the finite element model developed above. The simple FE analysis of a single CPS4 element clamped on one side and loaded on the other was executed. Then, results from the penalty frame model were compared with those from conventional Abaqus FE solution. Results are nondimensionalized by the Abaqus FE solution. The penalty parameter is: 7, =1000. ABAQUS Interface Interface with 7, z oo Vertical displacement of the tip 1.0000 1,0001 1.0000 This confirms that our FE model matches the one implemented in Abaqus. 102 4. 5 Penalty Parameter and Element Properties In this section the solutions obtained are investigated in order to find a relation between material and geometrical properties of the beam and the penalty parameter. 4.5.1 Axial Extension The exact displacement of the tip is given by the expression: W = —"—F (505) 766 The penalty parameter solution u”"""’”' differs from the exact one by the presence of an additional term —I:—. b 71 upenu/ty = a + L F (506) b h Y b 7, penalty exact The ratio u over u is evaluated in order to underline the dependences between the two solutions. a 1 penalty + _— F u _ b h Y b 7, mu — = 1+ (507) u a F yl b h Y Now, if the penalty parameter 7, is substituted by: h Y 71 = 13(7) (508) The ratio between the solutions becomes independent of material and geometrical properties of the element. 103 u penalty 1 = 1 + — exact u ,8 (509) The accuracy of the solution depends directly on the value assigned to the parameter ,8. 4.5.2 Vertical Flexure As in the previous section, the obtained flexure solutions are examined in order to find a relation between material and geometrical properties of the element and the penalty parameter. The exact theoretical solution for the vertical displacement of the tip of the element cannot match that from the Finite element formulation. The appropriate term of comparison is the outcome provided by the same element without using of the penalty interface frame. We call this solution w” : _ 4a(2a2 +b2)Q - b(az+2b2)hY I'Ii The ratio w”"’“’"" over w” is computed to be: 4a(2a2 +b2)Q Q WWW b(a2+2b2)hY+—bfl7., _1 (a2+2b2)hY w” = 4a(2a2+b2)Q — +4a(202+b2)Q b(az+2b2)hY =1+ ahY bth + 4(2a2 +b2)Q 2a(2a2 +b2)Q Substituting penalty parameter 7, with: (a2 +2b2)hY 40(202 + b2) I’lzfl' The previous ratio becomes dependent only on ,6. 104 (510) (511) (512) penalty 1 w” = 1 +1-3— (513) Thus, as stated before, the accuracy of the solution depends directly on the value assigned to the parameter ,6. 4.6 Building an interface Element Using Abaqus UEL So far, the way an interface element stiffness matrix is defined has been explained in detail. The evaluation of the matrices [G] has been widely described too. The stiffness matrix and generalized vector of unknown displacements associated with the interface element are: 6." —G.'-‘ 0 4' _G,“ G,“"+G‘2"" —G§’ q,, (514) o —G;‘ G; q; We now have the complete set of information for building an Abaqus User Element Subroutine, UEL [39], able to connect different meshes of 2-D quadrilateral elements. The subroutine receives all the necessary information about geometry and material properties of the two connected meshes from the input file. The parameters provided are: the nodes on the interface, thickness and Young’s Modulus of the two domains. Two version of the computer code for 2-D quadrilateral elements has been developed. One release has three equally spaced nodes along the interface whose displacement field is interpolated by natural splines. It is computationally very efficient, but it doesn’t guarantee accurate results in some particular configurations. This happens 105 since the continuity of the second derivatives at the connection between different interface elements is not enforced. The version of subroutine based on the general form of the spline interpolation, allows varying the number of points composing the interface element. It successfully connects long boundaries and many finite elements with only one interface element. It does always guarantee accurate results, but as the number of points grows the computational efficiency decreases very fast. For example, an interface element with 100 nodes has a stiffness matrix of 200x200, since the DOFs are two. The presence of one element with a big stiffness matrix, makes the global stiffness matrix less banded and strongly affects the time required by the solver. The UEL subroutine is written in Fortran 77, according to the requirements stated in Abaqus Manuals [39]. All the variables are defined to be in double precision format, in order to minimize numerical round-off errors as much as possible. In appendix A, a complete user manual for the interface element based on general form of the spline interpolation, is reported. 4.6.1 Automatic Choice of the Optimal Penalty Parameter One of the objectives of this study is to establish an automatic choice of the optimal penalty parameter. The investigations in section 4.5 constitute the basis for pursuing this goal. However, it should be underlined that we don’t need to get an exact value; we just need to know the right order of magnitude of the penalty parameter. In fact, even in the most complex FE analysis, there exists a range of values for this parameter for which the numerical outcomes change very little. This range can be equal to more than 12 orders of magnitude for simple analyses, but usually it doesn’t reduce to 106 less than 2 orders of magnitude in any situations. Taking into consideration the previous comments, we now recall the main results found in section 4.5. We saw that under conditions of axial load, the following expression of the penalty parameter would make the accuracy of the solution to depend only on the value assigned to the parameter ,6 : 7. 4(3) <51» a This expression is very simple, but still contains information that would be difficult to obtain, namely: the length of the element, a. Inside the UEL subroutine only the coordinates of the nodes on the interface are known. Moreover, along the interface the elements could have different lengths. Thus, taking advantage of the good practice in FE of using elements with aspect ratio not too distant from the square shape, we choose to approximate the dimension in the direction perpendicular tothe interface with the one parallel to it. According to the adopted symbols, the expression for 7, becomes: 71=fl'[Ib—Y-) (516) This simplification strongly affects the result regarding the flexural load case. The expression for 7, reduces to: 41412223519231} 1442—11 A final step is required in order to obtain a single expression for an automatic evaluation of the penalty parameter. It consists in eliminating the number 4 from the denominator of the last expression for 7,. This simplification is important since the 107 interface could have any orientation, so it is better to not make any distinction between horizontal and vertical displacement DOFs. In conclusion, we choose to make 7, to depend on the following parameter: hY 71—16(7) (518) Where the value of ,6 is usually assigned to be 1000. 4. 7 Numerical Results In this section the developed Abaqus User Element, UEL, for plane stress quadrilaterals, is employed to solve several problems. These FE analyses will test the precision and reliability of the proposed method. 4.7.1 Patch Test A two dimensional rectangular region is assumed to be composed of two domains. The domains are meshed independently and joined by one interface element. The geometrical configuration of the problem is shown in Figure 19. Other properties are: E, 2 E2 =1 and thickness =1. In discretizing the two subdomains, four-noded bilinear plane stress elements are used. Three different conditions are investigated: uniform axial strain in direction 1, uniform axial strain in direction 2 and uniform shear strain. This problem is intended to serve as a patch test for the interface element developed herein. If the continuity constraint along the inclined interface is not correctly imposed by the interface model, then discontinuity of displacements and strains across the interface would result. 108 The condition of a uniform strain in direction 1 is obtained by assigning a displacement equal to 1.6 mm to all nodes on the right vertical side and constraining the opposite side against movements in the same direction. The node at the bottom left corner is constrained against vertical motion. The results are in agreement with the exact solution to the number of significant digits available, i.e. axial strain is uniform in direction 1 and equal to 0.1 in all elements of the patch. The displacement distribution in direction 1 is reported in Figure 20. Similarly, the condition of a uniform strain in direction 2 is obtained by assigning a displacement equal to 1 mm to all nodes on the top horizontal side and constraining the opposite side against movements in direction 2. The node at the bottom left comer is constrained against horizontal motion. 6:21.80 Interface / 1 f / I A=10mm B=16mm L V Figure 19. Patch of incompatible 2D meshes connected by one interface element. 109 .SUUE-37 .2313-01 .462E-Ul .6922-01 .923E-01 .1543-01 .335Z-Ul ,6152-01 ,BISE-Ul .108t000 .231tvUD .354E~00 .477EQUD .EUUEvUU Figure 20. Patch of incompatible 2D meshes - Displacements distribution in the direction I (U 1)- VALUE .SOOI-37 .692t-U2 SIDE-01 BOSE-01 077E-01 3465-01 .615E-01 385E—01 lS4E-01 923E-01 692Z-01 .ASZE—Ol 2315-01 .UUUEoUD Figure 21. Patch of incompatible 2D meshes - Displacements distribution in the direction 2 ( U 2). 110 VALUE .915!-36 .ASDEODD 901:000 351:000 .BDZEoUD 2522.00 103E000 015:001 160:001 305E~01 .450£¢01 .S9St901 741:001 .BBSIOOI VALUE 549E-36 282E900 .563!~00 8455000 .126E000 .IDBI900 6895.00 971E000 .D2Slofll .153l~01 .2822401 .410E601 538E401 .666E901 (b) Figure 22. Patch of incompatible 2D meshes - Displacemenm distribution in the direction I (a) and in direction 2 (b) under uniform shear deformation. 111 The interface element results are in agreement with the exact solution to the number of significant digits available, i.e. axial strain is uniform in direction 2 and equal to 0.1 in all elements of the patch. The displacement distribution in direction 2 is reported in Figure 21. The condition of a uniform inplane shear strain is obtained by assigning displacements to the nodes on the boundaries such that horizontal sides form an angle of 30° with the axis 1 and vertical sides form an angle of 30° with the axis 2. The interface element results are in agreement with the exact solution to the number of significant digits available, i.e. shear strain is uniform in all elements of the patch. The displacement distributions in direction 1 and 2 are reported in Figures 22(a) and 22(b). 4.7.2 Isotropic 2D Cantilever Beam with a Vertical Interface A two-dimensional cantilever beam is meshed as composed by two domains. The domains are meshed independently and joined by one interface element. The geometrical and load configuration of the problem are represented in Figure 19. Others properties are: P=1, E, =E2 =1 and thickness=1. v1 P V: L3=1 V3 L,=2 L,=2 Figure 23. 2D Cantilever beam with vertical interface 112 Only one Interface element is adopted in this simple case, since more would not improve the solution accuracy. This element is identified by the vertical interface made of the three nodes V1, V2 and V3. For discretizing the two domains, four-noded bilinear elements are used. Two different load conditions are investigated: axial and bending loads. This problem may be considered another form patch test for the element developed herein. 4. 7.2.1 Axial Load Under the condition of a uniform axial load applied at the free end of the beam, the interface element results are in complete agreement with the exact solution to the number of significant digits available. Results from two different beam discretizations are reported in Table 7. Table 7. Tip axial displacement numerical results for beam extension. Mesh Solution A baqus 20X8 4.000 Abaqus with interface Element 10x8 Lefi/ 1010 Right 4-000 10x8 Lefl/ 10x4 Right 4.000 Classical - 4-000 It is important to study the stability of the solution obtained by the new method for various values of the nondimensionalized penalty parameter 7 / E . To run this test the automatic choice of the optimal parameter 7 , described in section 4.7, has been disabled. Results are produced in tabular and graphical form in Table 8 and Figure 24. It can be noticed that the solution is stable for a wide range of values of the nondimensionalized penalty parameter 7 / E . 113 Table 8. 2D Cantilever Beam under axial load — Mesh 10x8-10x4 Variation of the solution accuracy with the decimal LOG of the penalty parameter 012 3 4 5 6 7 8 9101112131415 6.00 4.20 4.02 4.00 4.00 4.00 4.00 4.00 4.00 4.00 4.00 4.00 4.00 4.02 4.13 7.42 tio Penalty/Exact 012 3 4 5 6 7 8 9101112131415 1.50 1.05 1.01 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.01 1.03 1.86 Uinterf./U Classical 2.0 41‘ 1.50—- - - 1o’eooeeoeeeoeo’ Luge/E) 0.5 T I r I I I r T I I I I V 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 , Figure 24. Stability of the solution for different values of 'y - axial load. 4. 7.2.2 Bending Load For the bending load case, the comparison is made to both the exact solution and a reference finite element solution. The coarse mesh of bilinear quadrilateral elements cannot exactly recover the classical solution. In Table 9 are results obtained from two different beam discretizations. The FE results, marked as Abaqus, are obtained from a traditional analysis using a compatible finite element model of the beam. The stability of the solution obtained is tested also in this load case for various values of the nondimensionalized penalty parameter 7 / E . Again the automatic choice of 114 the optimal parameter 7 was temporarily disabled. Results are produced in graphical form in Figure 25. The method proves very reliable and stable for bending loads. Table 9. Tip deflection numerical results for beam flexure. Mesh Solution Abaqus 20X2 258.800 20x4 260.000 20x8 260.400 Abaqus with interface Element 10x3 Lefi/ 10x?- Right 259-500 10x8 Lefi/ 10x4 Right 261.100 10x8 Lefi/ 10x8 Right 260.400 Classical - 265-5 Table 10. Tip deflection in a 2D beam under bending load — mesh 10x8-10x4 Variation of the solution accuracy with the decimal LOG of the penalty parameter 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 358.2 269.9 261.1 261.1 260.1 260.1 260.1 260.1 260.1 260.1 260.1 260.3 263.2 306.4 481.8 tio Penalty/Classical 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 1.349 1.016 0.983 0.983 0.979 0.979 0.979 0.979 0.979 0.979 0.979 0.980 0.991 1.154 1.814 . Vinterfjv Classical 2 o 0 e 1 :90 e"e* e e e'o'e'e'e 0 Log (WE) 0 I I fir I r r I r T r 7 U 1 o 1 72 3 4 5 776 877 ,9- 1° 3." -, 12‘ #13, 14: Figure 25. Stability of the solution for different values of 7 - bending load 115 4.7.3 Tension-Loaded Plate with a Central Circular Hole A plate with a central circular hole is subjected to a uniform load at its edges. This well studied elasticity problem allows one to investigate the capability of the interface elements in performing global/local analysis. Geometrical and load configurations are plotted in Figure 26. The plate height and width are, respectively, h = 100 and b = 200; the radius of the hole is r = 2.5 ; the distance of the interface from the center of the hole is Interface ---—J UN UI p-_-___——_-_ -— b=200 Figure 26. T ension-loaded plate with a central circular hole (not to scale). Taking advantage of symmetry, only one quarter of the plate is modeled. A fine mesh is applied in the area between the hole and the interface, while a much coarser mesh is adopted elsewhere. The complete finite element model is reported in Figure 27a. For a better view, the area near the hole is depicted in Figure 28a. A conventional FE model has been built with Abaqus in order to compare the maps of displacements and stresses, 116 see Figures 27b and 28b. The elements used to mesh the models are plane stress bilinear quadrilaterals, Abaqus CPS4. Five interface elements are employed along each segment of the interface. Horizontal U1 and vertical displacements U2 do not show any change in their values across the interface, as plotted in Figures 2% and 30a. Figures 31a and 32a demonstrate that only light discontinuities in the values of stresses at the interface are present. Results from the Abaqus conventional FE mode] are reported in Figures 29b- 32b. Comparing the maps of the displacements, it can be noticed that the interface elements does not affect at all their distribution. The maps of the stresses show just little differences among the models. It is important to remember that through the penalty formulation continuity of displacements has been enforced, but not continuity of stresses. Therefore, some discontinuities across the interface are expected and do not depend on inaccuracies of the penalty model. The finite element interface technology developed at NASA LaRC [6-11], which uses Lagrange multipliers to enforce the interface constraint conditions exactly, is not able to do any better in preserving continuity of stresses across the interface. 117 L. Figure 27. (a) Interface element and (b) conventional FE mesh. 0)) 118 (b) Figure 28. (a) Interface element and (b) conventional FE mesh - zoom of the deformed configuration. 119 U1 VALUE +6.956E—38 .717E+33 .543E+Ul +2.315E+Dl +3.087E+01 +3.859E+01 +4.6305+Ul +5 .402E+31 +6 .174E+Cl +6.946€+m +7.717E+Dl +8.489E+Bl +9.26lE+Dl +1.CC3E*02 + + P“ \l .717E+01 .489E+Ul .261E+01 .333E+UZ (b) Figure 29. (a) Interface element and (b) conventional FE solution - horizontal displacement distribution (U,). 120 U2 VALUE —1.533E+01 -l.415E+01 -1.297E+01 —1.179E+Ul —1.DGIE+UI ~9.432£+UO —8.253E+DU —7.U74E+OU -5.8958+UD -4.716E+UU —3.537E+UU ~2.3SBE+UU —1.179E+DD +4.160E-38 VALUE .533E+01 .415E+01 .297E+01 .179E+01 .061E+Dl .433E+OU .254E+UU .074E+UU .895E+UU .7IBE+UU .537E*DD .3588+OU .179E+UD .821E—38 (b) Figure 30. (a) Interface element and (b) conventional FE solution - vertical displacement distribution (U1). 121 VALUE .178E-02 .EBBE-Dl .0533-01 .4398-01 .8243'01 .1213+00 .359E+DD .5982+00 .537E+UU .075E+UD .314E+00 .SSZE+DD .7913+00 .UZBE+OD (a) VALUE .193E-02 .67DE—Dl .059E—01 .4488—01 .8378-01 .123Eo00 .3GIE+DU .GOOE+DD .839E+UO .07BE+UO .317E+00 .556E+OU .7BSE+DU 034E+DU (b) (a) Interface element and (b) conventional FE solution - stress distribution in the direction I (0'11). Figure 31. 122 VALUE .UIUE+UU .924E-01 .7523-01 .SBUE-Ul .407E-01 .2353-01 .063E-Dl .BBUE—Dl .1815-02 .5422-02 .626E—Dl .7998-01 .9713-01 .1432-01 VALUE .013E+UU .SSOE—Dl .775E—Di .6019—01 .426E-01 .251E-U1 .0762—01 .9013-01 .2638—02 .486E-02 .623E—01 .7983-01 .9735-01 .1482-01 (b) Figure 32. (a) Interface element and (b) conventional FE solution - stress distribution in the direction 2 (0'22). 123 tension hole. N; line axi 1.0, 0.5 . Figu“ The elasticity solution for an infinite plate with a central hole loaded in uniform tension says that the stress concentration factor, K, is equal to 3.0 at the upper edge of the hole. Moreover, it provides the equation for evaluating the change in the 0'll along the line axis X 2. Using this expression, we computed the stress variation along X 2 for our geometrical configuration and compared this result with our FE models. Figure 33 graphs the Elasticity exact solution against the penalty interface frame and the conventional FE ones. The results are all very near, validating the precision of the developed method. 3.0 . Elasticity Solution -~ Abaqus conventional 0 Interface Element ‘1.o x2 0.5 .» . , ‘ . -- - . - f . -. 7 ,- ,- ,- . 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 22.5 25.0 1 Figure 33. Stress Concentration Factor Along the axis 2. 124 the hole. mesh the at same Figure 3 model. ( interface pr€\'10u: confi (“1 income The problem has also been studied with a different mesh for the area containing the hole. We wanted to prove that in order to get accurate results it is not necessary to mesh the connecting finite element models with some of the nodes at the interface placed at same location. The area near the hole of the new finite element model is shown in Figure 34; the mesh of the second subdomain is unchanged with respect to the previous model. One interface element of the general form is employed along each segment of the interface. Comparison is made to the same reference Abaqus conventional FE model used previously. Horizontal UI and vertical displacements U2 plotted in Figures 35 and 36 confirm the remarkable accuracy of the method even in case of total nodal incompatibility at the interface. Figure 34. Interface element FE mesh #2 - zoom of the deformed configuration. 125 U1 VALUE «+1.94DE/37 ~7.717E+Ufl +1.543E+Dl +2.315E+Dl +3.087E+01 +3.659E+01 +4.63UE+01 +5.402E+Dl +6.174E¢Ul +6.946E+Ul +7.71?E+Ol +8.489E+01 +9.261E+01 +1.UDBE+UZ VALUE .SQBE737 .717E+UG .543E+01 .315E+01 .387E+01 '.859€+01 .630E+01 .432E+Dl .174E+01 .9463+01 .717E+Ul .489E+01 *9.261E+01 +1.0038+02 (b) Figure 35. (a) Interface element and (b) conventional FE solution - horizontal displacement distribution (U1). 126 Fig, UZ VALUE -1.532E+01 -1.41QE+01 ‘--1.297E+01 -1.1798+Ul .061E+Ul .43DE+UD .251E+00 .072E+00 .8945+UD .715E+UU .536E+UD .357E+UU .179E+UU .897E-3B VALUE .533E+01 .415E+01 .297E+01 .179E+01 .061E+01 .433E+DO .254E+00 .D74E+UU .895E+UU .716E+UU .537E+UU .358E+UU .179E+OO .BZiE-BB (b) Figure 36. (a) Interface element and (b) conventional FE solution - vertical displacement distribution (U z)- 127 4.7.4 Composite 2D Cantilever Beam with a Vertical Interface A 2D composite cantilever beam made of five materials is considered. The geometrical and loading configurations are the same as for the cantilever beam analyzed previously (see Figure 37). The entire left domain of the beam is made of one isotropic material having a Young's modulus equal to 2.9E+7 MPa and Poisson’s ratio 0.25. The right domain is composed of four isotropic layers of equal thickness. They have all the same Poisson ratio 0.25, but the Young‘s modulus E varies. Starting from the bottom, the layerwise E takes the values: 1E+6 MPa, 1E+7 MPa, 5E+6 MPa and 1E+8 MPa. A conventional finite element model (10x4-10x4) has been used as a baseline to test the results from the penalty interface method. V1 P V2 Q , .f . _3 L3=1 1V3 L,=2 L,=2 Figure 37. ZD Composite cantilever beam with vertical interface Under these conditions the interface is expected to undergo abrupt changes of slope. An important property of the interface elements is the automatic choice of the penalty parameters. In fact, the interface elements can join the four different layers to the left domain with the dissimilar stiffness. This behavior can increase the accuracy of the results, avoiding the possible corrupting effect of an overestimated penalty parameter. 128 The results are shown for two possible interface discretizations: 1 interface element and 4 interface elements. Axial and vertical displacements at the free end of the beam under axial load are plotted in Figures 38 and 39. Similarly, deflections for transversal load case displacements at the free end of the beam are reported in Figures 40 and 41 . Displacements from the interface element model compare very well to those from the conventional finite element model. Only in the transversal load case, a single interface element appears to slightly overestimate the vertical displacements. 129 1.0 , y —Abaqus Conventional 0.9 3 0.8 . - - - 4 Interface Elements 0.7 o 1 Interface Element 0.6 . I 0.5 0.4 - 0.3 0.2 . 0.1 0.0 - - -» - 0.E+00 4.E-03 6.E-03 8.E-03 1 .E-02 1 7 2.E-03 Figure 38. Beam extension along the thickness at the free end under axial load. 1 0 W . ....... 7 7 7 g ‘ y j -——Abaqus Conventional * . — — — 4 Interface Elements 0.7 - O 1 Interface Element 0.6 . 0.5 0.4 . 0.3 V 8.0E-03 8.5E-03 l l 6.5E-03 7.0E-03 7.5E-03 Figure 39. Beam vertical deflection along the thickness at the free end under axial load. 130 1.0 » 0 9 ’ ———Abaqus Conventional 08 . ‘ - - - 4 Interface Elements 1 * c 1lnterface Element 1 , 0.7 . 0.6 0.5 j 0.4 i 0.3 0.2 - 0.1 u 0.0 . »— »» -- . - - - — » -2.E-02 -1.E-02 0.E+00 1.E-02 2.E-02 3.E-02 4.E-02 A4_g _.__--#A A-._~_...._ ifizl Figure 40. Beam extension along the thickness at the free end under transversal load. 10 7 ’7 7 7’ ’7 ' 71.; _ __W; _' yog I ‘ —-Abaqus Conventional ‘ I - - - 4 Interface Elements 0.8: . w 0.7 l ----- 1 Interface Element 0.6 - 0.5 0.4 . 0.3 0.2 0.1 . \ 5 O Q § 00 L. A . . ,, . ._ .. .7 -77-. - i V.\ 9.. _L. W, 7. .7 . .7, W. 2 '- 1.15E-01 1.17E-01 1.19E-01 1.21 E-01 1.23E-01 1.25E-01 ___.. __~_ .2- - 2.... ___._ ._~_ ___. ____ ___. __ . _._J Figure 41. Beam vertical deflection along the thickness at the free end under transversal load. 131 4.7. 171110 Hi [h for Sam 1631 del. 4.7.5 Clamped-Clamped Asymmetric Beam In this section an asymmetric beam, clamped at both ends, is loaded with two concentrated loads applied at two different points of the mesh. This problem was studied in [4,13] to validate the proposed interface method. The geometry of the problem and the position of the two concentrated loads are shown in Figure 42. (0,1) (1,1) 0,0 ( ) a (1,0) b Figure 42. Clamped-clamped two-piece asymmetric beam. The data in parentheses represents the x and y coordinates of the points indicated The beam is meshed using two domains with different material properties. The ratio of the Young’s modulus in the domains a and b is 0.2. The domain a is discretized with an 8 by 8 mesh; two different meshes are adopted for the domain b. The mesh used for the conventional FE analysis is set to be 8 by 8, so that at the interface there is the same number of nodes. A different mesh is needed for discretizing domain b in order to test the interface model’s capabilities. An 8 by 12 mesh has been chosen. The two forces, having the same magnitude F=1 but different directions, produce a complex state of deformation with strong stress concentrations at the interface that can challenge the 132 interface element. Figures 43 and 44 compare the deformed configuration, produced by a conventional FEM, with that obtained using a FEM with four interface elements at the interface. No visible differences can be observed. The test is completed by a comparison between the in-plane displacements u and v along the interface as shown in Figures 45 and 46. The four interface elements connect two interfaces with a different number of nodes. Both displacement components are in satisfactory agreement with the conventional analysis. 111 Figure 43. Deformed configuration obtained through conventional Abaqus FE solution. *\ .2 l 17 K , ,- . .. ___ \ S 1‘“ M \i , 1 Figure 44. Deformed configuration obtained using interface elements approach. 133 1.0 0.9 I Interfacea 0.8 ~ 0 Interfaceb 0.7 f ,, - f g g 0.6 0.5 . 0.4 0.3 i 0.2 . 0.1 . 0.0-, _ w - -3.E-07 -2.E-07 -1.E-07 0.E+00 1.E-07 2.5-07 3.E-07 4.E-07 5.E-07 6.E-07? Figure 45. Axial displacements along the thickness at the interface —Abaqus Conventional I lnterfacea 0 Interfaceb 0.7 ** ” r r * ’ 0.5 7 0.4 ‘ 0.3 . 0.2 - 0.1 5 0.0 - ,_ ,, -5.2E-07 -4.8E-07 -4.4E-07 Figure 46. T ransversal displacements along the thickness at the interface 134 CHAPTER 5 PLATES 5.1 Finite Element Model In this chapter an interface element for plate elements is developed. A definition of the finite element model associated with the plate element is provided in this section. Plate elements based on the First Order Shear Deformation theory, or Mindlin plate theory has been chosen for our research. The degrees of freedom associated with this element are five for each node: u, v, w, 6,, (9),. However, the equations governing the in-plane displacements u and v are uncoupled from those governing the bending deflections w, 61,, 9,. It follows that all the results obtained from the study on the plane stress quadrilateral are valid for the in- plane displacements of the plate. Thus, we will focus our attention just on the bending DOFs w, 6,, 61,. The First Order Shear Deformation Theory (FSDT) is based on the displacement field of Mindlin: u=uO—26y v=v0+26lx (519) w = w0 The associated bending strains are: 89. 8, = —z—"- (520) ' 6x 135 8,22 5y 6‘ =0 rr=Q+Qfi 5y yr _—6j'+_ (521) (522) 6%) (an 6%) For a single isotropic plate with Young’s modulus Y, Poisson ratio v and thickness t; the laminate stiffness is defined by the following relations. Where: A3=th A: Yt2 l—v 136 VA 0 O 0 0 0 0 0 0 0 A(1‘V] () 0 0 2 0 D VD O 0 m) D 0 0 0 0 DC-V 2 awW {Qt-r ° Q, 0 A, 6. + 93 @r (526) (an (an (an 12(1—v2) 33% 6x gg_§"2 ' 5y 73V=fllo_,% ._ 0y ax ‘ _ 66V Zx— ax .299; I" 6x _ 66x 66) le,t'* ax 6y II MEN. 1=I a, = 2109,), NV ,- 9V = 21(0V ). N ,- We obtain the finite element model for bending. Where: 1161 M V,,]- M’ [K221 W] _[K'il [Kzilr [K331 137 _l 735:8” V‘. 6N 6N. KJ' =th{ ([aN’ -’ + 6N, ’ ' o, 5" 6x 6y 53’ l Assuming finite element interpolation of w, 6,, and 6V in the form: 5. w} (530) (531) (532) (533) (534) K32 = th{ [[gAyf—Lijdxdy} (535) K53: _th{ I[%N,]m@} (536) 3 ,aN 3 ,6N K52 = J. 24Yt 6N, j + Yt 2 6N, .l +th(Nle) dxdy (537) 0. (1+v) 6x 6x 12(1—v) 6y 6y , 3 .aN. 3 0N K0,?! m 2 0N, ,+[ Yt 101v, , dxdy Q“ 12(1—v) 0y ax 24(1+v) ax 3y (538) 3 5N 3 6N. K33: ( Y’ 2 5N3 '+[ Y’ )6”! ’+(th)N,NV dxdy (539) ' 0V 12(1—v) 0x 5x 24(1+v) ay ay Now, the previous defined stiffness matrix is computed for a rectangular plate element having the in plane dimension shown in Figure 47 and a thickness equal to t. A < v L Figure 47. Geometrical dimensions of the plate element In the following lines the [K] matrices composing the stiffness matrix are evaluated. 138 t k(_h~’+L2)zY 6hL(l+v) 12h L (1+y) k(h2 +1.2)17Y K11 k(—2h2+L2)t Y l2hL(l+w k(h2-2L2)IY t 12hL(l+y) K77 : 'tY (-212 12.112 (41:12.8) (-1.v)) 72hL(-1w) (1w) tY (1.212.112 (2111.2- :2) (-1.y)) 72hL(-14v) (14v) tr (21.212412 (221.2 :2) ( -1.v)) 144bL(—147) (14v) tY (4 L2 c2412 (41:12.3) ( -1w)) c y (L2 1:2-b2 (2:12-22) (-1.y)) 72111. (-1.y) (1w) r (412 c2412 (4212.12), (-194) 72bL(-1ov) (10V) t? (4 L2 c2412 (421.2“:2) (—1.v)) 144bL(-1¢v) (14v) tr (2 L2 c2412 (2kL2-t2) (-1.v)_) (ct—232 +19)” 10323332)” (of—23% ‘ 12hl.(l+v) I2hL(l+v) 12hL(l+v) 1(112 +1.2)! Y k(h2—2L2)t Y __1 (112.42): Y 6hL(l+v) 12hL(l+y) I2hL(1+y) ’3 ('32 7242 )3. Y k ('12 +32 )3 Y (“23132 )3 Y, l2hL(l+v) 6hL(l+v) 12hL(l+v) _ k ([12 +L2) (Y k(-2h2+L2)t Y k ((12 +1.2)th 12hL(I+v) l2hL(l+y) 6hL(I+y) 1 (540) tr (21.2 :24? (211.2 1.2) (—1w)) :1! (41.2 84.2 (4112.8) (4.») t 14451. (-1.v) tY (41.2 1:241? (411.2 144111. (~1w) t? (L2 c2-h2(21u.2- err (-212 c2412 (41:12.12) ( -1.y)) 72hL(-14v) (1w) 1441114 -1.y) (1w) 14411144") (1») ’ 721:1.(-1.v) (_ kLrY _ kLtY _ kLrX 12(1 +10 24 NH) 24 (I +r') _kl.tY _kLtY _kLtY 24(l+v) 12(l+v) [2(I+V) K12: kLzY kLzY kLzY 24(I+V) 12(l+v) 12(I+V) kLtY kLtY kLrY K l2(l+y) 24(/+V) 24(l+y) K33 = 'tY (1:2 (-2 22.4111? (—1.y)).L2c2 (-1.v)) 72hL( —lov) (14v) cr (4122 (1.2.21.2( -1.v)).1.2 12(- -1.y)) 144111. (-1¢v) (14v) :1! (2122 (c2okL2( (-1w)) -L2 1.2 (-1.y)) 14415L(—1.y) (1w) .2kL2( 1-.y)).1.2:2( -1.y)) 72hL( ~10V) (14?) cr (4122 (12.11.2( 1.y)).1.2 c2(- -1.y)) 1.44111. (-1¢V) (10v) tr (122(—2c2.411.2(-1.y)) 1.2 c2(-1w)) 72hL (-lov) (141!) _cr (112 (:2-2kL2 (4.13)) .L2 :2 (—1.y)) 72hL( -1¢V) (In!) car (2h2 (12.111.201.10) -L2 c2 (4.11)) _cr(h2(c2 r (41:2 (c240: L2 (-1.y)) .12 :2 (-14v)) cr (1.2 (.2 8.41:1? (-1.y)).1.2_ :2 (_-_1.y)) CY (112 (CZ-2kL2( -14Y))9L2t2( -1¢v)) c r (21:2 (1.2.): L2( (-1.y)) -L2 c2 ( 41.12)) " 144 11 L(-1:ov) (1w) 72hL(-lov) (1w) tr (1:2 72hL (~10Y) (107) 1“hL( -lov) (147) tY (2122 (c2.kL2(- 1.13))-12 1.2( 4.») 1441213 (~10V) (10v) tY (152(c2-21L2(-1:v)) L2c2(-1.y)) 72L1.(-1.y) (1w) CY (212 (-2£2¢43L2 (-10Y))ol12 62(-1ev)) 72514-1.» (1w) cr (4122 (8.111.“ -1.y)).1.2 12( -1.v)) 1“ b L (-14v) (14v) ( K13 K hth l2(l+v) _ hth 12(l+y) _ hkrY 24 (I H) hkrY 24 (1 +10 1“ h L (-1«rv) (14v) hkrY 12(l+v) _ hth [2(l+r') _ hth 24(l+v) hth 24 (l+r') 139 hth 24(I+V) _.hk!Y. 24(l+v) _thrf 12(l+v) hthV l2(l+y) 144bL( —1ov) (14v) tY (21.2 1.2.112 (21: 1.21:2) (4.") 144bL(-19v) (19v) tr (L2c2-112 (21:1.2-12) (-1450) 72bL( -1¢v) (14v) 1:? (-21.2 1:24:12 (41:12.12) (4.53)) 72151. (-1¢v) (1w) 1 (541) (Ivy) otz) (-19V)) (14v) :2) ( —1¢v)) (14y) _ kLtY ) 12(l+y) _ kLtYfi 24(I+V) 13L}! 24(l+v) kLtY‘ 12(I+V) ) (542) t Y (2122 (c24kL2(-1ov)) -L2 c2( -1.v)) 144111. (-1w) (143!) __c r (112 (c_2_-2 151.2( (3.51) 4.23.2 (71.5)) 725L(-1w) (1w) (:2 3’13 323-133)) 3,32 321-13») 721214-1410 (1w) ' " cr (41:2 (1.2.11.2 (-1.y)) .L2 :2 (-1.y)) 1“bL(-1+V) (10V) . . . _cr(112(c2-2kL2(-1w)):1.2c2(—1.y)) l 72hL(-1oy) (1w) cr (2112 (8.2L2( 1w))-L2 t2( 1.y)) 14421144451) (14v) cr (41.2 (12.12 (-1.y)).L2c2 (-1w)) 144hL (-10v) (14V) tr (1:2 (-2 c2.4kL2(-1w)) 4L2t2 (1.y)) 72hL(-1.v) (1w) 1 (543) _hth. ) 24(1+v) _ hktj- 24(I+V) _ hth 12(I+y) hthV /2(/+V) l (544) t 131’ t3Y(I—3y) 131’ :3Y(-1+3y)) 96(—I+r) 96 (4.42) 96—96 V 96 (-r+y2) (3 Y_(-I+3 y) :3 Y 13 Y (1—3 9) t3 y _ 96 (4+9?) 96—96 V 96 (4+?) 96(-1+v) K23 = 3 3 ,3 y r Y(—l+3 v) r3 Y r Y(l—3 9) 96—96 ,, 96 (4+9?) 96 (—1+y) 96 (4+9?) I3Y_(1-33’l t3)? r3Y(—I+3y) 13,15 _. ( 96(—1+;3) 96(—l+v) 96(—l +132) 96—96 V ) (545) 5. 2 Clamped Rectangular Plate Element Following the same path as in the Plane Elasticity case, the solution for a rectangular plate element has been developed. Figures 48 and 49 show loading and geometrical configuration. Analyzing a plate, nevertheless, requires the presence of 3 degrees of freedom for node. Having a total of 7 nodes, we should report a matrix with 21 rows and 21 columns, an effort that would be of little help in increasing our knowledge about the method. This since we would show nothing else that adding the same [G] matrices to the global stiffness matrix of the system three times, one for each DOF. In fact, it should be clear that the {G} matrixes depend only on the integration interval and on the approximation functions adopted. If all the DOFs are approximated by the linear Lagrange approximation functions and we adopt the same cubic splines, as before, all the [G] matrices are unchanged. In the following sections the results from our numerical study are reported. Two different penalty parameters are used: 7.. is associated to the transverse DOF w, instead 79 is related to the rotational DOF. 140 V1 1 1 1 2 + V2 b h i 9 1 0 V3 i 3 4 > Figure 48. Single plate element connected to a fixed penalty frame — T op view , la. Figure 49. Single plate element connected to a fixed penalty frame — Lateral view .< n “F‘— 5.2.1 Transversal Load We are interested in solving this problem for two simple load cases: transversal load and a bending moment applied at the free edge. We consider the transversal load first, as shown in Figure 50. The symbolic results obtained from the solution of the global system of equations associated with the problem in this study are quite complex expressions. However, keeping in mind that we only want to know the right order of magnitude of the penalty parameter, some reasonable simplifications can be applied. In this case, groups of terms, 141 whose contribution to the solution is small, have been eliminated. ii V Figure 50. Single plate element under transversal load applied at the free side gyms-f :_ _- __ .'-._‘-_'-_ i The solution obtained under these assumptions is proportional to that found in our study on the Timoshenko beam element. We should not be surprised by this achievement, since the Mindlin plate theory is an extension of the Timoshenko beam theory. F .- 3 2 3 2 ”3:“; L3 + L +—1—+‘L— Q= [LP—L ]+i+£— Q (546) i [h k}: h 7.... 79 3Y1 kCJA 7w 70 31’ ~—— (1 ) _ 12 2 ) L2 L L2 L (‘9..-).=(9.),E -—-—,——+— Q: -——+— Q (547) " ‘ th 79 21’] 70 21’ —~ _ 12 ) Where I is the second moment of inertia, 7V, is the penalty parameter associated to the DOF w, and 70 is the penalty parameter associated to the rotational DOFs 6V and 6V. 5.2.2 Bending Moment In this section we consider a concentrated bending moment applied at the tip, as plotted in Figure 51. As in the previous analysis, the symbolic results computed from the global system of equations are quite complex expressions. But, applying the same kind of simplifications, 142 groups of terms whose contribution to the solution is small have been eliminated. Figure 51 . Single plate element under bending moment applied at the free side Again, the solution obtained is proportional to the Timoshenko beam element solution for the analogous load case. 2 2 w2=w3§ —£—3———+—L— M=[-£—+-£]M (548) ”[1] 79 2Y1 79 L. 12 a i .. (6.). =(9..)'. E L +—‘- M=[i+—1—]M (549) ' ‘ ' ' Y1 70 5. 3 Penalty Parameter and Element Properties As for the plane stress quadrilateral, the obtained solutions are explored in order to find a relation between material and geometrical properties of the element and the penalty parameter. 5.3.1 Transversal Load The proper term of comparison for the achieved solution is the outcome provided by the same element without using the penalty interface frame. Appling the same type of 143 simplification used in on the penalty solution, the FE results are: ,.,, LSQ LQ L3 L w z + = -———+—— Q (550) t3h Y [(317 kGAH __ k— th 3Y[12] 2( ) (g). )“5 z- __I;_Q___ = [%]Q (551) 3 2y if 12 ’ is evaluated in order to underline the dependences between penalty The ratio w over w’ the two solutions. L3 L 1 L2 '1- —~—--——-+—— +~-+‘——-~ Q W pun) I} — 3Y1 kGA yw 76 r 1 1 L2 .. — = + + 552 w” L3 L L3 + L L3 + L ( ) —_—— + ___.___ —————— M __ ___.__ 3Y1 kGA Q 3Y1 kGA 7‘” 3Y1 kGA y, Now, if the penalty parameters 7W and 76 are substituted by: 2 “’3 141,; 3Y1 kG/i 2L2 ”‘3 _L:+___L__ 3Y1 kGA The ratio between the solutions becomes independent of material and geometrical (553) (554) properties of the element. penalty 1 1 1 ,..,.. =1+—+—=1+——- (555) w ' 2,3 2fl ,6 ["15 over (6),) is: malty The ratio (6y )p 144 2Y1 Setting: 2Y1 . = fl ' T We reach our goal. ( )pt'nully ‘ .- = 1+ -— 5.3.2 Bending Moment (556) (557) (558) The finite element solution computed without using the penalty interface frame is proportional to: W“: 2 L2 L_2 Y 13h 2Y1 12 pend/t 1' The ratio w over w'”” is computed to be: L2 L penal/y [E + —] M 2Y1 W . . z r y” = 1 + w/'[: £2— M L70 2Y1 Substituting penalty parameter 79 with: 145 (559) (560) (561) The previous ratio becomes dependent only on ,8. penalty W H, = 1 + -1— w The ratio (6),)Pmum over (9),)”; is: l L 1 ' malt - —*~ W" M (6.1),) I = _2¥I+Z’95 =1+ Y] (9,. )H: L M £79 ' _2Y1_ Setting: 79 = ,6 ° 211,1 we reach our goal. 6 perm/l)“ L17]— : 1 + _ (6...) 13 5.4 Building an interface Element for Plates (562) (563) (564) (565) (566) The Abaqus User Element Subroutine UEL for the plates doesn’t differ much from that for the 2-D quadrilateral plane stress elements. The subroutine for plate needs from the input file the same information as in the previous case. The parameters provided are: the nodes on the interface, thickness and Young’s Modulus of the two domains. Changes worth being mentioned are: (a) six DOFs per node are present; (b) the in- plane, transversal and rotational DOFs require different penalty parameters for optimal 146 behavior. The sixth DOF, t9: , corresponds to the rotation about the z-axis, often defined as the “drilling” degree of freedom. The statement in (b) results from our investigations on the penalty parameter. Six different penalty parameter are used, three for each domain. In particular: one, y“, is assigned to the in-plane DOFs (u, v), the second, y", , to the transverse DOF w and the third, 70 , is shared by all the rotational DOFs (6x, 6‘, 6: ). 5.4.1 Automatic Choice of an Optimal Penalty Parameter It is important to underline once more that we need only the right order of magnitude of the penalty parameter. Taking into consideration the earlier comment, we now recall the main results found in section 5.3. The penalty parameter related to the in-plane displacements, obviously, takes the same values assumed for the UEL dedicated to the 2-D quadrilateral plane stress element. We saw that under conditions of transversal load applied at the tip, the following value of the penalty parameter would make the accuracy of the solution to depend only on the value assigned to the parameter ,6. 2 w: - 567 7 fl _L: +__L__ ( ) 3Y1 kGA 7 —/3~ ”2 (568) " 21.1; 3Y1 kGA And: 2Y1 70 = :6 T (569) 147 Clearly, these requirements are in conflict and there are no simple relations between the two expressions for 79 to tell us which to choose. It has been discussed previously that the length of the element L is an unwanted presence in our expressions for 7. Thus, again we choose to approximate the dimension in the direction perpendicular to the interface with the one parallel to it. This means that in all the expressions it is assumed: L = h . 2 7‘“ _ fl 4h3 1 Yt th 2112 (570) ”zfl'4fi 1 __3_ + ____ Yt th And: t3Y 79 = flO—6- (571) When a bending moment is applied to the free side, in order to get the 79 independent from the element properties, we need: 7 = 13 ° -— (572) and: y. = fl g (573) We record another conflict in choosing 70- But in this case the first expression is clearly the one determining the bigger value of the parameter, so it is our choice. As we have just seen, afier the necessary substitution L = h , this expression becomes: 148 79 = fl -4— (574) At this point it is obvious that the two different conditions, transversal load and bending moment, could require different values of 70. A reasonable solution consists in choosing the parameter during the analysis using an IF clause. With all the necessary information about geometry and materials properties, the different values of ya can be compared. The bigger value will be selected in order to obtain: 6 penalty ( ) <1+-1— (575) (6)1115 — fl Table 11 summarizes all the relations among the penalty parameters and model characteristics. 149 Table 11. Expression for computing the penalty parameters in a plate element. All kind of loads Et ) 7a] fl.[-—bl—L 1 / Et ) yvl fl.(-ZLL I / 2 7 ’6 4b,3 + 1 WI Elt13 leltl 2b,2 fl'E 4b: +;—_] ya] bltl 1:101“ or ,6-5'65 £22 7112 fl [ b2 ) E4.) 7.2 fl[b—) 2 2 fl. wa 4b; + 1 E263 k202t2 2b,2 W ' 1 702 E26 1:202’2 5. 5 Numerical Results In this section the developed Abaqus User Element UEL for plate elements, which we will again call interface element, is employed to solve example problems. 150 5.5.1 Plate Bending Like the first analysis regarding the quadrilateral plane stress element, the following problem can be considered a patch test. A plate is divided in two domains discretized differently: 10 by 8 against 10 by 6. They are joined along the common interface through two Interface elements. Uniform bending moments, equal to 1 in magnitude, are applied at the opposite two edges. The four noded Abaqus shell element S4R is employed to mesh both domains. In-plane rigid body movements are avoided by applying the necessary constraints, moreover transverse displacements at the four comers are set to zero. Figure 52 plots deformed and original configurations of the plate. Others properties are: E1 = E2 = 1E6 , length = 4 , width = 1 and thickness = 0.1 . Figure 52. Original and deformed configurations of the plate This analysis would require the user to choose different orders of magnitude for the values of the penalty parameters related to transverse deflection w and rotational angles. However, the implemented ability of automatically choosing the appropriate penalty parameters for each DOF, saves this effort without sacrificing accuracy. Table 11 reports the results from the penalty interface FEM, obtained with ,6 = 151 1000 and ,6 = 10000, compared to the classical solution. Values shown are nondimensionalized by the exact solution. Table 12. Nondimensionalized results for the center deflection of the plate Masha Meshb w-fl=lE3 w'fl=1E4 Interface FEM 10x8 10x6 1.0004 1.0000 Outcomes confirm the reliability and the precision of the developed method. 5.5.2 Simply Supported Plate under Sinusoidal Load A simply supported plate is subjected to a load distributed over the surface according to the following expression: q = q0 singsinflfby (576) In which q0 represents the intensity of the load at the center of the plate. a and b are the plate dimensions, respectively, in the in-plane directions x and y. We assume the plate to be square with side dimensions equal to two, a = b = 2 . The value of q0 is taken to be 1. Taking advantage of geometrical and loading symmetries, only a quarter of the plate is analyzed. As shown in Figure 53, a fine mesh is used for a square near the center, while the remaining surface adopts a coarse mesh. The four noded Abaqus shell element S4R is employed to mesh both domains. Four Interface elements are used to connect the two domains. Others properties are: El = E2 = IE 6, V = 0.25 and thickness = 0.1 . 152 4y x % Figure 53. Penalty interface FE model of a flat plate. Figure 54 plots deformed and undeforrned configurations of the FEM seen from the bottom. Figure 54. Original and deformed mesh as seen from the bottom of the plate. 153 Figure 55. Distribution of traversal displacements w - Penalty interface from Figure 56. Distribution of traversal displacements w - Conventional FEM. 154 An analytical elasticity solution to this problem exists in the literature. It obviously will be our main term of comparison for validating our results. However, before performing this evaluation, we want to report the 2D maps for transverse displacement w produced by our interface element FEM and by a conventional 4 by 4 FE. Results are shown in Figures 55 and 56. These figures affirm that no discontinuities can be noticed across the interface. Moreover, the presence of the penalty interface frame doesn’t alter the proper distribution of displacements. Finally, the classical solution is compared to the interface element model. In Figure 57, the transverse displacements are reported along a straight line going from the middle of one side of the plate to the center. Since we are studying a quarter of the plate, this means that we are plotting w along one of the two sides of symmetry of the model. Even in the presence of a high gradient of deformation, the model behaves well. This further confirms our confidence that the developed method is both robust and accurate. 155 0.E+00 o 0.4 0.6 0.8 x 1.0 4.5-04 — —— Elasticity Solution 0 Interface Element -2-E-04 — Solution -3.E-04 . 4.13-04 . w -5.E-04 l Figure 57. Variation of the transverse displacement w along one of the two sides of symmetry of the model 156 CHAPTER 6 3D LINEAR BRICK ELEMENT In this chapter an interface element for three-dimensional solid 8-node linear brick elements is developed. Each of the eight nodes has three displacement components u, v and w; a generic configuration of the brick element is shown in Figure 58. Figure 58. T hree-dimensional solid 8-node linear brick element. 6. I A utomatic Choice of the Penalty Parameter The finite element model associated with this element can be found in several books [41-43] and is not difficult to derive. But, the task of computing and solving symbolically the system of equations associated with simple problems involving brick elements, overwhelms of the capabilities of the mathematical tools available to us. However, it seems reasonable to verify if the results obtained for the two-dimensional 4- 157 node quadrilateral element can be successfully applied to brick elements. It has been derived in section 4.6.1, a single expression for an automatic evaluation of the penalty parameter in the case of the quadrilateral element. The relation among the penalty parameter 7 and the properties of the finite element model, was found to be: Y 7 = fl-[flé—j (577) Where the Y is the Young’s modulus, b is the length of the element at the interface and while h is the thickness of the element. The same relation with a few obvious updates may be adopted for the automatic choice of the penalty parameter in the interface formulation for brick elements. Let’s consider the face of the brick element along the interface, see Figure 59. Face of the brick element along the interface Figure 59. Definition of a local coordinate system It is possible to construct a local coordinate system such that the x and y direction are in the plane identified by the face and the x direction correspond to the bottom side of the quadrilateral face. Then, h would be the side length along x and b the component along y of the other side departing from the origin. 158 Several numerical tests have shown that the relation (577) can correctly evaluate the penalty parameters. 6.2 Cubic Spline Interpolation Functions A 2-D interface element able to connect 3-D finite element meshes requires using two-dimensional interpolation functions [27-3 8]. Schumaker [30] and Prenter [31] proved that a two-dimensional version of the cubic splines, defined bicubic splines, can be constructed over a rectangular grid by taking the tensor product of one-dimensional splines. Their work also shows that the bicubic spline interpolation over a rectangular grid is smooth, having continuous second derivatives. In order to allow the connection of finite element meshes whose interface boundaries are not rectangular, the bicubic spline interpolation is mapped from a square region to the interface quadrilateral area. The only limitation is that the boundaries of the interface area must form a quadrilateral region. It is important to remember that by using several quadrilateral elements it is possible to mesh every kind of interface surface. Thus, by employing several interface elements it is possible to connect meshes of every type. For example, among the numerical results, we have been able to test the connection of two circular shafts of different diameters. The bicubic splines are constructed by taking the tensor product of one- dimensional splines in their general form described in section 4.2.], not the natural cubic spline over three points. Thus, the Abaqus user element subroutine associated with the 2- D interface element, allows the user to specify the desired number of nodes for the interpolation of the two-dimensional displacement field V. 159 6.3 Building an Interface Element for 8—Node Brick The Abaqus User Element Subroutine UEL for 8-node brick elements has some important differences with respect to the ones for the 2-D elements, but most of them are not visible to the user. As before, the parameters provided are: the nodes at the interface and Young’s Modulus of the materials composing the two domains. Nevertheless, a great amount of programming work has been required to deal with the much higher level of complexity brought by the additional dimension. In particular, a very difficult task was the integration of the matrices Gig, which requires dividing the integral over the interface surface in several integrals over smaller areas. These areas are the intersections of the bicubic spline mesh with the finite element meshes. As in the other cases, the UEL subroutine is written in Fortran 77. All the variables are defined to be in double precision format, in order to minimize numerical round-off errors as much as possible. 6. 4 Numerical Results In this section the developed Abaqus User Element, UEL, for 8-node brick elements is tested by several challenging problems. 6.4.1 3D Cantilever Beam with a Square Section A three-dimensional cantilever beam is meshed as composed by two domains. The domains are meshed independently, 5x5x10 and 4x4x10 elements, and joined by one interface element. The geometrical configuration of the problem is represented in Figure 60. Others properties are P =1000N and E, 2 E2 =1.MPa . 160 Figure 60. 3D Cantilever beam with a square section The only interface element used, is made by a uniform square grid of 36 points (6x6). For discretizing the two domains, solid 8-node linear brick elements are used. Two different load conditions are investigated: axial and bending loads. In order to compare the results obtained from the multidomain analysis, a conventional finite element model without interface was created for the entire beam and its mesh discretization, 5x5x20 brick elements, is that of the finer subdomain. These problems may be considered a form of patch test for the element developed herein. 6.4.1.1 Axial Load Under the condition of a uniform axial load applied at the free end of the beam, the interface element results are in complete agreement with the exact solution to the number of significant digits available, see Figures 61-65. There is no need to compare the 161 results with those obtained from a conventional mesh, since the solution to this linear problem is well known. Figure 61. Cantilever beam with a square section — Axial displacement U3 — 3D view +1.0003—35 +7.69ZE—Ud +1.53BE—D3 +2.308E-03 .3.0772-ua +3.8463—03 (A +4.6lSE-03 “ ‘5.3asz—03 +6.154E-03 .e.9zaarua +7.692E-U3 , .e.nezz—ua ' «9.2313—03 ‘1.uuus—uz Figure 62. Cantilever beam with a square section — Axial displacement U3 — Side view 162 AHA)“ 1-,/| (it Figure 63. Cantilever beam with a square section — Stress S33 Figure 64. Cantilever beam with a square section — Lateral displacement U] 163 Figure 65. Cantilever beam with a square section — Lateral displacement U2 6.4.1.2 Bending Load For the bending load case, the comparison is made to both the classical solution and a reference finite element solution. The mesh of 8-node linear brick elements cannot exactly recover the classical solution. In Table 13 the tip deflection for beam flexure is reported. The FE result, marked as Abaqus, is obtained from a traditional analysis using a compatible finite element model of the beam. Table 13. Tip deflection numerical results for beam flexure. Mesh Solution mesh Classical 164 VALUE .1745-36 .0953—01 .1B9E—01 .284E-01 .2388+00 .547E+DU .BS7E+00 .166E+00 .476E+00 .7BSE+00 .095E+UD .404E+00 .7138+00 Figure 67. Interface element model — Displacement U2 - (Classical deflection = 4. 0) 165 Figure 68. ABAQUS reference model — Displacement U3 Figure 69. Interface element model — Displacement U3 166 VALUE .BBBE+04 .9eea+ua .ueos+ua .1733+04 .267E+04 .360E+04 .5332+03 .saaa+03 .3sua+04 .267E+04 .1735+n4 .uaua+04 .9eas+o4 .8932+DL Figure 70. ABA QUS reference model — Stress S33 — (Classical = 6. 0E+4) VALUE —5.893E+04 . -a.9aaa+na “'_-—a.oeoa+ua -3.1732+04 -2.267E+04 -1.360E+04 —4.533£+03 ' +4.5332+03 .1.3eng+na +2.267E+Dd +3.1732+o4 .uaus+04 Figure 71. Interface element model - Stress S33 — (classical = 6. 0E +4) 167 Figures 66 to 69 demonstrate that the main components of displacements, U2 and U3, do not show any change in their values across the interface. The displacement field Ul is not reported since it is in the order of 1.E-l4, almost zero as expected. Only light discontinuities in the values of main stresses 033 at the interface are present, as plotted in Figures 70 and 71. 6.4.2 3D Cantilever Beam with a Quadrilateral Section The interface element for solid 8-node linear brick elements is able to connect finite element meshes whose boundaries at the common interface are not rectangular. To test this capability, a three-dimensional cantilever beam with a distorted quadrilateral section is meshed as composed by two domains. The domains are meshed independently, 5x5x10 and 4x4x10 elements, and joined by one interface element. The geometrical configuration of the problem is represented in Figures 72 and 73. Figure 72. 3D Cantilever beam with a quadrilateral section 168 The Young’s modulus of the materials composing the two domain of the beam are E, = E2 =1.MPa. A transversal load is applied in the 2-direction to the bottom four nodes at the tip of the beam and its total value is P =1000N . (-0.75,1.5) (090) Figure 73. Quadrilateral section of the beam. In Figure 73 are shown at the same time the two different discretizations applied to the section of the beam; 5x5 and 4x4 brick elements. Only one interface element is adopted; it is made by a uniform quadrilateral grid of 36 (6x6) points. For discretizing the two domains, solid 8-node linear brick elements are used. One load condition is reported: a bending load. Axial load results are not reported for brevity, but they are as good as in the previous case. In order to compare the results obtained from the multidomain analysis, a conventional finite element model without interface was created for the entire beam and its mesh discretization, 5x5x20 169 brick elements, is that of the finer subdomain. Figure 74 shows the deformed configuration of the mesh with the interface element. Figures 75-80 demonstrate that all components of displacements are continuous across the interface. As before, it is possible to notice limited discontinuities in the values of main stresses 0'33 at the interface, as plotted for the two lateral views in Figures 81 and 84. Figure 74. Deformed mesh 170 VALUE 792E-36 UTlE-Dl 141E—01 212E-01 2835—01 353E—01 424E—Ul 495E-Ul SESE—Ol 636E—01 071Eo00 178E000 ZBSE‘UU 392E+00 Figure 75. ABAQUS reference model — Displacement U2 VALUE 7915-36 DTUE-Ul 1415-01 leE—Dl 282E—01 352E~01 423E-01 493E—01 563E—Dl 634E—01 070E900 177E+DU ZBSE+DO 392E+00 Figure 76. Interface element model — Displacement U2 l7l VALUE 513E-35 387E—02 774E-02 160E—02 S47E—02 193E-01 4322—01 671E~01 909E-01 14BE—01 387E-01 625E—01 8645—01 lU3E-DI Figure 77. ABA QUS reference model - Displacement UI VALUE 513E—35 386E—02 7735—02 159E-02 5455-02 193E—01 432E-01 6705—01 909E-01 148E-01 .JBBE-Dl 625E-01 864Ev01 102E—01 Figure 78. Interface element model — Displacement U] 172 VALUE 11113-01 983E-02 fiSGE-DZ 729E-02 602E-02 748E-03 GEEK-02 779E-02 9062-02 033E-02 OISE-01 229E-01 441E-01 654E-Ul Figure 79. ABAQUS reference model — Displacement U3 Figure 80. Interface element model — Displacement U3 173 VALUE .208E904 .790E004 .372Es04 .53SE¢U3 355E403 174E#03 .0071403 .1BBEOD3 137Es04 .555I404 .973E404 .391E004 8092404 227Eo04 Figure 81. ABAQUS reference model — Lateral view I - Stress S33. VALUE 209E¢U4 .7903404 ,372!~04 ,536!~03 . 355E413 174E003 007E903 .188E+03 .137E904 55515904 .973h04 391E404 .809E004 .227E004 Figure 82. Interface element model — Lateral view I - Stress S33 174 VALUE -2 203E~04 —1 7902.04 ,., -1.372:.o4 ‘ ’. -9.sas:.ua ' -s.ass:.na , -1.174:.03 , g, .3 007:.03 'Z-o7 1332.03 .1 1312.04 91.555EOD‘ .1.97az.o4 , .2.391z.o4 .2.aogz.04 »3.227:.04 Figure 83. ABAQUS reference model - Lateral view 2 — Stress S33 VALUE .208E604 790Eo04 .37ZEoD4 S36Eo03 355E903 174E403 007E003 .188E003 131E904 .555E604 .973Eo04 .391E004 ,809EoD4 .227EOD4 Figure 84. Interface element model — Lateral view 2 — Stress S33 175 6.4.3 3D Shafts Assembly Connecting shafts of different diameters is a problem of interest to the industry. It would be very useful to the designers to be able to connect shafts of varying dimensions without having to remesh the each section of the model each time. The presented interface technology would be very useful in solving this kind of problem, allowing a faster search for the optimal design. Figure 85. Shaft assembly. For the test analysis the configuration shown in Figure 85 has been considered. 176 The larger shaft has a diameter of 2 m and an axial length of 5 m, while the slender one has a diameter of l m and an axial length of 10 m. The larger shaft has the one face encastred and is meshed with 480 solid 8-node linear brick elements. The slender shaft has the load applied to the free end and is meshed with 800 solid 8-node linear brick elements. In Figure 85 is possible to notice the difference between the meshes at the interface of the two domains. The circle of 1m in diameter identifying the common interface, is discretized with 64 elements for the first shafi versus 80 for the second slender shaft. For comparison purposes a reference conventional model without interface elements was constructed. The level of refinement of the mesh is comparable to that of the slender shaft in the multidomain analysis. Inside the common interface surface 17 quadrilateral areas are identified and an equal number of interface elements are used to enforce the continuity constraint between the incompatible meshes. Both shafts are made of isotropic material with Young’s modulus equal to E,=E,=1.MPa. 6.4.3.1 Bending Load A transversal bending load is applied to the bottom node at the tip of the beam, its direction is 2 and its value is P=100N . The comparison is made to the reference conventional model without interface elements. Figure 86 shows the original and deformed configurations of the mesh with the interface element. 177 ORIGINAL MESH Figure 86. Original and deformed mesh under bending load. Figures 87 to 92 demonstrate that all the components of displacement are continuous across the interface. Also in the stress maps for the reference conventional model and for the finite element assembly with the interface elements, it is hard to notice any difference at all; see Figures 93 and 98. 178 Figure 88. VALUE .145E-04 .6922—05 .929E-05 .167E—05 .4058-05 .643E-DS .8112-06 .BllE—U6 .643E-US .405E-05 .167E—05 .929E—05 .692E—US .14SE—04 ABAQUS reference model - Displacement U1 VALUE .145E-04 .6BBE-US .9263-05 .1653-05 .404E—DS .642E-DS .BUVE—DE .BD7E-06 .6428-05 .404E-05 .IGSE-DS .927E-05 .SBBE-OS .145E-U4 Interface element model - Displacement U1 179 VALUE -5.727E-02 —5.286E—02 —4.B4SE-DZ —4.QUSE—02 —3.965E—02 —3.524E—D2 ~3.084E—02 —2.643E-02 —2.203E—02 —1.762£—02 —1.322E—02 ‘B.811E—03 -4.dUSE-U3 .9458-38 Figure 89. ABA QUS reference model - Displacement U2 VALUE .736E—02 .295E—02 .BS3E—02 .412E—02 .971E-02 .53OE—02 .DBQE-DZ .6478-02 .206E-02 .7658—02 .324E—02 .8248—03 .412E-U3 .759E-3B Figure 90. Interface element model — Displacement U2 180 VALUE .847E-U3 .654E—03 .460E—D3 .266E—03 .073E—U3 .879E-D3 .BS7E-Dd .D79E-04 .7028-03 .BBSE—DB .DEBE-UB .2628-03 .4768-03 .67DE-U3 Figure 91. ABAQUS reference model - Displacement U3 Figure 92. Interface element model — Displacement U3 181 Figure 93. + H ++++ HHH N + N Figure 94. VALUE .929E+Ul .224E+01 .5198+01 .148E+00 .BQUE+01 .595E+01 .299E+01 .UDUE+02 .271E+02 .541E+02 .8123+02 .082E+02 .3538+02 .623E+02 ABAQUS reference model — Stress S11 VALUE .913E+01 .210E+01 .SD7E+01 .0393+UD .8993+01 .6023+01 .3058+01 .001E+UZ .271E+02 .5413+02 .BlZE+02 .082E+02 .3SZE+02 .6223+02 Interface element model -— Stress S11 182 Figure 95. ABAQUS reference model - Stress S22 VALUE -2.609E+02 —1.528E+02 —4.47sr+01 .+e.331a+01 +1.7142+02 7.2.7952+o2 , +3.87EE+02 V-+4.9sag+02 +6.037E+02 +7.117s+02 +8.1933+02 .2793+02 .036E+03 .1443+03 + + H w Figure 96. Interface element model — Stress S22 183 Figu. VALUE .498E+03 .268E+U3 .037E+U3 .068E+02 .763E+02 .458E+02 .153E+02 .151E+02 .456E+02 .761E+02 .066E+02 .037E+U3 .268E+03 .4983+03 Figure 97. ABAQUS reference model — Stress S33 VALUE .4293+03 .209E+03 .893E+02 .69SB+02 .4978+02 .298E+02 .100E+02 .09BE+02 .297E+02 .4953+02 .693E+02 .8913+02 .2093+03 .429E+03 Figure 98. Interface element model — Stress S33 184 6. 4.3.2 Torsion A torsion load is applied to the free end of the slender beam. Both shafts are made of isotropic material with Young’s modulus E =1.MPa and Poisson ratio v = 0.3. The comparison is made to the reference conventional model without interface elements. Figure 99 shows the deformed configuration of the mesh with the interface element. ’f/ j‘I‘ ”/ Figure 99. Deformed mesh under torsion load. Figures 100 to 113 shows that all the results from the multi-domain with interface elements analysis are in excellent agreement with those from the reference model. 185 Figure 100. Figure 101. .6488-02 .947E—02 ABAQUS reference model — Displacement U1 VALUE .948E~02 .648E-02 .349E—02 .049E—UZ .493E-03 .496E—U3 .499E-03 .4993—03 .496E-03 .493E-03 .049E-02 .349E-02 .6482-02 .SdflE-UZ Interface element model - Displacement U1 186 Figure 102. l l J:- «I + A Figure 103. VALU E .947E-02 .64BE-02 .348E—02 .049E-02 .4903-03 .494E—03 .49BE-03 .498E-03 .494E-U3 .490E-U3 .049E-02 .3488—02 .6483—02 .947E~02 ABAQUS reference model — Displacement U2 VALU E .948E-02 .648E-02 .349E-02 .0498—02 .493E-03 .496E-03 .499E-03 .499E-03 .496E-03 .493E—03 +1. .3498-02 .6488-02 .948E-02 0493-02 Interface element model — Displacement U2 187 VALUE .BdTE-US .486E—05 .1258—05 .764E—DS .403E-05 .UdZE-US .BUSE-OG .806E-06 .U42E—US 4D3E—05 .764E-05 .125E-US .QBGE-DS .847E—US Figure 104. ABAQUS reference model — Displacement U3 VALUE .847E-US .486E-DS .12SE—OS .7643—05 .403E—05 .UdZE-US .BDGE-Ub .BD6E-06 .UAZE—US .4038-05 .7643-05 .1253-05 .486E-05 .847E-05 Figure 105. Interface element model — Displacement U3 188 VALUE .115E+03 .7QDE+03 .464E+03 .1393+03 .13QE+02 .BBIE+02 .627E+02 .627E+02 .881E+02 .134E+02 +1.139E+03 .464E+03 .790E+03 .1158+03 ||lll mHHHN A H 7+ H 4. a 4. tn 1. H 4. H Figure 106. ABAQUS reference model — Stress S12 VALUE .1153+03 .7908+U3 .464E+03 .139E+03 .134E+02 .881E+02 .627E+02 .627E+02 .881E+02 .134E+02 .1398+03 .464E+03 .790E+U3 .115E+03 Figure 107. Interface element model —- Stress S12 189 VALUE .6053+03 .358E+03 .lilE+03 .64ZE+02 .173E+02 .704E+02 .2358+02 .235E+02 .7UQE+02 .173E+02 .642E+02 .111E+03 .3583+03 .605E+03 Figure 108. ABAQUS reference model — Stress S11 VALUE .605E+03 .BSBE+03 .1112+03 .GQZE+02 .173E+02 .7043+02 .235E+02 .2352+02 .704E+02 .173E+02 .6423+02 .1113+03 .3583+03 .GDSE+03 Figure 109. Interface element model — Stress S11 190 VALUE .605E+03 .3SBE+03 .111E+03 .6423+02 .17BE+02 .7048+02 .235E+02 .ZBSE+DZ .704E+02 .173E+02 .642E+02 .111B+03 .3SBE+U3 .605E+03 Figure 1 10. ABAQUS reference model — Stress S22 VALUE .605E+03 .358E+03 .111E+U3 .642E+UZ .173E+02 .704E+02 .235E+02 .2353+02 .704E+02 .173E+UZ .642E+02 .111E+U3 .3SBE+03 .605E+03 Figure 111. Interface element model — Stress S22 191 Figure 112. VALUE -2.026E+02 -1.714E+02 -"~1.403E+02 —1.091E+02 —7.7923+01 —4.67SE+01 —1.558E+01 +1.SSBE+01 +4.675E+01 +7.792E+01 +1.0913+02 +1.403E+02 " .+1.7143+02 +2.026E+02 ABAQUS reference model — Stress S33 Figure l 13. Interface element model — Stress S33 192 CHAPTER 7 INTRODUCTION TO THE MODELING OF DELAMINATION GROWTH Composite materials are widely used today because they enable high performance and low weight structures. However, due to the low tensile strength of the matrix, composite materials may experience delamination damage under service loading. This damage deteriorates load carrying capacity and can lead to premature fracture of the composite component. Extensive research has been performed in the attempt to better understand the major failure mechanisms for laminated composite delamination initiation and growth. The present chapter gives a brief review of the state of the art on the subject. The aspects considered first are the main causes of delamination. Then, computational techniques used in the literature to predict delamination initiation and growth are reported. 7. 1 Causes of Delamination In this section the main factors responsible for arising of interlaminar stresses able to originate delamination are examined. 7.1.1 F rec-edge stresses Delamination along the free edge of composite laminates has been studied since 193 the 1970s. Since then an important amount of work has been reported on the free edge problem in laminated composites, indicating that free edge delamination is attributed to the existence of interlaminar stresses near the free edges. These stresses occur due to the mismatch in engineering properties, i.e. mismatch in Poisson's ratio and coefficient of mutual influence between layers. In summary, three classes of interlaminar stress problems exist. Laminates of the type of [i0] exhibit only shear-extension coupling (no Poisson mismatch between layers), so 7;, is the only nonzero interlaminar stress. [0/90] laminates reveal only a Poisson mismatch between layers (no shear-extension coupling), so I” and 0': are the only nonzero interlaminar stresses. The mismatch in V“, give rises to interlaminar normal 0': and shear r“. stresses. Laminates that are combinations of i0, 0° and 90°, for example [tel/:02] laminates, show both shear-extension coupling and r and 0'z interlaminar stresses. Poisson mismatch between layers, so they have 2'3, ,2 Since high interlaminar shear is more likely to exist in the cases where the angle plies [i0] are adjacent to each other, one should try to avoid placing them together in a laminate. The other parameter affecting the interlaminar stresses is ply thickness. Thick plies tend to encourage higher interlaminar stresses thus causing early delamination. 7.1.2 Matrix cracks Delamination is expected to occur at the interface where the predominant interlaminar stress component is computed unless there is a perturbation in the state of stress, probably induced by damage, prior to the onset of delamination. In most cases the delamination propagates in the axial direction along the same interface. In some other 194 cases, the delamination in the course of propagation turns into the neighboring ply. Several studies showed that under applied compression, all delamination initiated prior to any transverse cracking, in fact, in most of the cases the specimen failed without revealing any transverse crack. Under tensile loading, however, the delamination in most laminates, especially those containing 90° plies, is preceded by a number of transverse cracks. Because of the presence of the transverse cracks, the location of delamination is not unique. The path of delamination along the axial direction varies widely and depends upon the size and location of transverse cracks, type of laminate, material system, etc. The formation of the delamination in [3:0/90]s and [0/i0/90]S laminates with 0 in the 10°- 45° range generally occurs in the following manner. Upon a further increase in applied load after transverse cracking, several delaminations in the form of axial cracks form at the mid-plane or interfaces between the +0/-0 layers or the -0/90 layers. Delaminations at the last two interfaces are nucleated by the transverse cracks. As the load increases, the delamination opens up further and either continues to run along the same interface or changes to other interfaces, normally at transverse cracks. At this time the delamination also propagates rapidly toward the middle of the specimen. A finite element study [44] of transverse cracking investigated the behavior of two types of cross-ply laminates, [90/0]s and [0/90]s subjected to an extension 80. It was found that when the load increases, numerous transverse cracks running across the thickness of the laminate perpendicularly and possessing an approximately uniform spacing along the length of the laminate will occur. These transverse cracks terminate perpendicular to the ply interface because of the stronger 0° ply and then tend to develop delamination cracks along the ply interface if the load increases further. The ply material 195 is graphite epoxy T300/5208. For the [90/0]s the effects of the thickness of 90° ply on the stress distribution along the x-axis has been studied. When the thickness of the 90° ply increases, the strength of singular stresses near the crack tip is enhanced. A similar dependence on the thickness of the 90° ply is found for the [90/0]s laminate. Therefore, the thickness of 90° ply has an important influence on the growth of the delamination crack. 7.1.3 Impact When a laminate is hit by an object or projectile, the material directly under is compressed and translates laterally in a time frame much less than that which is required for the overall response of the structure. The highly localized deformation gradient causes large transverse shear and normal stresses that can cause the damages to propagate and even failure of the laminate. Another effect of the impact is the creation of a compression stress wave, which travels from the impact surface through the thickness of the laminate. This wave is reflected from the back surface as a tension wave, which can cause failure at the first weak interface, resulting in chipping or splintering of the parts of the rear ply. Both internal stress waves and local out-of-plane deformations may initiate delamination at interfaces where there is major change in the angle between plies. The amount and type of the damage in the laminate depends upon the size, type and geometry of the laminate, impact energy and the loading on the laminate at the time of impact. At relatively low impact velocities, the laminate can respond by bending and failing either by shear resulting in delamination or flexural failure depending on short or long beam, respectively. At higher velocities somewhat different damage modes occur which may be caused by the combination of stress waves and the out of plane deformations. The 196 damage may result in delamination, fiber failure, matrix cracking, etc. If the velocity of the projectile is fairly high, the laminate acts as relatively rigid resulting in shear out and complete penetration of projectile. 7.1.4 Residual thermal stresses and moisture Residual thermal stresses, caused by temperature and moisture changes, are always present due to cool down of the laminate from the elevated curing temperature and these may have considerable influence on interlaminar stresses. 7. 2 Experimental Studies on Delamination Growth Extensive research has been reported on the delamination testing of various composite materials. The majority of the work has been concerned with the determination of the critical strain energy release rates G. The reason why researchers are focusing on G is that, according to the fracture mechanics approach, delamination growth depends upon the stress state of the crack tip which is governed by the stress intensity factors K], K” and K,” or the strain energy release rates G], G” and Gm. Thus, it is essential towards the understanding of the reported results, to define the strain energy release rates G and to describe the testing equipment used to determine it. The goal of this section is to address this issue. 7.2.1 Definition of the strain energy release rate G Energy rate analysis of the effects of flaws historically preceded crack-tip stress field analysis. The Griffith Theory [45] and later modifications by Irwin [46], termed the Griffith-Irwin Theory, made use of this approach [47-50]. Basically, these methods use 197 an energy balance analysis of crack extension. The total elastic energy made available per unit increase in crack surface area is denoted by G for the linear-elastic case. Physically, G may be viewed as the energy made available for the crack extension processes at the crack-tip as a result of the work from displacements of loading forces and/or reductions in strain energy in a body accompanying a unit increase in crack area. Alternatively, G can be regarded as a "generalized force" based on the potential energy change per unit forward displacement of a unit length of crack front, which results in G being defined as "the crack extension force". Following this line of argument, it is not difficult to show that for linear-elastic conditions: G = ___5U7‘(ANA) (578) 6A and G = 4.9211519 (579) 6A and G = 51-65 (580) 2 6A where Ur is the total strain energy in a cracked body with a crack area A. Ur is alternately expressed in terms of A and load point displacements, A,, or in terms of A and loads, P,-. In the last equation, C is the elastic compliance and the equation is written for a single loading force, but may be generalized for several forces. The G implied by the previous equations is the average value along a crack front weighted for the extent of crack extension involved for each increment of crack front in the three-dimensional sense. In two-dimensional situations. such as uniform extension of a straight-through crack in a thin plate subject to extension, G may be viewed as the value of a point quantity along the crack front. In several studies the critical value of Ge, required to cause the delamination 198 growth, is also defined as “delamination fracture toughness”. The Griffith criterion for fracture states that crack growth occurs if the total energy of the body remains constant or decreases as the crack length increases. In elastic solids, if W is the energy required for crack growth, then according to Griffith the necessary condition for crack growth can be expressed as: 02% (581) It is common to replace Egg/— by R the "crack resistance". It needs to be underlined that the stress intensity factors K and the strain energy release rates G are not independent parameters, as the work of Irwin proved: the “energy- balance” approach to the characterization of fracture is equivalent to the “critical stress intensity factor” one. That means there exists a relationship between the strain energy release rate and the stress intensity factor: G oc K 2 (582) The constant of proportionality in the previous expression is a function of the elastic constants of the materials. This relationship relates the crack-tip stress field and the “energy-balance” criterion for crack growth, which can be interpreted in terms of a critical K value that is required for the crack to enlarge. Generally, researchers prefer to adopt the maximum total strain energy release rate criterion, which predicts the delamination will grow when G reaches the critical value G6. Another commonly recognized assumption is that the delamination can be characterized by linear elastic fracture mechanics (LEFM). The fundamental assmnption of LEFM is that the behavior of the cracks, whether they grow or not, and how fast they 199 grow, is determined only by the stress field at the crack tip. 7.2.2 Test specimens for evaluating G Interlaminar damage mechanism of composite materials results generally from mode 1, mode II or mixed mode I+II loading. That is why the most adopted test specimens have been developed for studying crack propagation under the mentioned modes. A classification of the testing configurations can, consequently, be derived. In all these specimens, the aim is to measure critical strain energy release rate Gc during onset or extension of delamination. Specimens are calibrated by means of a relation between G and compliance C. 2 = £29 (583) 2b 6a where a is the crack length, b the width of the specimen, P the applied load and C is the compliance defined by the slope of the load P versus displacement 5 curve for the specimen. 5 C = — 584 P ( ) 7.2. 2. 1 Test specimens for Mode 1 Pure mode 1, opening mode, caused by interlaminar normal stresses, is one of the most common delamination mode. Several investigators, i.e. [51-57], agree in recognizing the double cantilever beam (DCB) as an accurate pure mode I specimen. A DCB, shown in Figures 114, is tested under displacement controlled conditions and P versus 6 curves are obtained for various crack lengths a which are then used to obtain C and (SC/6a. The value of Ge may then be determined using the critical load PC for 200 the initiation of delamination and the values of JC/da. Also, assuming DBC specimen to consist of two cantilever beams joined at the end of crack tip, the compliance may be expressed as: _ 8a3 bEh3 (585) where E is longitudinal modulus and h is the thickness of each beam. Substituting this into the general expression for G, one obtains: _:LP%: G _ 2ba (586) The beam relation can be used as long as the shear deformations can be ignored and the deflections are small, otherwise some corrections should be applied. 2h if ‘ I Figure 114. Double cantilever beam (DBC) specimen Only one other type of mode I specimen has been found [58] in the present review: the free edge delamination (FED). It seems FED specimen proposed provides a viable approach to determining delamination initiation, as it is relatively simple test specimen. However, because of the complex stress state in the free edge zone, the FED 201 specimen does not necessarily produce a pure mode I delamination. Therefore DCB is considered as most suitable to determine G/c. 7.2.2.2 Test specimen for Mode 11 To characterize the mode 11 delamination, shearing mode, the end notched flexure (ENF) and end-loaded split (ELS) specimens have been proposed, i.e. [51-57]. ENF specimens are projected as developing pure mode II behavior. The specimen is the same as the DCB but it is loaded in a three-point flexure which results in almost pure in plane shear delamination growth mode. 2h ‘l ‘ 1! ‘,_ 1’. Figure 115. E nd-loaded split (ELS) specimen For the mode II end-loaded split (ELS) test, the specimen is held at one end and loaded at the other, as shown in Figure 115. For a given deflection and curvature, the arms of the beam are subject to the same bending moments. 7. 2. 2.3 Test specimen for Mixed Mode I +11 The present study reported that the mixed mode delamination, i.e. [51-57], has been explored by use of the following tests: fixed-ratio mixed mode FRMM, mixed mode bending MMB, cracked-lap shear test CLS and modified notched flexure MENF. 202 A“ 41 V Figure 116. F ixed-ratio mixed-mode(FRMM) specimen. The fixed-ratio mixed-mode FRMM test containing a symmetrical crack is shown in Figure 116. For this symmetrically cracked FRMM test, the ratio of mode I to mode II loading is approximately constant throughout the test at 4:3; hence, the value of Gl/G” is 1.33. The mixed mode bending MMB specimens can provide a wide range of mixed mode I/II ratios by adjustment of the loading and lever fulcrum position. Using the cracked-lap shear test CLS specimen, the P-5 curves may be obtained for various crack lengths and 6C/6a determined. The substitution of PC and o'C/é‘a in equation for G give the value of Gc. CLS provides the total energy consisting of modes I and II. 7. 3 Delamination Criteria Initiation and evolution of the delamination may be predicted by a fracture mechanics approach or by interlaminar strength, introducing an interface constitutive law 203 between the layers. The two approaches are investigated separately in the next two sections. Since interlaminar stresses r r,: and 0': are responsible for delamination, the accurate study of delamination requires a complete 3D stress analysis. However, it is difficult to calculate interlaminar stresses analytically. Also, it should be recognized that the stresses at the delamination front are singular. That’s the reason why finite element analysis plays a dominant role in developing new criteria for predicting delamination. Consequently, the great majority of the works reviewed take advantage of the finite element method in some way. 7.3.1 Fracture Mechanics Based As it has been discussed previously the fracture mechanics approach uses strain energy release rate as a parameter for assessing delamination, initiation and growth. In fracture mechanics approach, strain energy release rate per unit area delaminated is evaluated and compared with Gc. Thus, one of the main efforts of the researchers is to correctly evaluate strain energy release rate. There exists few well-tested numerical methods for obtaining G and they will be explained in the next subsection. 7.3. 1.1 Numerical evaluation of the strain energy release rate Different methods for evaluating the strain energy release rate G have been adopted in conjunction with finite element analysis: the virtual crack closure technique (VCCT), the crack closure method (CCM) and the J-integral method. The J-integral method can be used to analyze both linear and nonlinear problems. However, this method can’t easily be applied to mixed-mode fracture problems. This is an important limitation 204 since delamination growth is a three-dimensional phenomenon with a general mixed- mode type of failure. The virtual crack closure technique (VCCT) does not have this restriction and, being simple and accurate, is almost universally adopted even if it can be used only for linear elastic problems. This choice is supported by several experimental results that well match numerical models based on the VCCT (i.e. Prathan and Tay [59]). The VCCT has been successfully used with both plate and 3D brick finite element models. Calculation of the Strain Energy Release Rate by V C C T In this section, classical definitions used in the calculation of the strain energy release rate by the VCCT method are reviewed. In finite element analysis, the total strain energy release rate could be calculated from two finite element solutions, one with a crack length a and another with a slightly different crack length a + An. Using the difference form of the previous equation gives: G = ___U...:a- U0 (587) Where U and U, are the total strain energies of the finite element models with crack a + Au lengths a+Aa and a, respectively. The value of G calculated from this relation is . . Aa a351gned as the stram energy release rate for a crack of length a + —2—. However, because of the difference approximation, the value of Aa needs to be very small. The formula is not attractive because (a) the computation involves differences of large numbers divided by a small number and consequently the resulting value is not accurate, and (b) this procedure involves two finite element runs. Hence, alternate forms of the formula for determining G have been proposed in the literature. These are based on the virtual crack 205 closure technique (VCCT) proposed by Irwin [46]. Elastic strain energy released during an incremental crack extension is equal to the work done in closing the incremental crack. When the crack is opened, the work done to close the crack opening could be written as: AW -—1—LM6'Azida (588) 2 A1? is the relative displacements between the crack surfaces along Aa when the crack is closed. On the basis of the crack closure equivalence, the available energy release rate G for a crack size a is expressed by: G = lim —1— Ma‘ A22 da (589) Aa—>0 2M Substituting in this equation the components of the surface stresses 6' and the corresponding relative displacements Afr leads to the components of G. In 2D problems modeled in the x-z plane, the mode-I and mode-II strain energy release rates are obtained from: . 1 G, = 1311—31555 fa,(x,0)w(x—Aa,0)dx (590) G - lim—l—fao' (x 0)u(x—Aa 0)dx (591) ll Arr—’0 2A0 x: -’ ’ Thus, the total strain energy release rate is: G = G, + G,, (592) Total In the previous expressions, w and u are, respectively, the crack opening and sliding displacements; 0', and 0'...- are the normal and shear stresses ahead of the crack tip. In 3D analyses, the strain energy release rates can be calculated using VCCT as: 206 Thu Wh . .1-‘+15'_1‘ Au G, = 131—11) 2Aa5y J; [L 0'_.(x,y,0) w(x — Aa,y,0) dx] dy (593) G — lim 1 [‘5] [Ma (x 0)u(x—Aa 0)dx]d (594) ll Arr—+0 2Aa6y y x: 9 y’ 9)” y G — 1' 1 ]i[ A" ( 0) ( ——A 0 dx d 595 III _5312102Aa5y y L 0y: xaya V x any» ) y ( ) Thus, the total strain energy release rate is: Glirtal : G] + G1] + Gm (596) Where w, u and v are the crack opening, sliding and tearing displacements, respectively. 6 y is the width of the crack front over which crack closure occurs, and 0', , 01,: and G,,, are the normal and shear stresses on the crack plane and ahead of the crack front. The crack closure representation is convenient for adaptation in numerical computation. That is another reason why several methods were proposed in the literature to calculate the strain energy release rates in finite element analyses using the VCCT. In the finite element method, continuous stress and displacement fields are approximated by their nodal values. Figure 117 (a) represents the finite element model near a crack-tip region. A crack of length a is shown with the crack tip at node 0. An incremental crack extension Aa is introduced replacing the crack tip node 0 with two separate nodes m and n as depicted in Figure 117 (b). For the new crack geometry the finite element solution for the nodal displacements (um,vm,wm) and (an, v", w") are found for nodes m and n, respectively. 207 Z a a E0 ‘P ._—+x y/ r s (a) m p x Y‘/// n (b) ll A a+Aa l‘ ' Figure 117. 2D FEM before crack opening (a) and after crack opening (b) The corresponding load F required to take the two bodies in contact for the ()9 delamination length a, is approximated to be the load F p obtained from the same analysis at the node p. Therefore, it is assumed that the load distribution at the interface near the delamination tip remains the same during delamination extension from a to a + Aa. This assumption requires the finite element mesh to satisfy certain conditions: the mesh is symmetric on the crack plane and about the crack front, the element size Aa at the crack front should be small and normality of the FE mesh is maintained near the crack front. These requirements are often easy to satisfy. Thus, the work required to close the crack opening is approximated by: 208 1 23,, AWE [E,(um—un)+Fv(vm—vn)+Fz(wm—wn)] (597) Where Fx, F, and F, are the components of the nodal forces at p. BI, is the equivalent width apportioned to node p: the average of the dimensions in the y direction of the two adjacent elements sharing the node p. Consequently, the energy release rate for the crack extension is approximated by: 1 I: ZACBP [17: (Wm _wn ):| (598) I G,, = 2MB. [1”; (u. — u. )1 (599) , 1 (II/I : flag-[F1 (vm — vn )] (600) In the literature, the virtual crack closure technique (V CCT) is employed in two ways: a) to evaluate the total strain energy release rate G across the width of the crack for various values of the crack dimensions, and b) to predict the crack growth. In the second case the total strain energy release rate G is computed for all the nodes in the original delamination front. The nodes with the largest values of G are released, simulating the advancement of the delamination front. Then, the strain energy release rate is recalculated along the new delamination front. The maximum value of G is found and the corresponding nodes are released. This procedure of determining G along the delamination front and subsequently releasing the corresponding nodes is repeated in order to simulate the progressive failure process. Calculation of the Strain Energy Release Rate by CCM The crack closure method CCM doesn’t differ much from the virtual crack 209 c105 11111 011 incr sep ..| 11111 101 Ct]; 1h( 011 H1 CCI of closure technique VCCT. However, since the forces and displacements of the system must be carried out at two separate configurations it requires two finite element analyses. That’s why it is much less used than the virtual crack closure technique VCCT. Figure 1 17 (a) represents the finite element model near a crack-tip region. A crack of length a is shown with the crack tip at node 0. Like in the VCCT approach, an incremental crack extension Aa is introduced replacing the crack tip node 0 with two separate nodes m and n as depicted in Figure 117 (b). For the new crack geometry the finite element solution for the nodal displacements (um,vm,wm) and (un,vn,wn) are found for nodes m and n, respectively. The crack opening is then collapsed by applying equal and opposite forces at nodes m and n such that they common displacements match those found earlier in 0 before it opens. Thus, the work required to close the crack opening is approximated by the same equation as in the VCCT: Awg—L[F,(um—u,,)+F,(vm—v,)+F,(w,,—w,,)] (601) P However, here 17,, F, and F_. are the components of the nodal forces at o. B], is the equivalent width apportioned to node p: the average of the dimensions in the y direction of the two adjacent elements sharing the node p. Consequently, the energy release rate for the crack extension is approximated by: G, = 2MB, [17, (W, — w, )] (602) 1 G,, 2M8. [F.. (u. — u. )1 (603) 1 -M[F.(vm —v,.)] (604) 7. 3. 1.2 Fracture mechanics based delamination criteria As reported previously, the virtual crack closure technique (V CCT) is very often adopted for studying delamination initiation and growth. In 1995 ku, Kao and Chang [60] conducted a set of experimental tests and finite element simulation in order to assess a delamination fracture criteria for composite material. For their investigation they employed the following specimens: the double- cantilever beam (DCB) test for mode I; the end-notched flexural (ENF) test for mode 11; cracked-lap shear test CLS and the modified end-notched flexure (MENF) tests for mixed mode I+II. The finite element analysis implementing the VCCT was necessary in order to separate the total value of G obtained by the experiments into the sum of each mode. The authors assumed that unlike cracks in homogeneous bodies, the interface cracks always induce opening, shearing and tearing mode fracture simultaneously for a single mode loading. Hence, unlike the unidirectional composite specimen, it would not be guaranteed that the DCB test may measure G, since the failure of DCB specimen may not only be determined by G,, but may also be influenced by G”. They had thought that before getting any experimental data from DCB test, one should only treat it like mixed mode failure test. But, results from the DCB tests confirmed that G, is far bigger than G11 and Gm, and the test is an almost pure mode I delamination. Same result was found about the ENF test for mode II delamination. According to our literature review, the most adopted [59, 61-70] mixed-mode interaction criterion is: 211 WI‘ OCC st co cat glu dc Cit it [EL] {21.) +[En] =Q” (605) Glc Gllc Gil/C Where the value of the parameters a and b is usually 1 or 2. If Q 21 , then delamination OCCUI'S. Kaczmarek, Wisnom, Jones [61] (1998) are among the authors computing the strain release rate implementing the VCCT method in finite element model. The computed values of the interlaminar fracture toughness Gc, critical value of G required to cause the delamination growth, has been used to study the free-edge delamination of glass/epoxy composite laminated beam loaded under bending. In 1998 Pradhan and Tay [59] applied a 3D finite element analysis to the study of delamination grth in notched composite laminates under compression. The crack closure method CCM was chosen as approach for computing the strain energy release rate G in the three delamination modes. Rinderknecht and Kroplin [63] in 1997, Wang and Raju [64] in 1996, Raju, Sistla and Krishnamurthy [65] in 1996 and Hitchings, Robinson and Javidrad [66] in 1994 developed similar finite element models for the analysis of delamination growth in composite plates. Obviously, the numerical models are different from each other, but the methodology for determining delamination initiation and growth is the same. The virtual crack closure technique (VCCT) is employed to compute the total energy release rate G, including contribution by all the three modes. 212 7.3.2 Interlaminar strength and mixed approaches In strength of material approach, local state of stress at the interface is compared with relevant strengths. Use of finite element analysis makes the strength of material approach attractive, as the stresses can be evaluated quickly and efficiently, and interlaminar stresses can be easily compared with the measured strengths. However, unless delamination initiation is the only concern, failure criteria usually combine strength of material features with fracture mechanics ones. Since the cohesive zone can still transfer load after the onset of damage, a softening model is required that describe how the stiffness is gradually reduced to zero after the interfacial stress reaches the interlaminar tensile strength. Reliable prediction of the softening behavior can be obtained relating the work of separation to the critical value of the strain energy release rate. It has been shown that if an individual interfacial location is considered, the area under the stress-relative displacement curve is equal the critical value of the strain energy release rate: G,, = f” a(6)d6 (606) where 6F is the relative displacement at failure and GC is the critical strain energy release rate for the corresponding delamination mode. Some stress-relative displacement models are illustrated in Figure 118. The material behavior is defined by the interfacial tensile strength 0",, the relative displacement at the on set of damage (30, the relative displacement at failure 6}. In single- mode delamination as the load is progressively increased, the relative displacement 6 between bottom and top FE mesh grows accordingly to the value of the penalty stiffness. When (5}) is reached the stress is at is maximum level, the interfacial tensile strength 0,. 213 For higher relative displacements the interface accumulates damages and its ability to sustain stress decrease progressively. Once 6 exceeds 6; the interface is fully debonded and it is no more able to support any stress. 0A r0 [in re 60 6”} 6” F 6 F 6 gp Figure 118. Interfacial constitutive model The delamination model just described is has been adopted by Reedy et al. [62] 1997, Davila et a1. [67] 2001, Mi et a1. [68] 1998, Chen et al. [69] 1999, and Alfano et a1. [70] 2001, Lammerant and Verpoest [71] 1996, Schellekens and Boerst [72] 1993 and Schipperen and Lingen [73] 1999, Moorthy and Reddy [74] 1999. Interface elements are introduced to connect the individual plies of a composite laminate, but the way this connection is realized can differ. Two groups can be identify: point to point interface elements acting like spring and connecting pairs of nodes, continuous interface elements connecting pairs of two or three dimensional finite elements. The work of Allix and Corigliano [75] 1999, Point and Sacco [76] 1996, Ladeveze [77] 1992, Bottega [78] 1983 clearly shares the same basic methodology, as 214 their references confirm. A constitutive law for the interface material, able to handle the delamination phenomenon, has been obtained starting from an adhesion model. It is based on the very simple physical idea of the behavior of the interface. In fact, it has been assumed that when two points on the surfaces in contact are in adhesion, their relative displacement must be zero. According to the work of Point and Sacco [75] the main feature of their delamination approach consists in the possibility of mathematically recovering fracture mechanics theory. Ladeveze [77] defined the damaged strain energy density in the interface layer as: 1 (02:12 (03212 02 02 E,,=— -,-+0 ~+ +0 32 + , (607) 2 k k(1-d,) k2(l—d,,) k1(1"'dm) where d, , (1,, and d ,,, are the internal damage variable correspondent to the associated delamination mode. Differentiation of ED with respect to each damage variable yields the associated thermodynamic forces Y, , Y,, and Y,,,. At this point a variety of damage evolution laws can be considered using a function Y (Y, , Y ,, ,Y ,,, ) , like: 1 rsr(Y/ +71YII 1'72””)2 (608) Y = sup Then, the internal damage variables are computed by means of a material function w( Y ). Delay effects are usually introduced in the function w(Y), because without their implementation the model would be accurate in the prediction of initiation but not in its evolution. The numerical procedure has been used for beam and plate problems. The results obtained using the finite element method seems to agree with the analytical solution. But, comparison with experimental results appears necessary to validate their approach. 215 The concept of interface model is also present in the work of Petrossian and Wisnom [79] 1998. They developed an interface element to model the resin-rich layers between plies, using plasticity theory. The element is able to behave elastically, yield under a changing mode ratio with either a perfectly plastic or work softening type response, and finally fail when a certain amount of plastic work has been done. According to the authors, when there is more then one possible delamination site the approach has a significant advantage over the conventional fracture mechanics methods such as the virtual crack closure technique, since it is not necessary to assume where or at what rates the cracks will propagate. Chakraborty and Pradhan [80] performed a fully 3D finite element analysis to study delamination at the interface GR/E and GL/E laminates with broken central plies. The interface between the broken and continues plies has been modeled with a resin rich layer. Their methodology for predicting delamination employs contemporarily strength and fracture mechanics approaches. The strength of material approach adopted is the most common one, the quadratic failure criterion: 67:: 2 fzr 2 273." 2 _ 171 1211171 ‘Q ‘6”) Where Z, X and Y are the transverse tensile strength, longitudinal shear strength and lateral shear strength, respectively. If Q 21, then delamination initiation occurs. This same criterion for identifying delamination initiation has been adopted by Joo and Sun [81] (1994). In this expression 5,: , I}. and "fa, are the average stress at the delamination 216 fr01 C12 01' SI dc ll front and are give by: t—fi g,QI ]=—.[:([ 0' asz ,2", (610) where x, is the critical distance over which the interlaminar stresses are averaged. The recommended value of x, is twice the ply thickness. As fracture mechanics approach the crack closure method CCM is used. The delamination is assumed to takes place if Q 21 or if the strain energy release rate G=Gc. Effects of important parameters like ply thickness, fibre angle, location of the break and stiffness of the interfacial resin has been studied. Delamination evolution is not explored. Ko, Lin and Chin [82] (1992), S. Mohammadi, Owen and Peric [83] (1998), Zhao, Hoa, Xiao and Hanna [84] (1999) and Chang and Springer [85] (1986) predicted delamination initiation by a strength approach very similar to the one of Chakraborty and Pradhan [80]. Here, the mathematical form changes to: _ 2 _ 2 _ 2 [fi {122] + E: zez (611) Z X Y If e 2 1, then delamination initiation occurs. In [82, 85] delamination evolution is not explored, while in [83] the crack growth is studied with the same methodology like in [62]. A three-dimensional finite-deformation cohesive element and a class of irreversible cohesive laws was developed by M Ortiz and A. Pandolfi [86] (1999) and A. Pandolfi, P. Krysl and M Ortiz [87] (1999). The element has some unique characteristics, 217 but the interfacial constitutive model is of the familiar form summarized in Figure 118. Moes, Dolbow and Belytschko [88] (1999) presented an innovative finite element able to account for crack growth inside the element. The numerical model is quite complicate and will not be described here. The approach adopted for delamination initiation and growth is based upon the maxrmumcrrcumferential stresscntenon Stresses. 6666611111166 by the classical fracture mechanics equation ”for the stress distribution at the crack tip. This finite element type of approach, also referred as singularity element approach, was also used by MA. Aminpour and K. A. Hosapple [89] (1991) and Jin Lee and Huajian Gao [90] (1995). 218 CHAPTER 8 A NOVEL INTERFACE TECHNOLOGY FOR MODELING DELAMINATION The interface element developed in the present study can be readily employed for the analysis of delamination growth in composite laminates. In addition, the present interface element approach would have several advantages over the conventional FE one. A special feature, which is useful in simulating delamination, is the ability of releasing desired portions of the interface smaller than the finite element length, determining the crack advance. This can be achieved by changing the extreme values of the interval of integration of the interface element or by reducing the value of the penalty parameter for that interval. Also, within each interface sub-region it is possible to evaluate forces at the interface and to reduce the value of the penalty parameter as needed. Thus, it is possible to overcome the limitation common to the delamination techniques found in the literature, which requires delamination growth to be simulated in a discretized form by releasing nodes of the FE model. Moreover, the ability of estimating the normal force at the interface, allows for the development of an accurate and mesh independent friction model. 219 8. 1 Damage Model The damage model implemented in the developed interface element is one that is commonly adopted [62, 67-74]. It mixes features of strength of materials approaches and fracture mechanics. A bilinear softening model has been implemented in the present model. see Figure 119. O Figure 119. Bilinear interfacial constitutive damage model. In single-mode delamination, as the load is progressively increased, the relative displacement 6 between the bottom and top FE mesh grows proportionally according to the value of the penalty stiffness y. When 60 is reached the stress is equal to the interfacial tensile strength 0",, the maximum stress level possible. For higher relative displacements the interface accumulates damage and its ability to sustain stress decreases progressively. Once 6 exceeds 6; the interface is fully debonded and it is no longer able to support any stress. If the load were removed after 60 has been exceeded but before 6; has been reached, the model would unload to the origin. For example, if after reaching point K the load is reduced, the model unloads along the line K0. If the load is reapplied, the stress 220 gn 1111 3C grows with the relative displacement along the same line K0. This behavior is obtained by an effective reduction in the penalty stiffness y. In the present model, a new parameter D is introduced in order to signify the damage accumulated at the interface: 0' =(1 -D)76 (612) Thus D is a damage parameter, whose initial value is zero. D starts growing when 6 2 60 and reaches the value 1 when 6 2 6p. The value of D is computed from geometry to be: = 6,,(6 —60) 6(6, _50) (613) 13(5) The interfacial constitutive model is entirely defined when two of the following properties are known: G,, a}, 60 and 6p. The following tWo relations exist among these parameters: G, = 530’ (614) a, = 94 (615) r The bilinear interface model is applied to a sub-region of the interface between the two meshes; the smaller is the length of each sub-region the higher is the accuracy of the prediction. A conventional implementation of the discussed damage technique requires the model to be applied along the length of one finite element, wherein the crack can advance in a discrete way only by failing one element at a time. Both limitations necessitate the use of a refined finite element mesh. Thanks to the special features of the 221 interface element previously presented, the damage evolution scheme in our model is effectively mesh-independent, wherein sub-regions of the interface much smaller than the finite element length can be released. Thus, the bilinear interface model can be applied to a desired fraction of the interface between the two meshes. If we divide the interface element into a given number of intervals n, each of them will obey the rules of the failure model independently from the others. This corresponds to changing the total potential energy of the system in the following way: 1 ” L, 2 ” L, 2 7: = 71"), +7rQ2 +§§(1_Di)7l iL,_,(V_ul) ds+,§,(l—Di)7ZIL1—1(V_u2) ds (616) where D,- is the damage parameter associated with the interval i, and the interval i is defined over the range ( L,.,, L,- ). L, is the value of the interface coordinate L at the end of the ith interval. The value of the relative displacement 6 is evaluated at the center of the interval i. By allowing crack advance in a more continuous way, a higher accuracy of the simulation can be achieved. To illustrate this concept, consider the two incompatible meshes shown in Figure 120. They are joined by two interface elements whose length equals that of five conventional elements of the top mesh and four elements of the bottom mesh. Two forces F applied at the tip are responsible for affecting a mode I stress field at the interface. The interface element at the crack tip is shown in Figure 120(a) as divided in 10 intervals. The intervals don’t need to match any of the nodes in the upper or lower mesh, but in this example some of them coincide. Moreover, the number of intervals in the interface 222 element is a parameter that can be changed as desired. 11 lJllllTld l..— ._..0.—.l _., [570% fl I 1 1 1 1 1 1 (b) ‘ 1 1 1 1 1 1 1 1 1 1i! (1 1 1 1 1 1 1 1 1 . ‘ L 1 1 1 F 1 1 r “ % £ 9 (d) 1 1 1 1 1 1 1 ‘. fl‘ 1 1 L 1 1 1 1 1 ’ 1.) 1 1 1 1 1 1 ~ .. Figure 120. Division of the interface element in intervals Following the progression in Figure 120, a simulation of the delamination growth is achieved by releasing portions (intervals) of the interface element. In Figure 120(b), the 223 first interval is failed. The portion of the interface element not released still applies its constraint to the lower and upper element next to the crack tip. In Figure 120(c) the second interval is failed, this determines the complete release of the element on the upper mesh near the crack tip. The element on the lower mesh moves downward too, but being still held in part, the movement is small. The next advance in the crack length, Figure 120(d), frees the lower element and only partially the upper. Figure 120(e) shows the effect of the releasing yet another interval. Similarly, when all the intervals are failed the FE model behaves as though the first interface element is not present. Then, the next interface element starts failing. In this example, intervals whose length is half of the smaller finite element extension have been used. Dividing the interface element into more intervals or reducing its length would obviously improve the accuracy of the model. 8.2 Mixed Mode Approach A novel approach for correctly simulating the delamination growth in mixed mode I+II conditions has been developed. Due to the peculiar finite element formulation of the penalty interface element, some degree of creativity was required. The damage models for the two modes needed to be linked according the quadratic failure criterion and the quadratic interaction criterion, maintaining at the same time independent damage parameters. The properties required to define the softening model are the following: the interlaminar tensile and shear strengths T and S, the penalty parameter y, the critical strain energy release rate G,,, and G,,,. At a given increment the finite element solution allows the computation of (Z and 6,, for an interval of the interface element; where 6: and 224 f0] 11 6, indicate, respectively, the mode I and mode II relative displacements. Using the following relations among stresses and relative displacements: 6... =3 (617) 7: 6.153 (618) 7:: 0. =26. (619) T...- = 7.5.. (620) 0': 2 TX: 2 _ 171151 -1 (.21, 2 2 [i] +[6‘ ] =1 (622) 6:0 6x0 , If the condition in (622) is not satisfied, no action needs to be taken. Otherwise, we assume the following ratio: at. Iual / a = C, (623) or tual —-°-l— \ 6:0 j K ( to be the same it was when the failure condition (622) was verified first, the point of damage onset. Figure 121 will help in clarifying the idea. At the increments when the left hand side of equation (622) is bigger than 1, the relative displacements 6, and 6; are higher than the minimum values necessary for damage onset. In Figure 121, where 63./6,0 and 6760 are the main axes, such a configuration could be represented by a generic point 225 A outside the quarter of circle defined by the quadratic failure criterion. 6/5x0 A (6/6x0)actual (ti/53.0) ’ 5z/5w > 0°“ (5/620)’ («z/6.01““ Figure 121. Graph 6/6w versus 61./do. We assume that when (622) was satisfied first the ratio (623) was the same. If the load steps are small, the assumption holds. Then, from geometry, it is possible to determine the value of the relative displacements 6; and 6; in F (here the superscript actual is omitted). 3;] 6 6; =6... (624) 226 i 6: i 6 a: = 5- 3° (625) The interfacial constitutive models need to be updated for both delamination modes I and II setting 6:0 = 6; and 65,, = 6; , see Figures 122 and 123. 0' 1 r T ‘ \ \ \ ‘ \ \ -, I I T! 5,. 61.; .9» t--------------- ----- V O °2 11$ Q .91 1., Figure 122. Updated interfacial constitutive model for mode I delamination 0' 41 IT‘\ I. ‘\ I ' \‘ l \ r : ‘ I s S' I I ‘s I ‘s s s I ‘s I ‘s I ‘\ I ‘s I ‘s | I \ I I s x s 7’ : ,\ I \\ 1 ~‘ 6 I s l I ‘~ ‘ I ‘s‘ l s A O ' I I 5.0 510 (21: 5.F Figure 123. Updated interfacial constitutive model for mode 11 delamination 227 The interlayer tensile strengths T and S are updated accordingly: 7" = y: .520 (626) S' = 7, ~6',0 (627) Note that the following inequalities hold. 6:0 S 6,0 , 6:0 5 6,0 , 6, 2 6:0 , 6, 2 6:0 (628) However, the modifications to the interfacial constitutive models do not end with change in the relative displacements at onset of damage. The quadratic interaction criterion predicts final failure to be reached when the following condition is verified: 2 2 [i] 4.02] :1 (629) 011' Gllc In a similar way as we did with the failure criterion, we assume the ratio between (G,,/G,,,Jmual and (GI/G,,)”mm not to vary as the work of separation growth, Figure 124. The point representing the actual situation of the system A, calculated using the strain energy release rate at the current increment, will be inside the domain defined by the quadratic interaction criterion. If not, the interval is totally failed. We assume that when (629) is satisfied the following ratio: = C, (630) will be unchanged. Research on the specimens commonly used for delamination studies shows that for a given configuration (geometry and loads) the ratio between the strain energy release rates for modes I and II, Gl/G”, do not change much during the entire test [56]. This fact provides a valid foundation to our assumption, in particular since the 228 EC 1h 10 equation (630) must hold only for the very limited advancement of the crack related to the failure of the portion of the interface element under consideration. Thus, it is possible to determine the value of the strain energy release rate G,’ and G,', in F. > GIl/GIIc 1 (GI/GM!) , (GI/Gurm “-0 (G/GIJ‘W’ (cl/Gr): Figure 124. Graph G,,/G,,, versus G,,/G,,. From geometry it is possible to determine the value of the (G,/G,c)' in F (here the 191 G,,. 2 2 (631) 11°] 1°] __ + _. G,,. G,,, In order to evaluate this expression we need to determine (G,/G,C), (actual value of G, superscript actual is omitted). Gl' = Glc divided by G,,), see Figure 125. 229 . T'-6.', 637963, oh: 2.} ___ -0 2. -l 5:0 .y: 53'} _ [62(1_D:)7:].65F G, = G,,. - (Area Triangle OBK) = 2 §:O__}:Z._§_;li__ [62(1_Dz)7:]§~:i__ —G’ T 2 2 6- = p ' = 1— —‘ 1 _ D, [le M 63,, ( *) 2 With the same methodology, G,', can be computed. (632) (633) (634) Now, the changes to the interfacial constitutive models can be completed for both delamination modes I and II setting G,',. = G,’ and G,',, =G;,. Figure 125 and 126 shows the final form of the interfacial constitutive models. As it can be noticed, the models have different penalty and damage parameters. In such a way, the maximum freedom is allowed for modeling modes I and II delamination. 0'. A T’ B K 3 ,2 G ,c ,,,, ',.," vvvvv 1(1-0322 L—"lig 1 6:2 .4" 4 O , . 620 2F Figure 125. Final interfacial constitutive model for mode I delamination 230 1?‘ > " O r ‘l I O 620 6,xF Figure 126. Final interfacial constitutive model for mode 11 delamination 8. 3 Energetic Approach for a Correct Distribution of Stress Initial numerical testing on the developed damage technique produced very good results for mode 11 delamination, but prediction of mode I behavior was sometimes inaccurate. In order to identify the causes of the problem, linear elastic fracture mechanics prediction of the stress distribution near the crack tip was compared to the one obtained with the interface model [49, 50, 91, 92]. Geometry and loads for the mode I case are shown in Figure 127(a), while original and deformed meshes of the finite element model are shown in Figure 127(b). The material has Young’s modulus 1 MPa and the thickness is 1 mm. Two meshes compose the finite element model; they are joined by one interface elements of 16 nodes along the expected crack growth path. Clearly the two meshes are compatible in this case. Normal stress G,, is obtained along the interface taking advantage of the peculiar abilities of the interface element in recovering the forces at the interface. 231 Interface element l l l l b=10 ‘ (a) (b) Figure 127. Fracture mechanics model. ”l -y 1 F—Oi—Fraictuire Mechanics Solution 1 1 El - D- - Interface Element 1‘ i — 9 — Improved Interface Element I I I I r I *I , 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6,0 6.5 7.0 7.5 Figure 128. Normal stress distribution along the entire interface. Figure 128 shows the stress distribution along the entire interface. Stress values in the two graphs are reported in a discrete form, they are averaged at the center of the 75 intervals in which the interface element is divided. 232 The fracture mechanics solution is computed according to the following 7217511+Sini§lsmiigiicosigi (635) K = f,(2.5,10).,/(2.57z) (636) equations: 0",,(K,6,s) = 2 f1(a,b)=1.122—O.231[%] 10.55%?) 3 , (637) —27.710[3] +30.382(£) b b f’a,,(1<,0.s)dr a',(x,,x,)= " x -x (638) 2 l Where 0', is the average value of the stress 03,. inside one interval. If we observe the stress distribution predicted by the Interface Element model, it is easy to notice an abnormal behavior near the crack tip. This behavior can be explained if we remember that in the adopted penalty method the continuity constraint is not enforced exactly. Instead, the integral of the square of the relative displacement must be a very small number. The inaccuracy observed in Figure 128 occurs only when very strong gradient of stresses are present at the interface. Mesh refinement would improve and eventually correct the results, but it is preferred that the model behave correctly for every possible discretization of the domain, especially since the accuracy of the delamination model strongly depends on a precise calculation of the stress and displacement distributions at the crack tip. The approach for overcoming the problem came from energetic consideration: it is possible to relate the fracture mechanics stress distribution to the interface model 233 solution by imposing that they must contribute equally to the total potential energy of the system. The energy associated with the continuity constraint imposed by the penalty parameter, can be expressed also in the following form: _ 1 ‘ 2 15(5) .. 3y I 5 ds (639) where 6 is the relative displacement of the connected interfaces and is a function of the position along the interface 6 = [u2(s) - u,(s)] = 6(5). Fracture mechanics predicts, both for Mode I and Mode II opening, the stress to vary near the crack tip according to an expression of the type: 0' = 0(3, K ,6) (640) For a given geometry and load applied to the system K does not change. Moreover, along the interface 6: 0. It follows: 3 0(3) = J; For an isotropic material, where s is the distance from the crack tip and C, is a constant. (641) In the penalty interface element approach, the stress at any location is a direct function of the penalty parameter and the relative displacement 6: (IQ/,6) =C2 ~76 (642) where C; is a constant. Now, we want replace the 6(s) distribution obtained from the finite element model with a new one, 6 ’(s), such that: E(6) = E(6') (643) If we choose 6 ’(s) to be: 234 5 '(s) = —-‘ (644) h It follows: 0(3) = 5 (645) J; This means that the stress would vary along the interface according to the fracture mechanics predictions. The constant C 3, in the expression for 6’(s), is computed by imposing the energetic constraint: E (6 ) = E (6 ') . 15(5):? j“ ’(5')2 ds (646) E(6)-l [’gz—ds—l C2[lns -lns] (647) —27<1 S 5‘27 3 2 1 E(6) is a known quantity that is computed once 6(s) is obtained from the finite element solution. The resulting improved stress distribution is plotted in Figures 128. This energetic approach for obtaining a correct distribution of stress near the crack tip is used only when the abnormal behavior is present and only for mode I delamination. An algorithm has been implemented in the code able to detect when a wrong stress distribution is present and automatically switch from the normal approach to the modified approach. In order to identify possible limits of the approach in its applicability to composite materials, the literature on the stress distribution near the crack tip when the materials of the connected meshes are different has been explored [93-97]. It has been found [97] that the change in the stress intensity factor for mode I with a variation of the Young’s modulus of the two materials is very limited. Even with a ratio 235 E ,/E2=70, where 1 and 2 marks the two different materials, the change in K, is less than 10%. In 2001 Alfano and Crisfield [70] have done a parametric study on the influence of the tensile strength 0', of the interface on the accuracy of the delamination growth prediction. The damage model they implemented in their interface element is the same as the one used by the present approach, the bilinear softening model previously described. They found that as long as the softening behavior is obtained relating the work of separation to the critical value of the energy release rate, moderate changes in a} don’t affect the results. Finally, a parametric study on our interface element has been performed. The exponent ,6 of the s, distance from the crack tip, in equations (641), (644) and (645) has been varied. The simulation of delamination growth for a DCB specimen, whose experimental results are in good agreement with our model for ,6 = -1/2, has been done with ,6 = ~1/4 and ,6 = -3/4 (according to [94], 0<,6 <-1). The load-displacement curves obtained with the different values of ,6, didn’t differ from the outcome for ,6 = -1/2. The limited influence of the ratio E,/Ez on the stress intensity factor for mode 1, plus the lack of need in getting a completely exact stress distribution, plus the scarce variation of the results with ,6; lead us to conclude that in the rare cases when the abnormal stress distribution is present and the materials at the interface are different, the proposed energetic approach should still be able to yield correct results. 236 8.4 Friction Model 8.4.1 Evaluation of the force at the interface For the 1-D case, Timoshenko beam element, a known relation between the Lagrange multipliers and the penalty parameters can be used to evaluate force at the interface. If the constraint C, that has been enforced is: C, (u,,u2,v)=(v—u,) (648) Defining u, ,,u,,,v, as the solutions derived from the penalty formulation, it is possible to verify that the Lagrange multiplier is related to the penalty parameter in the following way: zl, =7,-C,(u,,,u2,,v,)=y,-(v,—u,,) (649) Since the values of the Lagrange multipliers correspond to the force required to hold the two nodes together, the force at the interface can be easily calculated: F, =7, -C, (u,,,v,)= 7, ~(v, —u,,) For 2-dimensional elements, like 2D quadrilaterals and plates, the interface force is still simple to compute, but the relation between Lagrange multipliers and the penalty parameters need to be updated. In this case the term added to the Total Potential Energy of the system is: %y!(v—u)2ds (650) Where u and v are not nodal displacements, but displacements fields functions of s. As a consequence of this difference the Lagrange multiplier is now related to the penalty parameter by: 237 F=2=y-j(v—u)ds (651) A very simple verification of this relation is given. A 2D problem, analogous to that of section 3.1, has been analyzed in section 4.4.1. A rectangular plane stress element clamped to a fixed penalty frame along one side is loaded on the other with a uniform axial load F, Figure 17. The interface points V,, V,, V3 are fixed, so their displacements V,,V2,V3 are equal to zero. The solution derived by the finite element analysis for the displacements of the nodes of the element on the interface is: “99 ”11: b 7, (652) The integration required to compute the desired relation is expressed by: f(v—u)dy= h f{[T11(y)v. +T21(y)v2 +T3' (y)v3]—[N, (y)u9 +N,(y)u,,]}dy+ (653) i: {1712 Mm +732 (y)v2 +T32(y)v3]-[N1(y)u9 +N2 (y)un]}dy ‘- Where T ,k are the natural cubic spline interpolations functions of section 4.2.2, with the index k marking the first and second groups of functions valid, respectively, for the first and second intervals. The result of this integration is expressed by: 3bv, + 5bv2 + 3bv3 bu9 _ bu,, — (654) 16 8 l6 2 2 Substituting the achieved results in the relation for 11 it follows: b F b F ,1: - - d = - ——————— — =—F 655 71(V")8712b72brl () As predicted the correct value of the interface force is obtained. It is important to notice that we are not restricted to compute the force on the entire interface, but it can be evaluated for a portion of the interface length too. This can be performed by changing the extreme values of the interval of integration. Finally, computation of the interface force does not depend on the compatibility of the interface meshes. Even if the discretizations of the adjacent layers, or layer and skin-stiffener, are different we are still able to choose liberally the interface length along which to evaluate forces. 8.4.2 Implementation of a Friction Model The friction model can be applied to interface elements after complete failure or to interface elements whose only purpose is to avoid overlapping and enforce friction. For each of the intervals in which the interface is divided, the normal force at the interface F ,, can be calculated given the normal relative displacement 6,, : 1 F. = 57. 16.66 (656) The tangential force F, that the interface element needs to generate to simulate the friction phenomenon is: F. =%r.(1-D.) [6.61s =#F.. (657) Where ,u is the friction coefficient and 6, is the tangential relative displacement. The parameter D, that before failure was used as a damage parameter, is now employed as a scale factor able to reduce the value of the penalty parameter for the tangential DOF. 239 Assuming the value of D, that makes the equality (657) to hold true, the right amount of fi'iction will be generated by the interface element. The required damage parameter D,‘ related to the tangential relative displacement 6, can be determined from the following relation: D’—1- ”F" ' — 7.15.615 (658) 8.4.3 Numerical Test for the Friction Model A simple test has been performed to verify the reliability and the accuracy of the friction model implemented into the interface element. ov = l.E+6 Pa WIS. A\V w L = l. m Interface Element Figure 129. Boundary conditions and the geometry of the test for the friction model The loading, the boundary conditions and the geometry for the test problem are illustrated in Figure 129. Two finite element meshes are connected along the common boundary using one interface element. All the finite elements are rigid. The drawing is 240 three-dimensional for clarity, but the FE analysis was performed with a 2D mesh having the indicated thickness. During a quasi-static analysis of 400 load increments, a 1mm tangential displacement is gradually imposed to the nodes of the upper mesh. At the same time, a pressure load of 1 MPa is applied to the top the upper mesh and is kept constant through the entire analysis. According to the Coulomb law F, = ,u-N, since a coefficient of friction equal to ,u = 0.1 is assumed and the contact area is 1 m2, as outcome is expected a total reaction force at the interface nodes of: F, =,uF,,=/10'A=0.1-l.e6—N7-1.m2=l.e5N (659) m In Figure 130 is reported the reaction forces versus the tangential displacement produced by the FE analysis for the interface nodes 1, 2 and 3. #75000; *5 I I ' A I 4 ~ g ------- Node 2 ‘g , —-—--Node3 : 50000 -, ---------------------------------------------------------------- l '8 i 1 3 : °= 2 25000 i, 1 Horizontal displacement (mm) , O r r 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0l Figure 130. Test Friction — Graph F owe-Displacements. For each of the 400 increments in the displacement value, the reaction force at the central node is recorded to be a constant value of 0.5e5 N while on the other two nodes a 241 constant value of 0.25e5 N is obtained. The outcomes of this simple test, being accurate and constant through time, confirm the robustness of the adopted friction model. 8. 5 Numerical Results Results for two different double cantilever beam DCB specimens, one end-loaded split (ELS) specimen and one Fixed-Ratio Mixed Mode (FRMM) specimen are presented. These results are compared to measured data to assess the ability of the present damage model to simulate delamination growth. 8.5.1 Double Cantilever Beam Test #1 The loading, the boundary conditions and the geometry for the double cantilever beam DCB specimen are illustrated in Figure 131. The DCB test is recognized as an accurate pure mode I test. The properties assumed for the beam material are: E, ,=130 Gpa, E23: E 33 = 8 Gpa, G,, = 6 Gpa, V: 0.27. The properties of the DCB specimen interface are: G,, = 257 N/m and 0', = 48 Mpa. 22 8’ 11 L=150 V Figure 131 . Geometry and boundary conditions for the DCB test specimen. 242 DISPLACEMENT IAGNIPICATION FACTOR - 6. 67 DISPLACED HESH RESTART FILE - tbean STEP 1 INCREIENT 400 TIME COMPLETED IN THIS STEP 1.00 TOTAL ACCUKULATED TIME 1.00 ABAQUS VERSION: 5 8-1 DATE: 30-APR-2001 TIME: 113-01:51 Figure 132. Original and deformed (5mm opening) models of the DCB test specimen using a 120x2 mesh. Three finite element models of varying mesh refinement have been used to verify the capabilities of the present approach. In Figure 132, original and deformed models of the DCB test specimen are shown for the coarsest mesh. Two independent meshes compose the finite element models of the upper and lower part of the beam; they are joined by several interface elements. Each interface element connects eight finite elements: four on the bottom mesh and four on the top. Experimental results used to validate the present delamination approach have been reported by Chen and colleagues [69]. A plot of the reaction force as a function of the applied end displacement is shown in Figure 133. Results from the experiment are compared to analytical ones for all the finite element models: a) mesh 120x2, b) mesh 300x8, c) mesh 600x8. Note that the total number of elements in the meshes is indicated, for example two 120x1 meshes joined by interface elements compose the 120x2 mesh. The outcomes from 120x2 mesh present many local “bumps”. Mi, Crisfield, Davies and Hellweg [68] have analyzed this phenomenon, concluding that coarse meshes can induce these “false 243 instabilities”. As the mesh is refined, the predicted response agrees very well with that measured experimentally. As discussed previously, the damage technique implemented in our model allows portions of the interface, intervals, much smaller than the finite element length to be released. In Figure 134 convergence of the solution with the number of intervals is investigated for the 300x8 mesh. It can be noticed that as the number of intervals increases the accuracy of the results raise. Figure 135 illustrates the model behavior as function of the number of increments for the 300x8 mesh. Convergence of the solution to the experimental one increasing number of intervals is achieved. 14o. * ’ * - 1 - L ,- - ------ 120x2 mesh —\ E — - - 300x8 mesh w 120‘ § 600x8 mesh 1‘. . . . = 100« . .' -. .. - ' Expenmeflta' «9 '- ' ' ‘s - 1 ' = . 60- l 40- l 204 Opening displacement (mm) o l’ U V V I i 0 1 2 3 4 5 , Figure 133. F orce vs. displacement results of experimental and numerical DCB test 244 . «7: 9:3“ Etmxhtzk Fig 120; ' a 5“ * . ------ 2lntervals 3‘ E. 3 — — - 4 Intervals , L ,2 10,, ,. , 8 Intervals ’ s 0 1 ~ g a - 0 __ [Expatriate 1 E , l 1 N 301 ‘ g . l 60 - iiiiiiii l 40 . , l 20 1 Opening displacement (mm) , o r I r I 1 , 0 1 2 3 4 5 ' 1 Figure 134. F orce vs. displacement results of DCB test. Convergence of the solution with intervals number. Mesh 3 00x8. 160 - 140 « ' 120 . . e 100 “—1 a H 1 80-1 60‘ 40‘ 20j , Log (N. Increments) 0 u w 10 100 10001 Maximum reaction force (N) Figure 135. Convergence of the solution with the number of increments. Mesh 3 00x8. 245 8.5.1 RIO V perf1 Fig 108 the 38.) M cl 8.5.2 Double Cantilever Beam Test #2 To further verify the ability of the interface in accurately simulate delamination growth in composite materials, a second double cantilever beam DCB test was performed. L=100 ,1- V Figure 136. Geometry and boundary conditions for the DCB test #2 Material data and experimental results are taken from a different source [70] with respect to the first DCB test. The loading, the boundary conditions and the geometry for the double cantilever beam DCB specimen are illustrated in Figure 136. The DCB test is recognized as an accurate pure mode I test. The properties assumed for the beam material are: E,,=126 Gpa, E22=E33=7.5 Gpa, G,2=4.981 Gpa, v=0.281. The properties of the DCB specimen interface are: G,, = 263 N/m and 0', = 57 Mpa. A plot of the reaction force as a function of the applied end displacement is shown in Figure 137. A comparison of the experimental results with those from the finite element simulation further asses the robustness of the implemented delamination model. 246 F 1' 11 ”"2 120. — - -- -- — -~ 5.1 0 Experimental . 1 E 1 * M 100 , 1 ——Interface Element ' , u ,......1 l 3 1 1 Lg, , g 30 l O~ J . 1 I“) :1 N a: 60- . l . l 40- 20« Opening displacement (mm) 1 0 I I r r T I I TI I l 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 , 7 i V i 7 7, 7A 7 A v + L? L__ v i. _‘s v w ___.— _—7 .L__v._ ,..T_a Figure 137. Force versus displacement results of experimental and numerical DCB test #2 8.5.3 End-Loaded Split (ELS) Beam Test To characterize the mode 11 delamination, shearing mode, the end-loaded split (ELS) specimen is often used. The loading, the boundary conditions and the geometry for the end-loaded split (ELS) specimen are illustrated in Figure 138. The properties assumed for the beam material are: E, ,=l30 Gpa, E22: E33 = 8 Gpa, G, 2 = 6 Gpa, v: 0.27. The properties of the ELS specimen interface are: G,,. = 856 N/m and a", = 48 Mpa. As in the DCB test specimen models, two meshes compose the finite element models of the ELS specimen; they are joined by several interface elements. For the results shown, a 300x8 finite element mesh has been utilized. 247 h = 3.05 b=2.4 L=105 ‘5 Figure 138. Geometry and boundary conditions for the ELS test specimen. DISPLACEMENT HAGNIFICATION FACTOR - 0.389 DISPLACED MESH RESTART FILE I theal STEP 1 INCRENENT 300 TINE COMPLETED IN THIS STEP 1.00 TOTAL ACCUMULATED TIME 1.00 ABAQUS VERSION- 5 8—1 DATE: ll-NAY—ZUOI TIME: 12:55'01 Figure 139. Original and deformed (30mm tip displacement) models of the ELS test specimen for a 300x8 mesh. In Figure 139, original and deformed models of the ELS test specimen are shown. Experimental results used to validate the present delamination approach have been reported by Chen and colleagues [69]. A plot of the reaction force as function of the applied end displacement is shown in Figure 140. Convergence of the solution as the number of increments increases is demonstrated. These results indicate that the model is reliable and accurate. 248 2k tttks‘k F: «>2 Fig th ft Si 1601 7' H i if ii iii 2 l ------ 200 increments i : 14ol ‘ - - - 400 increments ’ ' - U ' \ , '5 | - - - - 800 Increments \\ . k. 120 1000 increments . \ . 'l . 1 g . of Experimental ,. \ .\\ ’ l ‘6 ‘\ . \ l g 100 ' .~ \ i a: 0 l 80 - l l 60 - l 40 - i 20 -1 Opening displacement (m) l. 0 r U I T V U i o 5 1o 15 20 25 30 i Figure 140. Force vs. displacement results of experimental and analytical ELS test for varying number of loading increments. 8.5.4 End-Loaded Split (ELS) Beam Test — With Friction Results from the end-loaded split (ELS) specimen test are satisfactory; however the model was not as accurate as in the case of the double cantilever beam DCB tests. This reduction in the quality of the prediction may be related to the absence of friction forces in the simulation. That is the reason why further testing as been done on the ELS specimen test taking into consideration friction forces. The loading, the boundary conditions, the geometry and the material properties are the same as in the previous ELS specimen test, the presence of friction being the only difference. initial FE analyses revealed a limited influence of the friction model on the results. A possible reason was identified in the presence of areas of highly concentrated 249 ~--il I r7 yet small normal forces at the end of beam where the load is applied and near the crack tip, while in between the normal force was zero. In these conditions the friction model might have been unable to produce the correct amount of friction force. Even if it is hard to say with certainty if the friction forces are really too small in this case to affect the results or if the model couldn’t perform well under such conditions, we decided to try to apply a small pressure load on the upper part of the specimen in order to produce a larger contact area. The pressure load is applied approximately in the region indicated in Figure 141, and it is chosen very small, 6 N, so that it shouldn’t affect directly the linear portion of the reaction force versus displacement graph. Figure 141. Original and deformed (3 0mm tip displacement) models of the ELS test specimen for a 3 00x8 mesh. Figure 142 shows that in the presence of friction forces, the prediction of the ELS specimen behavior agrees more closely to that observed in the experiments than when friction is ignored. 250 160; ._J N E, " 0 Q '~ 140 « / ’ \ e LL. / a / °§ 120 . /’ U / Q / \ / \r / 100 - / \ I z / / ‘0— / 80 - / 60 - ’ : ’ Interface Element - no frictlon 40 . . - - - Interface Element - with friction 0 Experimental 20 « ’ Opening displacement (mm) 0 I I 1 fi I I 3 o 5 1o 15 20 25 30 , Figure 142. Force vs. displacement results of experimental and numerical ELS test in presence of friction. 8.5.5 Fixed-Ratio Mixed Mode (FRMM) Test To test the mixed mode I+II delamination capabilities of the interface model, the Fixed-Ratio Mixed Mode (FRMM) specimen has been considered. The loading, the boundary conditions and the geometry for the FRMM specimen are illustrated in Figure 143. The properties assumed for the beam material are: E 1 ,=130 Gpa, E 22: E33 = 8 Gpa, G ,3 = 6 Gpa, v: 0.27. The properties of the ELS specimen interface are: G16 = 257N/m, 0”,. = 856N/m and 0', = 48 Mpa. As in the DCB and ELS test specimens models, two meshes compose the finite element models of the ELS specimen; they are joined by several interface elements. For the results shown, a 300x8 finite element mesh has been utilized. 251 r L=105 T . l V..- Figure 143. Geometry and boundary conditions for the FRMM test L. DISPLACEXENT KAONIEICATION PACTDR - 0.583 DISPLACED MESH RESTART FILE - than: STEP 1 INCBEIENT 1 TIME COMPLETED IN THIS STEP 2.220E-16 TOTAL ACCUHULATED TIME 0. ABAQUS VERSION- 5.8-1 DATE: 19-NOV-2001 TIRE: 22:57:12 Figure 144. Original and deformed (20mm tip displacement in absence of delamination) models of the F RMM test specimen. In Figure 144, original and deformed models of the ELS test specimen are shown. Experimental results used to validate the present delamination approach have been reported by Chen and colleagues [69]. A plot of the reaction force as function of the applied end displacement is shown in Figure 145. Though the failure load is underpredicted, the overall results fit well with the experimental results, proving that the novel mixed mode damage model implemented in the interface element is able to correctly predict delamination growth in mixed mode FE simulation. 252 :<\ =p~l.\.‘ --\.~‘.\-\§Q Fig 14o - a ,/ Interface Element . N I. . . i g 120 , ll - - - Analytical Linear Part ; It. / - g 0 Expenmental "5 100 . 7 i A ' i 8 E 80 - . 50 4 . . O 40 1 20 « Opening displacement (mm) ‘ 0 I' T I I . o 5 10 15 7 20 E Figure 145. Force vs. displacement results of experimental and numerical F RMM test. 253 CHAPTER 9 CONCLUSIONS In the present work, an effective and robust interface element technology has been developed using the penalty method. This approach overcomes the numerical difficulties associated with the existing methods based on Lagrange multipliers. Additionally, the present formulation leads to a computational approach that is very efficient and completely compatible with existing commercial software. Significant effort has been directed toward identifying those model characteristics (element geometric properties, material properties and loads) that most strongly affect the required penalty parameter, and subsequently to developing simple “formulae” for automatically calculating the proper penalty parameter for each interface constraint. A wide variety of problems has been investigated analytically and computationally in order to correctly identify the analytical functional form of the penalty parameter for each constraint within many classes of problems. The resulting interface element is particularly efficient for finite element modeling of composite structures. When the material properties vary across and/or along the interface, the present method is often much more accurate than adopting a unique value of the penalty parameter. The present interface element has been implemented in the commercial finite element code ABAQUS as a User Element Subroutine (UEL), making it convenient to test the behavior and accuracy of the interface element for a wide range of problems. 254 Additional capabilities were implemented in the interface element in order to simulate delamination growth in composite laminates. This entirely new technology is capable of joining and simulating crack growth between independently modeled finite element subdomains (e.g., composite plies). The interface element approach makes it possible to release sub-regions of the interface surface whose length is smaller than that of the finite elements, thereby allowing for a mesh-independent tracking of the crack front. A novel approach for correctly simulating the delamination growth in mixed mode l+II conditions has been developed. It successfully integrates widely used failure criteria in the peculiar penalty-based formulation of the interface element. A mesh-independent friction model has also been created; it can be applied to interface elements after complete failure or to interface elements whose only purpose is to avoid overlapping and enforce friction. Results for double cantilever beam DCB, end-loaded split (ELS) and fixed-ratio mixed mode (FRMM) specimens indicate that the method is capable of accurately predicting delamination growth for Mode 1, Mode II and Mixed Mode I+II. 255 APPENDICES 256 APPENDIX A USER MANUAL 257 Uelglue.f for 2-D Quadrilateral Elements The FORTRAN program uelglue.f is an ABAQUS User Element Subroutine (UEL) that calculates the finite element matrices associated with a penalty-based interface element able to connect incompatible meshes of 2-D quadrilateral linear elements. All of the necessary material and geometric data for the two connected meshes is located in the ABAQUS model input file (.inp). The required parameters are the nodes on the interface and the thickness and Young’s modulus of the two domains. Each penalty interface element is composed of a user specified number of equally spaced nodes. A patch test analysis is used to illustrate the basic format of the input file. Patch Test A two dimensional rectangular region is assumed to be composed of two domains. The domains are meshed independently and joined by one interface element. The geometrical configuration of the problem is shown in Figure 19. Other properties are: El = E2 :1 and thickness =1. In discretizing the two subdomains, four-noded bilinear plane stress elements are used. Three different conditions are investigated: uniform axial strain in direction 1, uniform axial strain in direction 2 and uniform shear strain. This problem was intended to serve as a patch test for the interface element. If the continuity constraint along the inclined interface is not correctly imposed by the interface model, then discontinuity of displacements and strains across the interface would result. The condition of a uniform strain in direction 1 is obtained by assigning a displacement equal to 1.6 mm to all nodes on the right vertical side and constraining the opposite side against movements in the same direction. The node at the bottom left comer 258 is constrained against vertical motion. The results are in agreement with the exact solution to the number of significant digits available, i.e. axial strain is uniform in direction 1 and equal to 0.1 in all elements of the patch. The displacement distribution in direction 1 is reported in Figure 20. Similarly, the condition of a uniform strain in direction 2 is obtained by assigning a displacement equal to 1 mm to all nodes on the top horizontal side and constraining the opposite side against movements in direction 2. The node at the bottom left comer is constrained against horizontal motion. The interface element results are in agreement with the exact solution to the number of significant digits available, i.e. axial strain is uniform in direction 2 and equal to 0.1 in all elements of the patch. The displacement distribution in direction 2 is reported in Figure 21. The condition of a uniform inplane shear strain is obtained by assigning displacements to the nodes on the boundaries such that horizontal sides form an angle with the axis 1 and vertical sides form an angle with the axis 2. The interface element results are in agreement with the exact solution to the number of significant digits available, i.e. shear strain is uniform in all elements of the patch. The displacement distributions in direction 1 and 2 are reported in Figures 22(3) and 22(b). Creating the input file The input file for the uniform axial strain loading condition is used to describe the command lines necessary to define the interface model. The complete input file is provided at the end of the section. In the following, all lines in the input file related to the interface element are reported in bold. For this problem five equally spaced nodes have been used to define the penalty 259 interface element. The coordinates of each of the five nodes composing the interface frame must be added to the node definition part of the input file. Note that the two subdomains do not have any common nodes. Note that the coordinates of the first and of the last nodes of the three groups of nodes must match. In the present case, the two groups of nodes (2, 101 , 20) and (4, 103, 28) must have the same coordinates. *HEADING Patch Test - Constant strain in direction 1 *NODE 1, .00000, .00000, .00000 2, 6.0000, .00000, .00000 3, 3.0000, .00000, .00000 4, 10.000, 10.000, .00000 5, 8.0000, 5.0000, .00000 6, .00000, 10.000, .00000 7, 5.0000, 10.000, .00000 8, .00000, 5.0000, .00000 9, 4.0000, 5.0000, .00000 20, 6.0000, .00000, .00000 21, 16.000, .00000, .00000 22, 11.000, .00000, .00000 23, 16.000, 10.000, .00000 24, 16.000, 2.0000, .00000 25, 16.000, 4.0000, .00000 26, 16.000, 6.0000, .00000 27, 16.000, 8.0000, .00000 28, 10.000, 10.000, .00000 29, 13.000, 10.000, .00000 30, 9.2000, 8.0000, .00000 31, 8.4000, 6.0000, .00000 32, 7.6000, 4.0000, .00000 33, 6.8000, 2.0000, .00000 34, 11.400, 2.0000, .00000 35, 11.800, 4.0000, .00000 36, 12.200, 6.0000, .00000 37, 12.600, 8.0000, .00000 101, 6.000, .00000, .00000 102, 7.000, 2.5000, .00000 103, 8.000, 5.0000, .00000 104, 9.000, 7.5000, .00000 105, 10.000, 10.000, .00000 260 i““*’_13 The element connectivity for each of the two subdomains does not require any input concerning the interface element. Note that the two subdomains are not connected. *ELEMENT, TYPE=CPS4, ELSET=BEAM 1, 1, 3, 9, 8 2, 3, 2, 5, 9 3, 8, 9, 7, 6 4, 9, 5, 4, 7 20, 20, 22, 34, 33 21, 22, 21, 24, 34 22, 33, 34, 35, 32 23, 34, 24, 25, 35 24, 32, 35, 36, 31 25, 35, 25, 26, 36 26, 31, 36, 37, 3O 27, 36, 26, 27, 37 28, 30, 37, 29, 28 29, 37, 27, 23, 29 The element section and material properties of the native ABAQUS elements are defined in the usual manner. *SOLID SECTION, ELSET=BEAM, MATERIAL=STEEL 1 *MATERIAL, NAME=STEEL *ELASTIC 1.E+O3, 0.3 The group of data related to the interface element begins with the *USER ELEMENT keyword: *USER ELEMENT, NODES=14, TYPE=U1, PROPERTIES=5, IPROPERTIES=3, COORDINATES=2, VARIABLES=4 1 , 2 These lines should be added to the .inp file exactly as they appear above, with the exception that the parameter NODE must be set equal to the total number of nodes at the interface. The total number of interface nodes is the sum of the number of interface nodes 261 from each of the two independent subdomains plus the nodes composing the interface element. The user should not change the parameters PROPERTIES, IPROPERTIES, COORDINATES and VARIABLES, since they are the same for every analysis. The numbers in the following line “1, 2” define which degrees of freedom are active for the interface element, and this line must not be changed. Note that if all of the interface elements in a mesh contain the same total number of nodes at the interface, then the above definition of *USER ELEMENT is used to define all such elements. For each interface element in the mesh that has a different number of total interface nodes, the *USER ELEMENT card should be repeated, with n incremented by one within the TYPE=Un parameter. For example, the second interface element definition will use TYPE=U2. The keyword *ELEMENT defines the element number and connectivity of the interface element: *ELEMENT, TYPE=U1, ELSET=FRAME1 100, 2, 5, 4, 101, 102, 103, 104, 105, 20, 33, 32, 31, 30, 28 The parameters TYPE and ELSET are required, where TYPE is set equal to the user element type (Un, where usually n=1) and ELSET identifies the set of interface elements sharing the same properties. The ELSET definition in the *ELEMENT card should match the ELSET definition in the *UEL PROPERTY card described below. The second line is required, as it contains the interface element number (100 in this case) followed by the list of nodes that compose the interface element. The list of nodes must 262 follow a particular order. Beginning with one of the two subdomain regions being connected, the nodes of this subdomain must be listed in order as they appear along the interface. It does not matter which “end” of this interface region contains the first node in the list, but the same direction chosen for the first subdomain interface region must be used for the subsequent definition of the interface element as well as for the second subdomain region. Next, the five equally spaced nodes of the interface element are defined in order according to the direction defined by the first subdomain interface region. Finally, the nodes of the second mesh are defined in order as they appear and in the same direction as defined in the first subdomain interface region. The keyword *UEL PROPERTY allows the specification of the main properties of the interface element: *UEL PROPERTY , ELSET=FRAME1 1.E-l-5, 1., 1., l.E+3, l.E+3, 3, 6, 5 ** BETA, 'rnru , THKZ , 21 , £2 , N, M, N93 The parameter ELSET is required and must match the definition in *ELEMENT. The second line contains the real and integer data for the interface element. The third line is an Optional comment line (indicated by **) that is used as an aid to remember the type of data required in the second line. The parameter BETA is a scaling factor for the penalty parameter. The value l.E+4 should work well for most analyses. Values of BETA smaller than l.E+3 often do not enforce compatibility sufficiently and are not recommended. Higher values are allowed thanks to presence of a round-off error control that automatically reduces the value of BETA if it is causing numerical problems. The meaning of the parameter BETA is clearly explained in the technical documentation on 263 my. the interface element. The next parameters are the thickness and Young’s modulus of the two subdomains, given in the same order as the nodes in *ELEMENT. In case of an orthotropic material, the Young’s modulus required is the one in the direction perpendicular to the interface. N and M are the number of nodes at the interface for each of the two subdomain meshes. NPS (Number Points Spline) is the number of equally spaced nodes at the interface used to define the penalty interface element. In the present example, the first interface has three nodes (2, 5, 4), the second interface has six nodes (20, 33, 32, 31, 30, 28) and five nodes (101, 102, 103, 104, 105) have been used to define the penalty interface element, so N=3, M=6 and NPS=5. Executing the program The Abaqus command line for performing an analysis with interface elements is the following: abaqus job=job_name user=ue1glue interactive The input file (job_name.inp), Abaqus include file aba _param. inc and the program in its compiled version uelglueo or as source code uelglue.f must be present in the directory where the analysis is performed. Example Input File *HEADING Patch Test - Constant strain in direction 1 *NODE 1, .00000, .00000, .00000 2, 6.0000, .00000, .00000 3, 3.0000, .00000, .00000 264 :3.“ --.-—--l 4, 10.000, 10.000, .00000 5, 8.0000, 5.0000, .00000 6, .00000, 10.000, .00000 7, 5.0000, 10.000, .00000 8, .00000, 5.0000, .00000 9, 4.0000, 5.0000, .00000 20, 6.0000, .00000, .00000 21, 16.000, .00000, .00000 22, 11.000, .00000, .00000 23, 16.000, 10.000, .00000 24, 16.000, 2.0000, .00000 25, 16.000, 4.0000, .00000 26, 16.000, 6.0000, .00000 27, 16.000, 8.0000, .00000 28, 10.000, 10.000, .00000 29, 13.000, 10.000, .00000 30, 9.2000, 8.0000, .00000 31, 8.4000, 6.0000, .00000 32, 7.6000, 4.0000, .00000 33, 6.8000, 2.0000, .00000 34, 11.400, 2.0000, .00000 35, 11.800, 4.0000, .00000 36, 12.200, 6.0000, .00000 37, 12.600, 8.0000, .00000 101, 6.000, .00000, .00000 102, 7.000, 2.5000, .00000 103, 8.000, 5.0000, .00000 104, 9.000, 7.5000, .00000 105, 10.000, 10.000, .00000 *ELEMENT, TYPE=CPS4, ELSET=BEAM 1, 1, 3, 9, 8 2, 3, 2, 5, 9 3, 8, 9, 7, 6 4, 9, 5, 4, 7 20, 20, 22, 34, 33 21, 22, 21, 24, 34 22, 33, 34, 35, 32 23, 34, 24, 25, 35 24, 32, 35, 36, 31 25, 35, 25, 26, 36 26, 31, 36, 37, 3O 27, 36, 26, 27, 37 28, 30, 37, 29, 28 29, 37, 27, 23, 29 *SOLID SECTION, ELSET=BEAM, MATERIAL=STEEL 1 *MATERIAL, NAME=STEEL 265 *ELASTIC l.E+3, 0.3 *USER ELEMENT , NODES=14 , TYPE=U1 , PROPERTIES=5 , IPROPERTIES=3 , COORDINATES=2 , VARIABLES=4 1, 2 *ELEMENT, TYPE=UI , ELSET=FRAMEI 100, 2, 5, 4, 101, 102, 103, 104, 105, 20, 33, 32, 31, 30, 28 *UEL PROPERTY, ELSET=FRAME1 1.E+5, 1., 1., l.E+3, l.E+3, 3, 6, 5 ** BETA, THKI, THKZ, E1, EZ, N, M, NPS *BOUNDARY 1, 1, 2, 6, 1, 1, 8, 1, 1, *STEP, PE *STATIC *BOUNDARY 21, 1, 23, 1 24, 1 25, 1 1 1 0.0 0.0 0.0 RTURBATION ‘ I ‘ ‘ ‘ l—‘l—‘l—‘l—‘l—J ‘ 26r I r 27, I 1, *NODE PRINT U, RF *RESTART, WRITE *END STEP ‘ H‘H‘HIJIAFA oxoxox0xoxm 266 REFERENCES 267 [1] [2] l3] [4] [5] [6] [8] REFERENCES Aminpour, MA. and K. A. Hosapple. "Finite element solutions for propagating interface crack with singularity elements", Engineering Fracture Mechanics, Vol. 39, No. 3, 1991, pp. 451-468. Z. Jinping and A. Huizu. "Stress analysis around holes in orthotropic plates by subregion mixed finite element method", Computers & Structures, Vol. 41, No. l, 1991, pp. 105-108. Farhat and F. X. Roux. "A method of finite element tearing and interconnecting and its parallel solution algorithm", Journal for Numerical Methods in Engineering, Vol. 32, 1991, pp. 1205-1227. C. Farhat and M. Gerardin. "Using a reduced number of Lagrange multipliers for assembling parallel incomplete field finite element approximations", Computer Methods in Applied Mechanics and Engineering, Vol. 97, 1992, pp. 333-3 54. Y. Maday, C. Mavriplis and A. Patera. "Non-conforming mortar element methods: Application to spectral discretizations", NASA CR-181729, ICASE Report No. 88-59, 1988. J .B. Ransom, S.L. McCleary and M. A. Aminpour. "A New Interface Element for Connecting Independently Modeled Substructures", AIAA Paper, Number 93- 1503, 1993. J .M. Housner, M. A. Aminpour, C.G. Davila, J .E. Schiermeier, W.F. Stroud, J.B. Ransom, and RE. Gillian. "An Interface Element for Global/Local and Substructuring Analysis", presented at the MSC World Users’ Conference, Los Angeles, CA, May 8-12, 1995. M. A. Aminpour, J. B. Ransom, and S. L. McCleary. "A Coupled Analysis Method for Structures with Independently Modeled Finite Element Subdomains", 268 [9] [10] [11] 112] [13] [14] 115] [161 117] International Journal for Numerical Methods in Engineering, Vol. 38, 1995, pp. 3695-3718. J.B. Ransom. “Interface Technology for Geometrically Nonlinear Analysis of Multiple Connected Subdomains”, AIAA Paper No. 97-1298, 1997. M. A. Aminpour and T. Krishnamurthy. "A Two-Dimensional Interface Element for Multi-Domain Analysis of Independently Modeled Three-Dimensional Finite Element Meshes", AIAA paper no. 97-1297, pp. 1853-1861, 1997. M. A. Aminpour, T. Krishnamurthy and T. D. Fadale. “Coupling of Independently Modeled Three-Dimensional Finite Element Meshes with Arbitrary Shape Interface Boundaries”, AIAA paper no. 98-2060, pp. 1853-1861, 1998. M. A. Aminpour, Stephane Pageau and Youngwon Shin. “An Alternative Method for the Interface Modeling Technology”, AIAA paper no. 2000-1352, pp. 1-13, 2000. Maenghyo Cho and Won Bae Kim. "A Coupled Finite Element Analysis of Independently Modeled Substructures by Penalty Frame Method", AIAA paper no. 98-2061, pp. 3025-3032, 1998. Carlos A. Felippa. "Error Analysis of Penalty Function Techniques for Constraint Definition in Linear Algebraic System", International Journal for Numerical Methods in Engineering, Vol. 11, pp. 709-728, 1977. Carlos A. Felippa. "Interactive Procedures for Improving Penalty Function Solutions of Algebraic System", International Journal for Numerical Methods in Engineering, Vol. 12, pp. 821-836, 1978. Carlos A. Felippa. "Penalty-Function Interactive Procedures for Mixed Finite Elements Formulation", International Journal for Numerical Methods in Engineering, Vol. 22, pp. 267-279, 1986. G. Prathap. "Locking, Rank and Singularity of Penalty-Linked Stiffness Matrix and Consistency of Strain-Field", Computers & Structures, Vol. 52, No. 1, pp. 35- 39, 1994. 269 EFT-”131m ~ I .4'" 0»: [18] [19] [20] [21] [22] [23] [24] [25] [261 [27] Ahmed K. Noor and Jeanne M. Peters. "Penalty Finite Element Models for Nonlinear Dynamic Analysis", AIAA Journal, Vol. 24, No. 2, pp. 312-320, 1986. Ahmed K. Noor, M. Asce and Jeanne M. Peters. "Penalty Finite Element Formulation for Curved Elastica", Journal of Eng. Mechanics, Vol. 110, No. 5, pp. 694-712, 1984. E. Haugeneder. "A New Penalty Function Element for Thin Shell Analysis", International Journal for Numerical Methods in Engineering, Vol. 18, pp. 845- 861, 1982. Ted Belytschko and Mark 0. Neal. "Contact-Irnpact by the Pinball Algorithm with Penalty and Lagrangian Methods", International Journal for Numerical Methods in Engineering, Vol. 31, pp. 547-572, 1991. G. F. Carey and R. Krishnan. "Convergence of Iterative Methods In Penalty Finite Element Approximation Of The Navier-Stokes Equations", Computer Methods In Applied Mechanics And Engineering, Vol. 60, pp. 1-29, 1991. J. N. Reddy. "On The Accuracy and Existence of Solutions to Primitive Variable Models of Viscous Incompressible Fluids", Int. J. Engineering Sci, Vol. 16, pp. 921-929, 1978. Henry Bertin and Hiroyuki Ozoe. "Technique for Rapid Convergence of the Penalty Finite-Element Method with a Modified Galerkin Scheme and Its Application to Natural Convection", Numerical Heat Transfer, vol. 10, pp. 311- 325, 1986. E-M. Salonen. "An Iterative Penalty Function Method in Structural Analysis", International Journal for Numerical Methods in Engineering, Vol. 10, pp. 413- 421, 1976. R. Codina, M. Cervera And E. Orate. "A Penalty Finite Element Method for Non- Newtonian Creeping Flows", International Journal for Numerical Methods in Engineering, Vol. 36, pp. 1395-1412, 1993. A. W. Al-Khafaji and J. R. Tooley. "Numerical Methods in Engineering Practice", HRW Series in Mechanical Engineering, 1986. 270 [23] [29] [30] [31] [32] [33] [34] I35] [36] [37] [381 [39] [40] M. L. James, G.M. Smith and J. C. Wolford. "Applied Numerical Methods for Digital Computation", Harper & Row, 1985. Arthur Sard and S01 Weintraub. "A Book of Splines", Wiley, 1971. Larry L. Schumaker, “Spline functions: basic theory”. New York: Wiley, c1981. P. M. Prenter. “Splines and variational methods”, New York: Wiley, 1975. Elaine Cohen, Richard F. Riesenfeld and Gershon Elber. “Geometric modeling with splines: an introduction”. Natick, Mass: AK Peters, c2001. C. de Boor, K. Hollig, S. Riemenschneider. “Box splines”. New York: Springer- Verlag, c1993. Pierre-Jean Laurent, Alain Le Méhauté, Larry L. Schumaker. “Curves and surfaces”. Boston: Academic Press, c1991. J. H. Ahlberg, E. N. Nilson. “The theory of splines and their applications”. New York, Academic Press, 1967. ‘ Gerald E. Farin. “Curves and surfaces for computer aided geometric design: a practical guide”, Boston: Academic Press, 1993. Paul Dierckx. “Curves and surfaces fitting with splines”. New York: Clarendon Press, 1993. H Spath. “Two dimensional spline interpolation algorithms”, Wellesley, Mass: AK Peters, 1995. Abaqus Manuals, Hibbit, Karlsson & Sorensen Inc., 1999. Jonathan B. Ransom. “On Multifunctional Collaborative Methods in Engineering Science”, Langley Research Center, Hampton, Virginia, NASA/TM-2001-211046, 2001. 271 [41] [42] [43] [44] [45] [46] [47] [48] J. N. Reddy. “An introduction to the finite element method”. Highstown, NJ, McGraw-Hill, 1993. R. D. Cook, D. S. Malkus, M. E. Plesha. “Concepts and applications of finite element analysis”. New York: Wiley, 1989. K-J Bathe. “Finite element procedures”. Englewood Cliffs, NJ, Prentice Hall, 1996. W.Chen and S. Yang. "Multilayered Hybrid-Stress Finite Element Analysis of Composite Laminates with Delamination Crack Originating from Transverse Cracking", Engineering Fracture Mechanics, Vol. 54, No. 5, pp. 713-729, 1996. A. A. Griffith. "The Phenomena of Rupture and Flow in Solids", Phil. Trans, A221, 1921, pp. 163-198. G. R. Irwin "Analysis of stresses and strains near the end of a crack traversing a plate", J. of Applied Mechanics, Transaction of ASME, pp. 361-364, 1957. H. Tada, P. C. Paris and G. R. Irwin. "The stress analysis of cracks handbook", ASME Press — Third Edition, New York, NY, USA, 2000. T. Ireman, J .C. Thesken, E. Greenhalgh, R. Sharp, M. Gadke, S. Maison, Y. Ousset, F. Roudolff, A. La Barbera. "Damage propagation in composite structural elements - coupon experiments and analyses", Composite Structures, v 36, n 3-4, Nov-Dec, p 209-220, 1996. D. Brock. "The practical use of Fracture Mechanics", Kluwer Academic Publishers, Boston, MA, USA, 1989. A. R. Luxmoore and D. R. Owen "Numerical Methods in Fracture Mechanics", Proceedings of the First International Conference — Swansea, UK, 1978. F. Ozdil and L. A. Carlsson. "Characterization of mode I delamination growth in glass/epoxy composite cylinders", J. of Composite Materials, Vol. 34, No. 5, pp. 398-419, 2000. 272 [52] [53] [54] [55] [56] 157] [58] [59] [60] [61] N. S. Choi, A. J. Kinloch and J. G. Williams. "Delamination fracture of multidirectional carbon-fiber/epoxy composites under mode I, mode 11 and mixed-mode I/II loading", J. of Composite Materials, Vol. 33, No. 1, pp. 73-100, 1999. M. Kenane, M.L. Benzeggagh. "Mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites under fatigue loading", Composites Science and Technology, v 57, n 5, May, p 597-605, 1997. F. Ducept, P. Davies, D. Gamby. "An experimental study to validate tests used to determine mixed mode failure criteria of glass/epoxy composites", Composites - Part A: Applied Science and Manufacturing, v 28A, n 8, p 719-729, 1997. Chyanbin ku, C.J. Kao, L.E. Chang. "Delamination fracture criteria for composite laminates", Journal of Composite Materials, v29, n15, p 1962-1987, 1995. A.J. Kinloch, Y. Wang, J.G. Williams, P. Yayla. "Mixed-mode delamination of fiber composite materials", Composites Science and Technology, v 47, n 3, p 225- 237, 1993. N. J. Pagano. "Interlaminar response of composite materials", Composite Materials Series, Elsevier, Editor: N. J. Pagano, Vol. 5, 1989. A. C. Garg. "Delamination — A Damage Mode in Composite Structures", Engineering Fracture Mechanics, Vol. 29, No. 5, pp. 557-584, 1988. SC. Pradhan and T. E. Tay. "Three-Dimensional Finite Element Modelling of Delamination Growth in Notched Composite Laminates Under Compressive Loading", Engineering Fracture Mechanics, Vol. 60, No. 2, pp. 157-171, 1998. Chyanbin ku, C.J. Kao, L.E. Chang. "Delamination fracture criteria for composite laminates", Journal of Composite Materials, v29, n15, p 1962-1987, 1995. K. Kaczmarek, M.R. Wisnom, M.I. Jones. "Edge delamination in curved (O/i45)s glass-fibre/epoxy beams loaded in bending", Composites Science and Technology, v 58, n 1, Jan, p155-161,1998. 273 [62] [63] [64] [65] [66] [67] [68] [69] [70] ED. Jr. Reedy, F .J. Mello, T.R. Guess. "Modeling the initiation and grth of delaminations in composite structures", Journal of Composite Materials, v 31, n 8, p 812-831,1997. S. Rinderknecht, B. Kroplin. "A computational method for the analysis of delamination growth in composite plates", Computers and Structures, v 64, n 1-4, Jul-Aug, p 359-374, 1997. J. T. Wang and I. S. Raju. "Strain Energy Release Rate Formulae for Skin- Stiffener Debond Modeled with Plate Elements", Engineering Fracture Mechanics, Vol. 54, No. 2, pp. 221-228, 1996. LS. Raju, R. Sistla and T. Krishnamurthy. "Fracture Mechanics Analyses for Skin-Stiffener Debonding", Engineering Fracture Mechanics, Vol. 54, No. 3, pp. 371-385, 1996. D. Hitchings, P. Robinson and F. Javidrad. "Finite Element Model for Delamination Propagation Composites", Computers & Structures, Vol. 60, No. 6, pp. 1093-1104, 1996. C. Davila, P. P. Camanho and M. F. de Moura. "Mixed-Mode Decohesion Elements for Analyses of Progressive Delamination", 42"“ AIAA/ASME/ASCE/A HS/ASC Structural Dynamics and Materials Conference, Seattle, Washington, paper AIAA-01-l486, 2001. Y. Mi, A. Crisfield, A. 0. Davies and H. B. Hellweg. "Progressive Delamination Using Interface Elements", .1. of Composite Materials, Vol. 32, No. 14, pp. 1246- 1272, 1998. J. Chen, M. Crisfield, A. J. Kinloch, E. P. Busso, F. L. Mattehws and Y. Qiu. "Predicting Progressive Delamination of Composite Materials Specimens via Interface Elements", Mechanics of Composite Materials and Structures, Vol. 6, pp. 301-317, 1999. G. Alfano and M. A. Crisfield. "Finite element interface models for the delamination analysis of laminated composites: mechanical and computational issies", International Journal for Numerical Methods in Engineering, Vol. 50, pp. 1701-1736, 2001. 274 [71] [72] [73] [741 [75] [76] I77] [78] [79] [30] L. Lammerant, I. Verpoest. "Modelling of the interaction between matrix cracks and delaminations during impact of composite plates", Composites Science and Technology, v 56, n 10, Oct, p 1171-1178, 1996. J. C. J. Schellekens and R. DE Boerst. "A Non-Linear Finite Element Approach for the Analysis of Mode-I Free Edge Delamination in Composites", Int. J. of Solids Structures, Vol. 30, No. 9, pp. 1939-1953, 1993. J. H. A. Schipperen and F.J. Lingen. "Validation of two-dimensional calculation of free-edge dalamination in laminated composites", Composite Structures, Vol. 45, pp. 233-240, 1999. C. M. Dakshima Moorthy and J. N. Reddy. "Recovery of Interlaminar Stresses and Strain Energy Release Rates in Composite Laminates", Finite Elements in Analysis and Design, Vol. 33, pp. 1-27, 1999. O. Allix and A. Corigliano. "Geometrical and interfacial non-linearities in the analysis of delamination in composites", 1. J. of Solids and Structures, Vol. 36, pp. 2189-2216, 1999. N. Point and E. Sacco. "A Delamination Model for Laminated Composites", Int. J. of Solids Structures, Vol. 33, No. 4, pp. 483-509, 1996. P. Ladeveze. "A Damage Computational Method for Composite Structures", Computers & Structures, Vol. 44, No. 1/2, pp. 78-87, 1992. W. J. Bottega. "A growth law for propagation of arbitrary shaped delaminations in layered plates", Int. J. of Solids Structures, Vol. 19, No. l 1, pp. 1009-1017, 1983. Z. Petrossian, Michael R. Wisnom. "Prediction of delamination initiation and growth from discontinuous plies using interface elements", Composites - Part A: Applied Science and Manufacturing, v 29, n 5-6, p 503-515, 1998. D. Chakraborty and B. Pradhan. "Effect of ply thickness and fibre orientation on delamination initiation in broken ply composite laminate", J. of Reinforced Plastic and Composites, Vol. 18, No. 8, pp. 735-758, 1999. 275 [81] [82] [83] [84] [85] [86] I87] [88] [89] [90] J. W. Joo and C. T. Sun. “Failure criterion for laminates governed by free edge interlaminar shear stress”, Journal of Composite Materials, v. 26, n. 10, pp 573- 586, 1994. Chu-Cheng Ko, Chien-Chang Lin, Hsiang Chin. “Prediction for delamination initiation around holes in symmetric laminates”, Composite Structures, v 22, n 4, p 187-191, 1992. S. Mohammadi, D. R. J. Owen and D. Peric. "A Combined Finite/Discrete Element Algorithm for Delamination Analysis of Composites", Finite Elements in Analysis and Design, Vol. 28, pp. 321-336, 1998. J. Zhao, S. V. Hoa, X. R. Xiao and 1. Hanna. "Global/Local Approach Using Partial Hybrid Finite Element Analysis of Stress Fields in Laminated Composites with Delamination Under Bending", J. of Reinforced Plastic and Composites, Vol. 18, No. 9, pp. 827-843, 1999. F. Chang and G. S. Springer. "Strength of fiber reinforced composite bends", Journal of Composite Materials, v. 20, pp 30-45, 1986. M Ortiz and A. Pandolfi. “Finite-deformation irreversible cohesive element for three-dimensional crack-propagation analysis”. International Journal for Numerical Methods in Engineering, Vol. 44, pp. 1267-1282, 1999. A. Pandolfi, P. Krysl and M Ortiz and. “Finite element simulation of ring expansion and fragmentation: The capturing of length and time scale through cohesive models of fracture”. International Journal of Fracture, Vol. 95, pp. 279- 297, 1999. N. Moes, J. Dolbow and T. Belytschko. "A finite element method for crack grth without remeshing", Int. J. Numer. Meth. Engng., Vol. 46, pp. 131-150, 1999. MA. Aminpour and K. A. Hosapple. "Finite element solutions for propagating interface crack with singularity elements", Engineering Fracture Mechanics, Vol. 39, No. 3, 1991, pp. 451-468. Jin Lee and Huajian Gao. “A hybrid finite element analysis of interface cracks”. International Journal for Numerical Methods in Engineering, Vol. 38, pp. 2465- 2482, 1995. 276 [91] [92] [93] [94] [95] [96] [97] [98] I99] [100] [101] AS. Krausz and K. Krausz. "Fracture kinetics and crack growth", Kluwer Academic Publishers, Boston, MA, USA, 1988. M. H. Aliabadi and D. P. Rooke. "Numerical Fracture Mechanics", Kluwer Academic Publisher, 1991 . R. El Abdi. “A special finite element for analysis of the singularity in a bimaterial containing a crack perpendicular to the interface”. Engineering Fracture Mechanics, Vol. 39, No. 6, pp. 1061-1065, 1991. S. H. Staab and T. C. Chang. “A finite element analysis of the finite width interface mixed mode bimaterial fracture problem", Computers & Structures, Vol. 19, No. 5/6, pp. 879-884, 1984. K. Saito, S. Araki, T. Kawakami and I. Moriwaky. “Global-local finite element analysis of stress intensity factor for a crack along the interface of a two phase material", Proceedings of the [CCM-10, Whistler, BC, Canada, Vol. I, August, 1995. Q. Yang and Q. Qin. “Numerical simulation of cracking process in dissimilar media", Composite Structures, Vol. 53, pp. 403-407, 2001. K. Machida. “Stress intensity factor of three-dimensional interface crack under mixed-mode loading by finite element analysis”. Key Engineering Materials, Vol. 145-149, pp. 589-594, 1998. F. Ozdil and L. A. Carlsson. "Characterization of mode 11 delamination growth in glass/epoxy composite cylinders", .l. of Composite Materials, Vol. 34, No. 4, pp. 274-298, 2000. F. Ozdil and L. A. Carlsson. "Characterization of mixed mode delamination in glass/epoxy composite cylinders", J. of Composite Materials, Vol. 34, No. 5, pp. 420-441, 2000. Aditi Chattopadhyay, Changho Nam and Dan Dragomir-Daescu. "Delamination modeling and detection in smart composite plates", J. of Reinforced Plastic and Composites, Vol. 18, No. 17, pp. 1557-1572, 1999. R. M. Jones. "Mechanics of composite materials", Taylor & Francis Inc., Philadelphia, PA, USA, 1999. 277 [102] [103] [104] [105] [106] [107] [108] [109] [110] [111] [112] A. S. D. Wang. "Free-edge delamination in layered composites- A historical review of the past 30 years", ASMEI998, MD-Vol. 84, pp. 223-226, 1998. S. H. Narayan and J. L. Beuth. "Designation of mode mix in orthotropic composite delamination problems", Int. J. of Fracture, Vol. 90, pp. 383-400, 1998. Narendra V. Bhat, Paul A. Lagace. "Analytical method for the evaluation of interlaminar stresses due to material discontinuities", Journal of Composite Materials,v28,n3,p190-210,1994. Wan-Lee Yin. "Simple solution of the free-edge stresses in composite laminates under thermal and mechanical loads", Journal of Composite Materials, v. 28, n. 6, pp 573-586, 1994. Jerry Zhiqi Wang, Darrell F. Socie. "Failure strength and damage mechanisms of e-glass/epoxy laminates under in-plane biaxial compressive deformation", Journal of Composite Materials, v 27, n 1, pp 40-58, 1993. S. F. Muller de Almeida and G. M. Candido. "Effect of the free edge finishing on the tensile strength of carbon/epoxy laminates", Composite Structures, v 25, pp 287-293, 1993. W. Becher. "Closed-form solution for the free-edge effect in cross-ply laminates", Composites Structures, v 26, pp 39-45, 1993. D. Bruno and A. Grimaldi. "Delamination Failure of Layered Composite Plates Loaded in Compression", Int. J. of Solids Structures, Vol. 26, No. 3, pp. 313-330, 1990. W. J. Bottega. "A growth law for propagation of arbitrary shaped delaminations in layered plates", Int. J. of Solids Structures, Vol. 19, No. 11, pp. 1009-1017, 1983. W. J. Bottega and A. Maewal. "Delaminations buckling and growth in laminates", J. of Applied Mechanics, Vol. 50, pp. 184-189, 1983. J. D. Whitcomb and I. S. Raju. "Superposition Method for Analysis of F rec-Edge Stresses", J. of Composite Materials, Vol. 17, pp. 492-507, 1983. 278 [113] [114] [115] [116] [117] [118] [119] [120} 1121] [122] [123] S. S. Wang and I. Choi. "Boundary-layer effects in composite laminates: Part 1 — F rec-edge singularities", J. of Applied Mechanics, Vol. 49, pp. 541-548, 1982. S. S. Wang and I. Choi. "Boundary-layer effects in composite laminates: Part 2 — Free-edge stress solution and basic characteristics", J. of Applied Mechanics, Vol. 49, pp. 549-560, 1982. E. Altus, A. Rotem and M. Shmueli. "Free-Edge Effect in Angle Ply Laminates — A New Three Dimensional Finite Difference Solution", J. of Composite Materials, Vol. 14, pp. 21-30, 1980. R. L. Spliker and S. C. Chou. "Edge Effects in symmetric composite Laminates: importance of satisfying the traction-free-edge condition", J. of Composite Materials, Vol. 14, pp. 2-20, 1980. N. J. Pagano and E. F. Rybicki. "On the significance of effective modulus solutions for fibrous composites", .1. of Composite Materials, Vol. 8, pp. 214-228, 1974. R. Byron Pipes and N. J. Pagano. "The influence of stacking sequence on laminate strength ", .1. of Composite Materials, Vol. 5, pp. 50-57, 1971. R. Byron Pipes and N. J. Pagano. "Interlaminar stresses in composite laminates under uniform axial extension", J. of Composite Materials, Vol. 4, pp. 538-548, 1970. H. S. Kim, M. Cho and G. Kim. "Free-edge strength analysis in composite laminates by the extended Kantorovic method", Composite Structures, Vol. 49, pp. 229-235, 2000. M. Schulze and W. D. Nix. "Finite Element Analysis of the Wedge Delamination Test ", Int. J. of Solids and Structures, Vol. 37, pp. 1045-1063, 2000. Y. Ousset. "Numerical simulation of delamination growth in layered composite plates", 1. J. of Mechanics and Solids, Vol. 18, pp. 291-312, 1999. J. Macklerle. "Delamination of Composites - Finite Element and Boundary Element Analyses -— A bibliography (1995-1997)", Finite Elements in Analysis and Design, Vol. 30, pp. 243-252, 1998. 279 [124] [125] [126] [127] [128] [129] H. T. J. Yang and C. C. He. "Three-Dimensional Finite Element Analysis of Free Edge Stresses and Delamination Composite Laminates", J. of Composite Materials, Vol. 28, No. 15, pp. 1394-1412, 1994. K. N. Shivakumar and I. S. Raju. "An Equivalent Domain Integral Method for Three-Dimensional Mixed-Mode Fracture Problems", Engineering Fracture Mechanics, Vol. 42, No. 6, pp. 935-959, 1992. X. Lu and D. Liu. "Finite Element Analysis of Strain Energy Release Rate at Delamination Front", J. of Reinforced Plastic and Composites, Vol. 10, pp. 279- 292, 1991. E. J. Barbero and J. N. Reddy. "Modelling of Delamination in Composite Laminates Using a Layer-Wise Plate Theory", Int. J. of Solids Structures, Vol. 28, No. 3, pp. 373-388, 1991. S. El-Sayed and S. Sridharan. “Predicting and tracking interlaminar crack grth in composites using a cohesive layer model”. Composites Part B: engineering, Vol. 32, pp. 545-553, 2001. N. E. Jansson and R. Larsson. “A damage modelfor simulation of mixed mode delamination growth”. Composites Structures, Vol. 53, pp. 409-417, 2001. 280 IIIIIIIIIIIIIIIIIIIIIIIII 11 11,1111 '1 L206 3 1293